A sparse direct solution method based on data out-of-core storage

Abstract

The present application relates to the technical field of data processing, and especially relates to a sparse direct solution method based on data kernel external storage, which comprises the following steps: obtaining a coefficient matrix of a high-dimension sparse linear equation group, analyzing a non-zero element distribution mode and a matrix topological structure to establish a task dependency relationship; detecting a system memory capacity and a disk performance parameter to set a slice size parameter, dividing a wavefront matrix into data slices with balanced floating-point operation amounts to execute a left-looking numerical decomposition; creating a message queue and a monitoring mechanism to realize asynchronous disk storage; reading data from the disk to execute a numerical update cycle; and establishing a back substitution solution memory work area to generate a numerical solution by using a previous generation back substitution algorithm. The present application utilizes disk storage to relieve memory pressure and supports large-scale finite element simulation application.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to a sparse direct solution method based on data kernel external storage. Background Technology

[0002] Existing sparse direct solution techniques based on off-core storage are specifically numerical methods for solving large-scale sparse linear equation systems. Their design is motivated by the challenge of insufficient memory capacity due to the enormous matrix dimensions in real-world problems. The core technology involves storing the sparse matrix and its padding elements generated during decomposition on an external disk. A block-based strategy divides the matrix into multiple sub-blocks, loading only the necessary sub-blocks into memory for processing, while retaining the remaining data on disk. The solution process optimizes the access order based on the matrix structure, employs efficient input / output management algorithms to reduce the number of data exchanges between disk and memory, and combines padding reduction preprocessing techniques to control computational complexity, thereby achieving stable solutions to large-scale problems in limited memory environments.

[0003] Existing sparse direct solution techniques based on off-core memory suffer from the following technical pain points: In industrial CAE simulations, when using the finite element method to analyze large-scale complex models, the numerical decomposition process of sparse matrices generates tens of times more data than the original matrix due to non-zero element filling, causing computational memory requirements to far exceed the physical memory capacity of a personal computer. Taking the simulation of automotive chassis structural strength as an example, after the model mesh is refined, equations with millions of degrees of freedom are generated. When using existing off-core solvers for numerical decomposition, the maximum wavefront matrix requires 20% of the decomposition factor memory, making it impossible for a 64GB memory platform to complete the solution task. This forces companies to purchase expensive servers or interrupt the simulation process, severely restricting product development efficiency and design iteration speed. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a sparse direct solution method based on data kernel external storage. This invention solves the technical problem that the surge in memory demand caused by non-zero element filling in sparse matrix numerical decomposition leads to insufficient computer physical memory, thus restricting large-scale finite element simulation applications.

[0005] To solve the above-mentioned technical problems, the specific contents of the present invention are as follows:

[0006] This invention provides a sparse direct solution method based on data kernel-external storage, comprising:

[0007] Step 1: Obtain the coefficient matrix of the high-dimensional sparse linear equation system to be solved, analyze the non-zero element distribution pattern and matrix topology of the coefficient matrix, establish the task dependency relationship of the multi-wavefront numerical decomposition, and generate the wavefront matrix corresponding to the task dependency relationship.

[0008] Step 2: Detect system memory capacity and disk performance parameters, calculate and set the slice size parameters based on memory capacity and disk performance parameters, divide each wavefront matrix into multiple data slices with balanced floating-point operation according to task dependencies and slice size parameters, perform left-look numerical decomposition calculation on each data slice, and generate decomposition factor data for the data slice.

[0009] Step 3: Create a message queue to be input / output and an input / output completion status monitoring mechanism. Store the decomposition factor data into the message queue to be input / output, extract the decomposition factor data from the message queue to be input / output and perform asynchronous disk storage. Release the corresponding memory buffer after the disk storage is completed by the input / output completion status monitoring mechanism.

[0010] Step 4: Read the required decomposition factor data from the disk, perform subsequent wavefront matrix numerical updates according to task dependencies, and repeat the process from patch division to numerical update. Monitor the loop execution status through an input / output completion status monitoring mechanism until all wavefront matrix decompositions are completed.

[0011] Step 5: Establish a back-substitution solution memory workspace. Read the complete stored decomposition factor data from the disk and load it into the back-substitution solution memory workspace. Perform numerical calculations on the decomposition factor data using the previous-back substitution algorithm in the back-substitution solution memory workspace. Combine the calculation results of each step to generate a numerical solution for a high-dimensional sparse linear system of equations.

[0012] Furthermore, in the sparse direct solution method based on data kernel external storage of the present invention, step 2 includes:

[0013] Obtain the dimension information and unprocessed principal component index number of the current wavefront matrix from the task dependencies;

[0014] The range of the principal index of the current region is calculated using the set region size parameters and the unprocessed principal index number.

[0015] The wavefront matrix is ​​divided into multiple data regions based on the range of the principal index.

[0016] Calculate the floating-point computation cost for each data segment and perform a balance adjustment.

[0017] Furthermore, the sparse direct solution method based on off-core storage of the present invention calculates the primary index range of the current region from the unprocessed primary index number obtained from the task dependency relationship and the region size parameter based on the system memory capacity and disk performance parameter settings, including:

[0018] Obtain the unprocessed principal index number q and the total number of principal elements N from the task dependencies of the current wavefront matrix;

[0019] Read the area size parameter S calculated based on system memory capacity and disk performance parameters;

[0020] The terminating principal index p of the current region is determined by the formula p = q + min( floor(S), N - q ), where floor represents the floor function.

[0021] The primary index range [q, p] is used as the calculation range for the current region.

[0022] Furthermore, in the sparse direct solution method based on data kernel external storage of the present invention, step 3 includes:

[0023] Create a producer-consumer model, which includes a producer role computation thread and a consumer role input / output thread;

[0024] Initialize the pending input / output message queue and the completed message queue. The two queues adopt a first-in-first-out (FIFO) management strategy.

[0025] The producer role calculation thread adds the decomposition factor data to the tail of the input / output message queue;

[0026] The consumer role's input / output thread retrieves the decomposition factor data from the head of the pending input / output message queue;

[0027] The consumer role's input / output thread performs disk storage of the decomposition factor data;

[0028] The consumer role's input / output thread adds disk write completion status information to the completed message queue.

[0029] Furthermore, in the sparse direct solution method based on data kernel external storage of the present invention, step 3 further includes:

[0030] Atomic variables are created and used to implement synchronous control for adding and deleting messages in the pending input / output message queue and the completed message queue;

[0031] Create and initialize global state variables, and use global state variables to pass input and output operation status signals between the producer role calculation thread and the consumer role input and output thread;

[0032] When the producer's computation thread and the consumer's input / output thread access the message queue, a spin-wait mechanism is used to ensure thread-safe queue access.

[0033] Furthermore, in the sparse direct solution method based on data kernel external storage of the present invention, step 4 includes:

[0034] Based on the task dependencies established in step 1, determine the subsequent wavefront matrix to be updated;

[0035] Read the decomposition factor data stored in step 3 from the disk;

[0036] Create a numerical update memory buffer and load the read decomposition factor data into the numerical update memory buffer;

[0037] A multi-wavefront numerical update algorithm is adopted, which uses the loaded decomposition factor data to perform numerical update calculations on the subsequent wavefront matrix.

[0038] Furthermore, in the sparse direct solution method based on data kernel external storage of the present invention, step 2 further includes:

[0039] The overall memory requirements for numerical decomposition are estimated based on the dimension of the coefficient matrix and the non-zero element distribution pattern.

[0040] Based on the estimated total memory requirements and the detected available physical memory capacity, the adjustment value of the region size parameter is calculated through an optimization algorithm;

[0041] Based on the detected disk input / output performance characteristics, an optimized value for the asynchronous storage memory buffer size is calculated using a performance optimization algorithm.

[0042] Furthermore, in the sparse direct solution method based on data kernel external storage of the present invention, step 2 further includes:

[0043] Before performing the left-look numerical decomposition calculation in step 2, the mathematical properties of the coefficient matrix are detected by a matrix analysis algorithm;

[0044] When a symmetric positive definite matrix is ​​detected, the Cholesky decomposition algorithm is used to perform numerical decomposition.

[0045] When an asymmetric matrix is ​​detected, the LU decomposition algorithm is used to perform numerical decomposition.

[0046] During the numerical decomposition calculation process, both real and complex number data types are processed simultaneously.

[0047] Furthermore, in the sparse direct solution method based on data kernel external storage of the present invention, step 4, monitoring the loop execution state through an input-output state monitoring mechanism, includes:

[0048] The status monitoring mechanism collects disk storage operation status signals through input and output;

[0049] Based on the monitored input and output operation status signals, a dynamic scheduling algorithm is used to dynamically schedule the execution tasks of the producer role computing threads.

[0050] The scheduling of computational tasks and disk input / output tasks is carried out in parallel, completing the overlap of computation and input / output operations.

[0051] Furthermore, in the sparse direct solution method based on off-core storage of the present invention, step 5, performing numerical calculations on the factorization data in the back-substitution solution memory working area using a previous-back substitution algorithm, includes:

[0052] Based on the solution order determined by the task dependencies, read the complete stored decomposition factor data from the disk in reverse order;

[0053] Load the read decomposition factor data into the back-substitution solution memory work area;

[0054] According to the calculation order required by the previous generation and back generation algorithm, the previous generation calculation and back generation calculation steps of the previous generation and back generation algorithm are executed in sequence.

[0055] By merging the results of previous and back-substitution calculations, a numerical solution to a high-dimensional sparse linear system of equations is generated.

[0056] Beneficial effects of this invention;

[0057] This invention effectively solves the technical problem of soaring memory requirements caused by non-zero element filling in sparse matrix numerical decomposition by using a sparse direct solution method based on data kernel storage. It utilizes disk space to store large-scale data, loading only the data currently necessary for computation into memory, significantly reducing peak memory usage. It establishes task dependencies through multi-wavefront numerical decomposition to guide the computation order, and sets fragment size parameters based on system memory capacity and disk performance parameters to achieve balanced floating-point computation by partitioning data into segments. Combined with a producer-consumer model, it implements asynchronous disk storage and an input / output completion status monitoring mechanism, allowing computation and storage operations to overlap and improve overall efficiency. It employs adaptive algorithms, such as Cholesky decomposition or LU decomposition, to optimize numerical stability based on matrix mathematical properties. It iteratively executes numerical updates and back-substitution solutions to ultimately generate numerical solutions for high-dimensional sparse linear equations, thereby supporting large-scale finite element simulation applications in a limited memory environment, reducing hardware costs and improving computational performance. Attached Figure Description

[0058] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments will be briefly introduced below. Obviously, for those skilled in the art, other drawings can be obtained based on the drawings without creative effort.

[0059] Figure 1 The flowchart illustrates a sparse direct solution method based on external storage of data kernels, as provided in an embodiment of the present invention.

[0060] Figure 2 The wavefront matrix based on region left-view decomposition process provided in the embodiments of the present invention.

[0061] Figure 3An asynchronous I / O process framework based on a producer-consumer model is provided for embodiments of the present invention. Detailed Implementation

[0062] To make the technical solution of the present invention clearer, the present invention will be clearly and completely described below with reference to specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. The present invention provided by various embodiments will be described in detail below with reference to the accompanying drawings. To better understand the purpose of the present invention, the present invention will be described in further detail below.

[0063] Please see Figure 1 The present invention provides a sparse direct solution method based on data kernel-external storage, comprising:

[0064] Step 1: Obtain the coefficient matrix of the high-dimensional sparse linear equation system to be solved, analyze the non-zero element distribution pattern and matrix topology of the coefficient matrix, establish the task dependency relationship of the multi-wavefront numerical decomposition, and generate the wavefront matrix corresponding to the task dependency relationship.

[0065] Step 2: Detect system memory capacity and disk performance parameters, calculate and set the slice size parameters based on memory capacity and disk performance parameters, divide each wavefront matrix into multiple data slices with balanced floating-point operation according to task dependencies and slice size parameters, perform left-look numerical decomposition calculation on each data slice, and generate decomposition factor data for the data slice.

[0066] Step 3: Create a message queue to be input / output and an input / output completion status monitoring mechanism. Store the decomposition factor data into the message queue to be input / output, extract the decomposition factor data from the message queue to be input / output and perform asynchronous disk storage. Release the corresponding memory buffer after the disk storage is completed by the input / output completion status monitoring mechanism.

[0067] Step 4: Read the required decomposition factor data from the disk, perform subsequent wavefront matrix numerical updates according to task dependencies, and repeat the process from patch division to numerical update. Monitor the loop execution status through an input / output completion status monitoring mechanism until all wavefront matrix decompositions are completed.

[0068] Step 5: Establish a back-substitution solution memory workspace. Read the complete stored decomposition factor data from the disk and load it into the back-substitution solution memory workspace. Perform numerical calculations on the decomposition factor data using the previous-back substitution algorithm in the back-substitution solution memory workspace. Combine the calculation results of each step to generate a numerical solution for a high-dimensional sparse linear system of equations.

[0069] This invention provides a sparse direct solution method based on off-core data storage, which achieves efficient solution of high-dimensional sparse linear equation systems through five interrelated technical steps. First, the coefficient matrix of the equation system to be solved is obtained. By analyzing the non-zero element distribution pattern and matrix topology, the task dependencies required for multi-wavefront numerical decomposition are established. These task dependencies determine the computational order and data dependencies between the various wavefront matrices, providing a logical foundation for subsequent calculations.

[0070] After establishing task dependencies, the system memory capacity and disk performance parameters are checked, and the slice size parameters are calculated and set based on the results. According to the task dependencies and slice size parameters, each wavefront matrix is ​​divided into multiple data slices with balanced floating-point computational load. Left-look numerical decomposition is performed on each data slice to generate decomposition factor data. Slicing ensures balanced computational load across all data slices, improving computational efficiency.

[0071] A pending input / output message queue and an input / output completion status monitoring mechanism are created. The generated factorization data is stored in the pending input / output message queue. Factorization data is retrieved from the pending input / output message queue and asynchronous disk storage operations are performed. The disk storage completion status is monitored through the input / output completion status monitoring mechanism. Upon detection of disk storage completion, the corresponding memory buffer is released, achieving dynamic management of memory resources.

[0072] The system reads the required decomposition factor data from disk and performs subsequent wavefront matrix numerical updates based on task dependencies. It iteratively executes the steps from partitioning to numerical updates, monitoring the loop execution status through an input / output completion state monitoring mechanism until the decomposition of all wavefront matrices is complete. This iterative execution mechanism enables large-scale computation.

[0073] A dedicated back-substitution solution memory workspace is established. Complete decomposition factor data stored on disk is read and loaded into the back-substitution solution memory workspace. Within this workspace, a previous-back-substitution algorithm is used to perform numerical calculations on the decomposition factor data. By combining the results of each step, a numerical solution to the high-dimensional sparse linear equation system is finally generated. The back-substitution solution process utilizes the complete decomposition data stored on disk, ensuring the accuracy and reliability of the calculation results.

[0074] The entire technical solution effectively addresses the memory capacity limitation in solving large-scale sparse linear equation systems through the coordinated operation of key steps such as multi-wavefront numerical decomposition, asynchronous disk storage, cyclic numerical updates, and back-substitution. The various steps are tightly linked through task dependencies and data flows, achieving optimal computational efficiency and resource utilization.

[0075] Specifically, in the sparse direct solution method based on data kernel external storage of the present invention, step 2 includes:

[0076] Obtain the dimension information and unprocessed principal component index number of the current wavefront matrix from the task dependencies;

[0077] The range of the principal index of the current region is calculated using the set region size parameters and the unprocessed principal index number.

[0078] The wavefront matrix is ​​divided into multiple data regions based on the range of the principal index.

[0079] Calculate the floating-point computation cost for each data segment and perform a balance adjustment.

[0080] In the sparse direct solution method based on off-core storage of the present invention, step 2 involves the partitioning and numerical decomposition calculation of the wavefront matrix. This step is based on the task dependency relationship established in step 1, and the fragment size parameters are set by the system memory capacity and disk performance parameters to achieve efficient data processing. Specifically, step 2 includes obtaining the dimension information and unprocessed principal component index number of the current wavefront matrix from the task dependency relationship. This information is directly derived from the task dependency relationship established by the non-zero element distribution pattern and matrix topology structure of the analytical coefficient matrix in step 1, ensuring the accuracy and consistency of data acquisition.

[0081] The primary index range of the current region is calculated using the set region size parameters and the unprocessed primary index number. The region size parameters are calculated based on system memory capacity and disk performance parameters, while the unprocessed primary index number is extracted from task dependencies. The primary index range is determined through mathematical operations, thus providing a precise calculation basis for subsequent partitioning. The primary index range calculation process dynamically adjusts using the region size parameters and the unprocessed primary index number to ensure the rationality and adaptability of the calculated range.

[0082] The wavefront matrix is ​​divided into multiple data regions based on the principal component index range. This division process uses the boundaries determined by the principal component index range to segment the wavefront matrix into multiple independent data blocks, each corresponding to a specific computational task. The partitioned data regions maintain structural integrity and computational independence, facilitating subsequent parallel processing and data management.

[0083] The floating-point computation load for each data segment is calculated and balanced. The floating-point computation load is calculated based on the matrix dimension and computational complexity of the data segment, and an algorithm is used to evaluate the computational load of each segment. The balancing process dynamically adjusts the segment size or computation order based on the floating-point computation results, bringing the computational load of each data segment closer to equilibrium, thereby improving overall computational efficiency and resource utilization. This balancing adjustment achieves a uniform distribution of computational tasks, avoids resource bottlenecks, and supports the stable solution of large-scale sparse linear equation systems.

[0084] Step 2 forms a continuous technical chain among its sub-steps, from data acquisition and calculation scope determination to area partitioning and load balancing. These tightly linked logical steps collectively support the smooth execution of multi-wavefront numerical decomposition. The entire process relies on task dependencies and system parameters to optimize the allocation of computing resources, providing high-quality decomposition factor data for subsequent steps.

[0085] Specifically, the sparse direct solution method based on off-core storage of the present invention calculates the primary index range of the current region from the unprocessed primary index number obtained from the task dependency relationship and the region size parameter based on the system memory capacity and disk performance parameter settings, including:

[0086] Obtain the unprocessed principal index number q and the total number of principal elements N from the task dependencies of the current wavefront matrix;

[0087] Read the area size parameter S calculated based on system memory capacity and disk performance parameters;

[0088] The terminating principal index p of the current region is determined by the formula p = q + min( floor(S), N - q ), where floor represents the floor function.

[0089] The primary index range [q, p] is used as the calculation range for the current region.

[0090] In the sparse direct solution method based on data kernel-external storage of this invention, the process of calculating the principal component index range of the current region achieves efficient data partitioning based on the structural information provided by the task dependency relationship and the system parameter settings. First, the unprocessed principal component index number and the total number of principal components of the current wavefront matrix are obtained from the task dependency relationship. These parameters are directly derived from the matrix topology analysis and non-zero element distribution analysis results established in step one, providing basic data support for the index range calculation. The unprocessed principal component index number represents the position of the principal component to be processed, and the total number of principal components reflects the overall dimension of the wavefront matrix.

[0091] The system reads the partition size parameter, calculated based on system memory capacity and disk performance parameters. This parameter is generated by detecting system hardware characteristics and optimizing configuration, reflecting the resource allocation strategy for memory and disk collaboration. As a dynamic adjustment factor, the partition size parameter directly affects the partition granularity and computational load distribution.

[0092] The terminating principal index of the current region is determined through mathematical operations. The calculation process combines the unprocessed principal index numbers, the total number of principal elements, and the region size parameter for comprehensive evaluation. The size parameter is processed using a floor function to ensure that the index value is an integer, and the minimum value is compared with the remaining number of principal elements to prevent index out-of-bounds errors. Determining the terminating principal index achieves precise partition boundary delineation.

[0093] The calculation range of the current slice is defined by the pivot index range, which defines the processing boundary of the data slice and guides subsequent numerical decomposition operations. The index range starts from the unprocessed pivot number and ends at the calculated termination index, forming a continuous processing interval. The explicit definition of the calculation range provides a standard basis for wavefront matrix partitioning and supports balanced adjustment of floating-point operation load.

[0094] The entire calculation process establishes a complete derivation chain for the principal component index range through the close integration of four stages: parameter acquisition, size reading, index determination, and range definition. Each stage relies on the structural information provided by task dependencies and the optimized results of system parameter settings to achieve partitioned calculation and resource condition matching. The calculation of the principal component index range lays the foundation for wavefront matrix partitioning and promotes the efficient execution of subsequent numerical decomposition.

[0095] Specifically, in the sparse direct solution method based on data kernel external storage of the present invention, step 3 includes:

[0096] Create a producer-consumer model, which includes a producer role computation thread and a consumer role input / output thread;

[0097] Initialize the pending input / output message queue and the completed message queue. The two queues adopt a first-in-first-out (FIFO) management strategy.

[0098] The producer role calculation thread adds the decomposition factor data to the tail of the input / output message queue;

[0099] The consumer role's input / output thread retrieves the decomposition factor data from the head of the pending input / output message queue;

[0100] The consumer role's input / output thread performs disk storage of the decomposition factor data;

[0101] The consumer role's input / output thread adds disk write completion status information to the completed message queue.

[0102] In the sparse direct solution method based on off-core storage of the present invention, step 3 involves the implementation of asynchronous disk storage operations. A producer-consumer model is used to efficiently manage the data writing process, thereby optimizing memory usage and disk I / O performance. Specifically, step 3 includes creating a producer-consumer model, where the producer role is undertaken by a computation thread, responsible for generating factorization data, and the consumer role is undertaken by an input / output thread, responsible for executing disk storage tasks. This role division achieves logical separation of computation and storage, supporting parallel processing.

[0103] Initialize the pending input / output message queue and the completed message queue. Both queues employ a first-in, first-out (FIFO) management strategy to ensure the sequentiality and fairness of data processing. The pending input / output message queue is used to temporarily store the factorization data to be written to disk, while the completed message queue records the disk storage completion status information. The queue initialization process allocates necessary memory resources and sets up synchronization mechanisms.

[0104] The producer thread adds the factorization data to the tail of the input / output message queue. This addition operation follows the queue's first-in, first-out (FIFO) principle, ensuring data enters the queue in the order it was generated, thus preventing data loss or out-of-order delivery. The producer thread maintains the integrity of the queue while adding data, ensuring reliable data retrieval for the consumer thread.

[0105] The consumer input / output thread retrieves factorization data from the head of the message queue. The retrieval operation follows a first-in, first-out (FIFO) strategy, ensuring that the earliest data entering the queue is processed first. The consumer thread actively monitors the queue status and immediately retrieves data when the queue is not empty, reducing waiting time.

[0106] The consumer-side input / output thread performs disk storage of the factorization data. This storage operation writes the data to disk for persistence, leveraging the large disk capacity to alleviate memory pressure. The storage procedure is asynchronous, allowing computation threads to continue execution without blocking, thus improving overall throughput.

[0107] The consumer's input / output thread adds disk write completion status information to the completed message queue. This status information includes a storage success flag and a data location index, facilitating subsequent querying and verification. Updates to the completed message queue notify other components of the storage progress, supporting collaborative work.

[0108] Step 3, through the collaboration of the producer-consumer model, achieves pipelined processing of computation and storage, reducing I / O wait time and improving the efficiency of large-scale data processing. Each sub-step is interconnected, forming a closed loop from data generation to storage completion, thus achieving rational utilization of system resources and reliable data persistence.

[0109] Specifically, in the sparse direct solution method based on data kernel external storage of the present invention, step 3 further includes:

[0110] Atomic variables are created and used to implement synchronous control for adding and deleting messages in the pending input / output message queue and the completed message queue;

[0111] Create and initialize global state variables, and use global state variables to pass input and output operation status signals between the producer role calculation thread and the consumer role input and output thread;

[0112] When the producer's computation thread and the consumer's input / output thread access the message queue, a spin-wait mechanism is used to ensure thread-safe queue access.

[0113] In the sparse direct solution method based on off-core data storage of this invention, the additional implementation of step 3 involves synchronization control and state management mechanisms in a multi-threaded environment. These mechanisms provide reliable concurrent execution guarantees for the producer-consumer model. Atomic variables are created, and synchronization control is implemented for adding and deleting messages in the pending input / output message queues and completed message queues. As underlying synchronization primitives, atomic variables ensure the atomicity and visibility of concurrent access operations to the message queues, preventing data inconsistency caused by multiple threads simultaneously modifying the queue structure. Synchronization control covers the two key operations of adding and deleting messages in the queues, ensuring the order and integrity of data addition by the producer thread and data retrieval by the consumer thread.

[0114] A global state variable is created and initialized to transmit input / output operation status signals between the producer's computation thread and the consumer's input / output thread. This global state variable acts as a communication medium between threads, storing current input / output operation status information, including disk storage progress, queue load, and error status codes. The status signal transmission employs a memory visibility guarantee mechanism, enabling the producer thread to promptly perceive the execution status of the consumer thread, and the consumer thread to accurately obtain the data generation progress of the producer thread, supporting bidirectional state synchronization.

[0115] When producer threads (computation threads) and consumer threads (input / output threads) access the message queue, a spin-wait mechanism is employed to ensure thread-safe queue access. This mechanism achieves thread blocking by continuously checking resource availability. When the message queue is inaccessible, the thread remains actively waiting instead of being immediately suspended, reducing context switching overhead. Thread-safe queue access guarantees that multiple producer and consumer threads can concurrently operate on the message queue while maintaining data consistency and the correctness of the operation sequence.

[0116] These synchronization and communication mechanisms together form the basic architecture for multi-threaded collaboration. Atomic variables provide low-level synchronization guarantees, global state variables achieve high-level state coordination, and the spin-wait mechanism optimizes resource access efficiency. These mechanisms work together to form a complete thread-safe solution, supporting the stable operation of the producer-consumer model in high-concurrency environments and providing reliable technical support for asynchronous disk storage operations.

[0117] Specifically, in the sparse direct solution method based on data kernel external storage of the present invention, step 4 includes:

[0118] Based on the task dependencies established in step 1, determine the subsequent wavefront matrix to be updated;

[0119] Read the decomposition factor data stored in step 3 from the disk;

[0120] Create a numerical update memory buffer and load the read decomposition factor data into the numerical update memory buffer;

[0121] A multi-wavefront numerical update algorithm is adopted, which uses the loaded decomposition factor data to perform numerical update calculations on the subsequent wavefront matrix.

[0122] In the sparse direct solution method based on off-core data storage of the present invention, step 4 involves the numerical update process of subsequent wavefront matrices. This process performs iterative calculations based on task dependencies and decomposition factor data stored on disk to achieve an asymptotic solution of large-scale sparse linear equation systems. Specifically, step 4 includes determining the subsequent wavefront matrices to be updated based on the task dependencies established in step 1. The task dependencies are derived from the non-zero element distribution pattern of the coefficient matrix and matrix topology analysis, defining the computational order and data dependencies between wavefront matrices, and guiding the selection and scheduling of the update order.

[0123] The decomposition factor data stored in step 3 is read from the disk. In step 3, the decomposition factor data is written to the disk via asynchronous disk storage operations. Here, the read operation retrieves the required data based on the storage location and data type, supporting subsequent calculations. The read process utilizes disk I / O channels and file system interfaces to efficiently obtain persistent data.

[0124] A numerical update memory buffer is created and the read decomposition factor data is loaded into it. The numerical update memory buffer provides dedicated memory space for update calculations. The load operation copies data from disk into memory, reducing I / O overhead during the calculation process. Buffer management dynamically allocates and releases based on memory capacity and data size.

[0125] A multi-wavefront numerical update algorithm is employed, using the loaded decomposition factor data to perform numerical update calculations on subsequent wavefront matrices. The algorithm utilizes the decomposition factor data to perform matrix operations and numerical calculations, updating the element values ​​of the wavefront matrices. The algorithm executes matrix multiplication and addition operations to generate updated matrix data. The update calculation is repeated cyclically until all wavefront matrices have been processed.

[0126] Step 4 guides the update order through task dependencies, provides input support through disk data reading, optimizes data access through memory buffers, and executes core calculations using the multi-wavefront algorithm, forming an efficient and accurate numerical update process. The sub-steps are logically interconnected, supporting cyclic processing of large-scale data, ultimately completing the decomposition of all wavefront matrices.

[0127] Specifically, in the sparse direct solution method based on data kernel external storage of the present invention, step 2 further includes:

[0128] The overall memory requirements for numerical decomposition are estimated based on the dimension of the coefficient matrix and the non-zero element distribution pattern.

[0129] Based on the estimated total memory requirements and the detected available physical memory capacity, the adjustment value of the region size parameter is calculated through an optimization algorithm;

[0130] Based on the detected disk input / output performance characteristics, an optimized value for the asynchronous storage memory buffer size is calculated using a performance optimization algorithm.

[0131] In the sparse direct solution method based on data kernel-external storage of the present invention, step 4 involves the numerical update process of subsequent wavefront matrices. This process is based on the task dependency relationship established in step 1, and realizes the asymptotic solution of large-scale sparse linear equation systems through a loop execution mechanism. Specifically, step 4 includes determining the subsequent wavefront matrices to be updated according to the task dependency relationship established in step 1. The task dependency relationship is derived from the non-zero element distribution pattern of the coefficient matrix and the matrix topology analysis, which clearly defines the calculation order and data dependency relationship between each wavefront matrix, providing logical guidance for the update operation.

[0132] The decomposition factor data stored in step 3 is read from the disk. In step 3, the decomposition factor data is persistently saved to the disk medium through asynchronous disk storage operations. The read operation retrieves the required data blocks based on the storage location information and data type identifier, providing input data support for numerical updates. The read process utilizes disk I / O channels and file system interfaces to achieve efficient data retrieval and transmission.

[0133] A numerical update memory buffer is created and the read decomposition factor data is loaded into it. This buffer allocates dedicated memory space for update calculations. The loading operation copies data blocks from disk to the memory working area, reducing direct I / O access overhead during computation. Buffer management dynamically allocates and releases memory based on system memory capacity and data block size, optimizing memory resource utilization.

[0134] A multi-wavefront numerical update algorithm is employed, using loaded decomposition factor data to perform numerical update calculations on subsequent wavefront matrices. Based on numerical analysis principles, the algorithm utilizes decomposition factor data for matrix operations and numerical computation to update the values ​​of each element in the wavefront matrix. The algorithm execution process includes matrix multiplication, matrix addition, and numerical correction operations, generating updated matrix data blocks. The update calculation adopts a cyclic execution mode until all wavefront matrices have been processed, and the execution status is tracked through a monitoring mechanism.

[0135] Step 4 guides the update order through task dependencies, provides input support through disk data reading, optimizes data access through memory buffers, and executes core calculations using the multi-wavefront algorithm, forming a complete and efficient numerical update process. The sub-steps are tightly integrated, supporting iterative processing of large-scale data and ultimately completing the decomposition of all wavefront matrices, laying the foundation for subsequent back-substitution solutions.

[0136] Specifically, in the sparse direct solution method based on data kernel external storage of the present invention, step 2 further includes:

[0137] Before performing the left-look numerical decomposition calculation in step 2, the mathematical properties of the coefficient matrix are detected by a matrix analysis algorithm;

[0138] When a symmetric positive definite matrix is ​​detected, the Cholesky decomposition algorithm is used to perform numerical decomposition.

[0139] When an asymmetric matrix is ​​detected, the LU decomposition algorithm is used to perform numerical decomposition.

[0140] During the numerical decomposition calculation process, both real and complex number data types are processed simultaneously.

[0141] In the sparse direct solution method based on off-core storage of the present invention, the additional implementation of step 2 involves matrix mathematical property detection and adaptive decomposition algorithm selection. This process is performed before the left-look numerical decomposition calculation, and improves computational efficiency and numerical stability through intelligent algorithm selection. Specifically, step 2 includes detecting the mathematical properties of the coefficient matrix through a matrix analysis algorithm before performing the left-look numerical decomposition calculation in step 2. The matrix analysis algorithm is based on linear algebra and matrix theory, analyzes the symmetry, positive definiteness, and special structural features of the coefficient matrix, and generates a mathematical property detection report, providing a decision basis for subsequent algorithm selection.

[0142] When a symmetric positive definite matrix is ​​detected, the Cholesky decomposition algorithm is used to perform numerical decomposition. The Cholesky decomposition algorithm is designed for the special mathematical properties of symmetric positive definite matrices. It uses matrix decomposition techniques to decompose the coefficient matrix into the product of a lower triangular matrix and its transpose. This algorithm has the advantages of low computational cost and numerical stability, and is suitable for handling mathematical models with symmetric structures.

[0143] When an asymmetric matrix is ​​detected, the LU decomposition algorithm is used to perform numerical decomposition. The LU decomposition algorithm is applicable to general matrix structures, decomposing the coefficient matrix into the product of lower triangular and upper triangular matrices. It supports partial pivoting to improve numerical stability. This algorithm has the characteristics of strong versatility and wide applicability.

[0144] During the numerical decomposition calculation, both real and complex number data types are processed simultaneously. The calculation process supports mathematical operations in both the real and complex number fields, including addition, subtraction, multiplication, and division of real numbers and addition, subtraction, multiplication, and division of complex numbers, as well as conjugate operations. Data type processing maintains calculation accuracy and algorithm consistency, meeting the numerical calculation needs of different application scenarios.

[0145] The entire technical solution intelligently selects the most suitable decomposition algorithm through pre-matrix characteristic detection. Cholesky decomposition optimizes the computational efficiency of symmetric positive definite matrices, while LU decomposition ensures the numerical stability of asymmetric matrices. It also accommodates both real and complex data types, forming a complete adaptive numerical decomposition system. Each step is logically interconnected: matrix analysis provides the basis for algorithm selection, algorithm execution provides the results for numerical decomposition, and data type processing ensures computational accuracy, collectively improving the performance and reliability of the sparse direct solution method.

[0146] Specifically, in the sparse direct solution method based on data kernel external storage of the present invention, step 4 of the loop execution state monitoring mechanism through input and output includes:

[0147] The status monitoring mechanism collects disk storage operation status signals through input and output;

[0148] Based on the monitored input and output operation status signals, a dynamic scheduling algorithm is used to dynamically schedule the execution tasks of the producer role computing threads.

[0149] The scheduling of computational tasks and disk input / output tasks is carried out in parallel, completing the overlap of computation and input / output operations.

[0150] In the sparse direct solution method based on off-core storage of the present invention, the input / output completion status monitoring mechanism in step 4 involves the coordinated work of disk storage operation status monitoring and dynamic scheduling of computing tasks. This mechanism achieves overlapping execution of computing and input / output operations through status signal acquisition and task scheduling. Specifically, step 4 includes acquiring disk storage operation status signals through the input / output completion status monitoring mechanism. The input / output completion status monitoring mechanism continuously monitors the execution progress and completion status of disk storage operations. The acquired operation status signals include information such as disk write progress, data block integrity flags, and error status codes. The status signals reflect the current actual execution status of disk storage.

[0151] Based on the monitored input and output operation status signals, a dynamic scheduling algorithm is used to dynamically schedule the execution tasks of the producer role computing threads. The algorithm analyzes system load and resource utilization based on real-time collected operation status signals, intelligently adjusting the task allocation and execution order of the producer role computing threads. The scheduling process considers disk storage progress and computing resource availability, prioritizing the scheduling of computing tasks associated with completed disk storage operations to reduce data waiting time.

[0152] The scheduler coordinates the execution timing of computation tasks and disk I / O tasks, allowing computation and I / O operations to overlap. The scheduler synchronizes the execution of computation threads and I / O threads, ensuring that computation tasks execute synchronously during disk I / O operations. Parallel execution leverages the concurrent processing capabilities of multi-core processors; computation tasks process ready data while I / O tasks perform new disk storage operations, achieving time overlap between computation and I / O operations.

[0153] The entire monitoring and scheduling process forms a closed-loop control system. Status monitoring provides real-time system feedback, dynamic scheduling makes decisions and adjustments based on feedback, and parallel execution maximizes resource utilization. Each link is interconnected through status signal transmission and algorithm decision-making. The input-output status monitoring mechanism provides data support for dynamic scheduling, and the dynamic scheduling algorithm provides task allocation schemes for parallel execution. Ultimately, this achieves an efficient overlap of computation and input-output operations, significantly improving the overall performance of large-scale sparse direct solutions.

[0154] Specifically, in the sparse direct solution method based on off-core storage of the present invention, step 5, which involves performing numerical calculations on the decomposition factor data in the back-substitution solution memory working area using a previous-back substitution algorithm, includes:

[0155] Based on the solution order determined by the task dependencies, read the complete stored decomposition factor data from the disk in reverse order;

[0156] Load the read decomposition factor data into the back-substitution solution memory work area;

[0157] According to the calculation order required by the previous generation and back generation algorithm, the previous generation calculation and back generation calculation steps of the previous generation and back generation algorithm are executed in sequence.

[0158] By merging the results of previous and back-substitution calculations, a numerical solution to a high-dimensional sparse linear system of equations is generated.

[0159] In the sparse direct solution method based on off-core storage of the present invention, step 5 involves the process of performing numerical calculations on the decomposition factor data in the back-substitution solution memory working area using a previous-back-substitution algorithm. This process is based on the solution order determined by task dependencies, reading data from the disk and performing calculations to generate the final numerical solution. Specifically, step 5 includes reading the completely stored decomposition factor data from the disk in reverse order according to the solution order determined by task dependencies. The task dependencies originate from the multi-wavefront numerical decomposition task dependencies established in step 1, defining the calculation order and data dependencies between wavefront matrices. Reverse order reading enables back-substitution solution to process data backward from the last wavefront matrix, satisfying the data dependency requirements. The read operation retrieves the decomposition factor data stored in steps 3 and 4 from the disk, which includes the decomposition results of all wavefront matrices.

[0160] The decomposition factor data is loaded into the back-substitution solution memory work area, which is a memory area specifically established in step 5 to store and process the data required for back-substitution. The loading operation transfers data blocks from the disk to the memory work area for efficient numerical computation. The memory work area is managed based on the system memory capacity and data size, and memory resources are dynamically allocated to avoid memory overflow.

[0161] Following the computational order required by the previous-backward substitution algorithm, the previous-generation and backward substitution steps of the algorithm are executed sequentially. The previous-backward substitution algorithm is a numerical computation method. The previous-generation calculation involves solving the lower triangular system from front to back, while the backward substitution calculation involves solving the upper triangular system from back to front. The computational order is determined by the algorithm requirements to achieve numerical stability and accuracy. The execution process includes matrix-vector multiplication and linear algebra operations, and supports the processing of real and complex number data types.

[0162] The results of previous and back-substitution calculations are combined to generate a numerical solution to a high-dimensional sparse linear system of equations. The merging operation combines the results of previous and back-substitution calculations into a complete solution vector, which represents the final numerical solution of the high-dimensional sparse linear system of equations. This solution vector can be used for subsequent analysis or output. The result generation process is based on the principle of numerical computation, ensuring the correctness and accuracy of the solution.

[0163] The entire process of step 5, through the coherent flow of data reading, loading, calculation, and result merging, achieves the complete execution of back-substitution solution. The sub-steps are logically connected, the task dependencies guide the data reading order, the memory workspace optimizes the calculation efficiency, the previous-generation and back-substitution algorithms ensure numerical accuracy, and finally generate a reliable numerical solution.

[0164] This invention addresses the technical problem of soaring memory requirements caused by non-zero element padding in sparse matrix numerical decomposition by employing a sparse direct solution method based on off-core data storage. This method stores large-scale data on disk and loads only the data currently necessary for computation into memory, effectively alleviating physical memory capacity limitations.

[0165] This invention first obtains the coefficient matrix of a high-dimensional sparse linear system of equations, analyzes the distribution pattern of non-zero elements and the matrix topology, and establishes the task dependencies of the multi-wavefront numerical decomposition method. This step determines the computational order and data dependencies, laying the foundation for subsequent partitioning processing.

[0166] The system detects system memory capacity and disk performance parameters, and sets partition size parameters accordingly, dividing each wavefront matrix into multiple data partitions with balanced floating-point computational load. The partitioning strategy is optimized based on memory and disk performance, ensuring each partition size adapts to available memory and avoids memory overflow. Left-look numerical decomposition calculations are performed on each data partition to generate decomposition factor data, with memory usage controlled throughout the decomposition process.

[0167] A message queue for pending input / output and an input / output completion status monitoring mechanism are created to achieve asynchronous disk storage. After factorization data is stored in the message queue, it is asynchronously written to disk, and the memory buffer is released after storage is complete. This mechanism allows computation and storage to overlap, reducing memory usage time.

[0168] The system reads the required decomposition factor data from disk and performs subsequent wavefront matrix updates based on task dependencies. It iteratively executes the partitioning and update steps, tracking progress through a status monitoring mechanism until all wavefront matrix decompositions are completed. This cyclical processing allows large-scale problems to be completed in stages.

[0169] A back-substitution solution memory workspace is established. Complete factorization data is read from disk, and numerical computation is performed using a previous-back substitution algorithm to generate the final numerical solution. The back-substitution solution utilizes disk data, loading it into memory only when necessary, thus minimizing memory requirements.

[0170] By using partitioned computation, asynchronous storage, and cyclic updates, this invention shifts memory-intensive operations to disk, enabling the solution of large-scale sparse linear equations in a limited memory environment and supporting large-scale finite element simulation applications.

[0171] Specific embodiments of the present invention: This embodiment demonstrates the application of a sparse direct solution method based on off-core storage in practical industrial CAE simulation. This embodiment achieves efficient solution of a million-degree-of-freedom sparse linear equation system generated in automotive chassis structural strength simulation on a computer platform with 64GB of memory.

[0172] Hardware and environment configuration;

[0173] The hardware environment implemented in this embodiment is as follows: operating system: Windows 10, memory capacity: 64GB DDR4, processor: Intel Core i5-10400F CPU @ 2.90GHz (6 cores), storage device: 1TB NVMe SSD solid-state drive, numerical calculation library: self-developed solver based on this method.

[0174] This invention first obtains the coefficient matrix of a high-dimensional sparse linear equation system generated from the simulation of the structural strength of an automobile chassis. The coefficient matrix has a dimension of 1,200,000 × 1,200,000 and 18,650,000 non-zero elements.

[0175] By analyzing the non-zero element distribution pattern and matrix topology of the coefficient matrix using matrix analysis algorithms, it was found that the matrix is ​​a symmetric positive definite matrix. The task dependencies of the multi-wavefront numerical decomposition method were established, determining that a total of 358 wavefront matrices need to be processed, with the largest wavefront matrix having a dimension of 85,000 × 85,000.

[0176] The system's memory capacity (64GB) and disk performance parameters (SSD sequential write speed 3.2GB / s, read speed 3.5GB / s) were detected. Based on the memory capacity and disk performance parameters, the partition size parameter S=16,384 was calculated and set.

[0177] Based on task dependencies and slice size parameters, each wavefront matrix is ​​divided into multiple data slices with balanced floating-point operations. For the largest wavefront matrix with dimensions of 85,000×85,000, the slice division is calculated according to the formula p=q+min(floor(S), Nq), resulting in a total of 6 data slices. The floating-point operations of each slice differ by no more than 5% after equalization adjustment.

[0178] Left-look numerical decomposition is performed on each data patch using the Cholesky decomposition algorithm (which detects symmetric positive definite matrices) to generate decomposition factor data for the data patch. The decomposition process also handles real number data types.

[0179] Create a producer-consumer model, including 6 producer role computation threads and 1 consumer role input / output thread. Initialize the pending input / output message queue and the completed message queue, both using a first-in, first-out (FIFO) management strategy.

[0180] The producer's computation thread adds the factorization data to the tail of the input / output message queue, while the consumer's input / output thread retrieves data from the head of the queue and performs asynchronous disk storage. Atomic variables are created and synchronization control is implemented for queue operations, using a spin-wait mechanism to ensure thread-safe queue access.

[0181] The disk storage completion status is monitored through an input / output completion status monitoring mechanism, and the corresponding memory buffer is released after storage is complete. In this embodiment, the asynchronous storage memory buffer size is optimized to 2GB, effectively overlapping computation and I / O operations.

[0182] Determine the subsequent wavefront matrix to be updated based on task dependencies, and read the stored decomposition factor data from disk. Create a numerical update memory buffer (8GB in size) and load the decomposition factor data.

[0183] A multi-wavefront numerical update algorithm is employed, using the loaded decomposition factor data to perform numerical update calculations on subsequent wavefront matrices. The algorithm iteratively executes the steps from patch partitioning to numerical update, and a state monitoring mechanism monitors the loop execution status through input / output completion, dynamically scheduling computational thread tasks.

[0184] This embodiment completes 357 value update cycles, ultimately completing all wavefront matrix decompositions.

[0185] Establish a back-substitution solution memory workspace (12GB in size), and read the complete stored decomposition factor data from the disk in reverse order according to the solution order determined by the task dependencies.

[0186] The read decomposition factor data is loaded into the back-substitution solution memory workspace, and the previous generation calculation and back-substitution calculation steps are executed sequentially according to the calculation order required by the previous generation and back-substitution algorithm. The results of the previous generation calculation and back-substitution calculation are merged to generate a numerical solution to the high-dimensional sparse linear equation system.

[0187] Table 1. Detailed test cases and results information

[0188]

[0189] Table 1 provides detailed test results and use case information. It can be seen that, compared with MKL Pardiso, this invention successfully reduced peak memory by an average of about 7%, and by a maximum of 30%; and reduced computation time by an average of about 38%, and by a maximum of 56%.

[0190] The technical features of this invention are explained below:

[0191] The multi-wavefront method is a numerical approach for directly solving sparse matrices. Its construction process is based on the analysis of the non-zero element distribution pattern and matrix topology of the coefficient matrix, establishing dependencies between tasks through elimination trees or dependency graphs in graph theory. The multi-wavefront method decomposes a large-scale sparse matrix into multiple smaller, denser sub-matrices called wavefront matrices, which are processed in an order determined by the dependencies. In this invention, the multi-wavefront method guides the entire numerical decomposition process. Step 1 establishes task dependencies by analyzing the coefficient matrix; steps 2 and 4 perform partitioning and numerical updates based on these dependencies. The data processing involves the generation, partitioning, and factor calculation of the wavefront matrix, with the output being structured decomposed factor data and the updated wavefront matrix.

[0192] The producer-consumer model is a concurrent programming model. Its construction involves creating two types of thread roles: producer threads are responsible for generating data, and consumer threads are responsible for processing data. They communicate through a shared buffer or message queue. In this invention, this model is applied to the asynchronous disk storage process in step 3. The producer thread, responsible for computation, places the generated decomposition factor data into a pending input / output message queue, while the consumer thread, responsible for input / output, retrieves data from the queue and performs a disk write operation. The data processing process includes data queue management, synchronization control, and status monitoring. The output is a persistently stored decomposition factor file and a status signal indicating the completion of the operation.

[0193] Cholesky decomposition is a matrix factorization algorithm for symmetric positive definite matrices. Its construction process is based on the symmetric positive definite property of a matrix, decomposing it into the product of a lower triangular matrix and its transpose. In this invention, Cholesky decomposition is used in the left-look numerical decomposition calculation in step 2. This algorithm is invoked when the coefficient matrix is ​​detected to be symmetric positive definite through matrix analysis. The data processing includes numerical decomposition calculations for data regions, and the output is the lower triangular factor matrix of that region.

[0194] LU decomposition is a matrix factorization algorithm applicable to general square matrices. The construction process decomposes the matrix into the product of a lower triangular matrix and an upper triangular matrix, typically requiring a pivoting operation to ensure numerical stability. In this invention, LU decomposition is used in the left-look numerical decomposition calculation in step 2. This algorithm is invoked when the matrix analysis algorithm detects that the coefficient matrix is ​​asymmetric. The data processing includes decomposition calculation of the data region and pivot selection, with the output being the lower triangular factor matrix and the upper triangular factor matrix of that region.

[0195] The pre- and back-substitution algorithm is a numerical algorithm for solving trigonometric linear equation systems. The construction process includes two stages: pre-substitution calculation and back-substitution calculation. Pre-substitution calculation solves the lower trigonometric system, and back-substitution calculation solves the upper trigonometric system. In this invention, the pre- and back-substitution algorithm is applied to the back-substitution solution stage in step 5 to generate the final numerical solution of the high-dimensional sparse linear equation system. The data processing includes reading the decomposition factor data from disk, loading it into the memory workspace, and executing the pre- and back-substitution calculation steps. The output is the numerical solution vector of the original equation system.

[0196] The coefficient matrix of a high-dimensional sparse linear equation system is a matrix containing a large number of zero elements and with extremely high dimensionality, typically originating from scientific and engineering computational fields such as finite element analysis and computational fluid dynamics. The construction process generates the matrix by discretizing partial differential equations or physical models, and the distribution of non-zero elements exhibits a specific pattern. In this invention, the coefficient matrix serves as input data. Step 1 establishes task dependencies by analyzing the distribution pattern of non-zero elements and the matrix topology, providing a basic data structure for subsequent numerical decomposition. The data processing includes matrix reading, pattern analysis, and structure resolution, with the output being structured matrix information and non-zero element distribution characteristics.

[0197] The multi-wavefront method is a numerical approach for directly solving sparse matrices. Its construction process is based on graph theory analysis of the coefficient matrix, establishing computational order relationships between tasks by eliminating trees or dependency graphs. In this invention, the multi-wavefront method guides the entire solution process. Step 1 establishes task dependencies by analyzing the matrix topology, and step 4 performs numerical updates based on these dependencies. The data processing involves the generation, sorting, and computation scheduling of the wavefront matrix, and the output is a structured task dependency graph.

[0198] The segment size parameter is a partitioning parameter calculated based on system memory capacity and disk performance. The construction process involves detecting hardware performance characteristics and using optimization algorithms to calculate this parameter. In this invention, the segment size parameter is used for data segmentation in step 2. The calculation formula p = q + min(floor(S), Nq) determines the partition boundaries based on the unprocessed principal index and the total number of principal elements. The data processing includes parameter calculation and dynamic adjustment, and the output is an optimized segment partitioning scheme.

[0199] Left-look numerical decomposition is a matrix factorization calculation method. Its construction process is based on matrix block decomposition technology, performing local decomposition operations on data regions. In the technical solution of this invention, left-look numerical decomposition is applied to the data region processing in step 2 to generate decomposition factor data. The data processing process includes matrix block calculation and local decomposition operations, and the output result is the decomposition factor matrix of the data region.

[0200] The producer-consumer model is a concurrent programming paradigm. Its construction involves creating threads that act as both producers and consumers, exchanging data through a shared buffer. In this invention, this model is used for asynchronous disk storage in step 3. The computation thread acts as the producer, generating data, while the I / O thread acts as the consumer, processing the storage. The data processing includes thread creation, queue management, and task coordination, resulting in a highly efficient data pipeline processing architecture.

[0201] Asynchronous disk storage is a non-synchronous data persistence technology. Its construction process separates computation and storage operations, enabling them to execute in parallel. In this invention, asynchronous disk storage is implemented in step 3, where factorization factor data is written to disk using a message queue mechanism. The data processing includes data serialization, disk writing, and status monitoring, with the output being a persistently stored factorization factor file.

[0202] The input / output completion status monitoring mechanism is a performance monitoring system that collects and analyzes disk operation status signals during its construction. In this invention, this mechanism is used in steps 3 and 4 to monitor the completion of storage operations and release memory buffers. The data processing includes status acquisition, signal analysis, and resource management, and the output is real-time I / O operation status information.

[0203] The numerical update loop is an iterative computation process built upon task dependencies, repeatedly executing data reading and matrix update operations. In this invention, the numerical update loop is implemented in step 4, reading data from the disk and updating the subsequent wavefront matrix. The data processing includes loop control, data loading, and matrix operations, with the output being the updated wavefront matrix data.

[0204] The back-substitution solution memory workspace is a specially allocated memory space, dynamically allocated and managed according to the back-substitution solution requirements during its construction. In this invention, the back-substitution solution memory workspace is created in step 5 and is used to store the decomposition factor data read from the disk. The data processing process includes memory allocation, data loading, and space management, and the output is an optimized memory working environment.

[0205] The pre- and back-substitution algorithm is a numerical method for solving trigonometric equation systems. The construction process includes two stages: pre-substitution calculation and back-substitution calculation, handling lower trigonometric and upper trigonometric systems respectively. In this invention, the pre- and back-substitution algorithm is applied to the back-substitution solution stage in step 5 to generate the final numerical solution. The data processing includes matrix-vector multiplication and linear solving operations, and the output is the numerical solution vector of the original equation system.

[0206] Cholesky decomposition is a decomposition algorithm for symmetric positive definite matrices. The construction process decomposes the matrix into the product of a lower triangular matrix and its transpose. In this invention, Cholesky decomposition is used for the numerical decomposition in step 2, and this algorithm is employed when a symmetric positive definite matrix is ​​detected. The data processing includes matrix decomposition and factorization, and the output is a lower triangular factor matrix.

[0207] LU decomposition is a general matrix factorization method that constructs a matrix by decomposing it into the product of lower triangular and upper triangular matrices. In this invention, LU decomposition is used for the numerical decomposition in step 2, and this algorithm is employed when an asymmetric matrix is ​​detected. The data processing includes partial pivoting and matrix decomposition, and the output is lower triangular and upper triangular factor matrices.

[0208] Dynamic scheduling algorithms are resource allocation strategies that adjust the task execution order in real time based on the system state during the construction process. In the technical solution of this invention, the dynamic scheduling algorithm is used in step 4 to adjust the execution of computing tasks according to the I / O operation state. The data processing process includes state analysis and task scheduling, and the output result is an optimized task execution sequence.

[0209] Spin-and-wait is a thread synchronization technique that coordinates threads by repeatedly checking resource availability during its construction. In this invention, the spin-and-wait mechanism is used for message queue access in step 3 to ensure thread-safe queue operations. The data processing process includes resource detection and thread synchronization, and the output is a safe multi-threaded access environment.

[0210] Atomic variables are synchronization primitives in concurrent programming, providing uninterrupted read and write operations to ensure data consistency during their construction. In this invention, atomic variables are used for queue operation synchronization in step 3, enabling atomic access to shared resources. The data processing process includes atomic operations and memory synchronization, and the output is shared data with guaranteed consistency.

[0211] Global state variables serve as a centralized storage for system state information, and their construction process provides a mechanism for state sharing and communication between threads. In this invention, global state variables are used for producer-consumer communication in step 3 to transmit I / O operation status signals. The data processing process includes state updates and signal transmission, and the output is a unified view of the system state.

[0212] Floating-point load balancing is a load optimization technique that involves calculating and adjusting the computational load of each computing unit. In this invention, floating-point load balancing is used for data partitioning in step 2 to achieve a basic balance in computational load across partitions. The data processing process includes computational load calculation and partition adjustment, and the output is a load-balanced distribution of computational tasks.

[0213] Message queues are an inter-process communication mechanism that uses a first-in, first-out (FIFO) strategy to manage data transmission during their construction. In this invention, a message queue is used for asynchronous storage in step 3, including a queue of pending input / output and a queue of completed processes. The data processing includes queue management and message passing, and the output is an ordered data transmission channel.

[0214] Matrix analysis algorithms are tools for detecting mathematical properties, and their construction process is based on linear algebra theory to analyze the mathematical properties of matrices. In the technical solution of this invention, the matrix analysis algorithm is used for preprocessing in step 2 to detect the symmetry and positive definiteness of the matrix. The data processing process includes property detection and classification judgment, and the output is a report of the mathematical properties of the matrix.

[0215] Task dependencies are constraints describing the computational order, and the construction process determines the dependencies between computational steps by analyzing the matrix structure. In this invention, task dependencies are established in step 1, guiding the computational order of the entire solution process. The data processing includes dependency analysis and relationship establishment, and the output is a structured task scheduling graph.

[0216] The memory buffer is a temporary data storage area, and its construction process dynamically allocates memory space according to computational needs. In the technical solution of this invention, the memory buffer is used for numerical updates in step 4 and back-substitution in step 5, storing data read from the disk. The data processing process includes buffer allocation and data caching, resulting in an efficient data access environment.

[0217] Performance optimization algorithms are system tuning methods that calculate optimal parameter configurations by analyzing hardware characteristics. In this invention, the performance optimization algorithm is used for parameter setting in step 2, calculating optimized values ​​for the slice size and buffer size. The data processing includes performance analysis and parameter optimization, and the output is optimized system configuration parameters.

Claims

1. A sparse direct solution method based on data kernel external storage, characterized in that, include: Step 1: Obtain the coefficient matrix of the high-dimensional sparse linear equation system to be solved, analyze the non-zero element distribution pattern and matrix topology of the coefficient matrix, establish the task dependency relationship of the multi-wavefront numerical decomposition, and generate the wavefront matrix corresponding to the task dependency relationship. Step 2: Detect system memory capacity and disk performance parameters, calculate and set the slice size parameters based on memory capacity and disk performance parameters, divide each wavefront matrix into multiple data slices with balanced floating-point operation according to task dependencies and slice size parameters, perform left-look numerical decomposition calculation on each data slice, and generate decomposition factor data for the data slice. Step 3: Create a message queue to be input / output and an input / output completion status monitoring mechanism. Store the decomposition factor data into the message queue to be input / output, extract the decomposition factor data from the message queue to be input / output and perform asynchronous disk storage. Release the corresponding memory buffer after the disk storage is completed by the input / output completion status monitoring mechanism. Step 4: Read the required decomposition factor data from the disk, perform subsequent wavefront matrix numerical updates according to task dependencies, and repeat the process from patch division to numerical update. Monitor the loop execution status through an input / output completion status monitoring mechanism until all wavefront matrix decompositions are completed. Step 5: Establish a back-substitution solution memory workspace. Read the complete stored decomposition factor data from the disk and load it into the back-substitution solution memory workspace. Perform numerical calculations on the decomposition factor data using the previous-back substitution algorithm in the back-substitution solution memory workspace. Combine the calculation results of each step to generate a numerical solution for a high-dimensional sparse linear system of equations.

2. The sparse direct solution method based on data kernel external storage according to claim 1, characterized in that, Step 2 includes: Obtain the dimension information and unprocessed principal component index number of the current wavefront matrix from the task dependencies; The range of the principal index of the current region is calculated using the set region size parameters and the unprocessed principal index number. The wavefront matrix is ​​divided into multiple data regions based on the range of the principal index. Calculate the floating-point computation cost for each data segment and perform a balance adjustment.

3. The sparse direct solution method based on data kernel external storage according to claim 2, characterized in that, The primary index range of the current region is calculated using the unprocessed primary index number obtained from the task dependency relationship and the region size parameter based on the system memory capacity and disk performance parameters. Obtain the unprocessed principal index number q and the total number of principal elements N from the task dependencies of the current wavefront matrix; Read the area size parameter S calculated based on system memory capacity and disk performance parameters; The terminating principal index p of the current region is determined by the formula p = q + min( floor(S), N - q ), where floor represents the floor function. The primary index range [q, p] is used as the calculation range for the current region.

4. The sparse direct solution method based on data kernel external storage according to claim 3, characterized in that, Step 3 includes: Create a producer-consumer model, which includes a producer role computation thread and a consumer role input / output thread; Initialize the pending input / output message queue and the completed message queue. The two queues adopt a first-in-first-out (FIFO) management strategy. The producer role calculation thread adds the decomposition factor data to the tail of the input / output message queue; The consumer role's input / output thread retrieves the decomposition factor data from the head of the pending input / output message queue; The consumer role's input / output thread performs disk storage of the decomposition factor data; The consumer role's input / output thread adds disk write completion status information to the completed message queue.

5. The sparse direct solution method based on data kernel external storage according to claim 4, characterized in that, Step 3 also includes: Atomic variables are created and used to implement synchronous control for adding and deleting messages in the pending input / output message queue and the completed message queue; Create and initialize global state variables, and use global state variables to pass input and output operation status signals between the producer role calculation thread and the consumer role input and output thread; When the producer's computation thread and the consumer's input / output thread access the message queue, a spin-wait mechanism is used to ensure thread-safe queue access.

6. The sparse direct solution method based on data kernel external storage according to claim 5, characterized in that, Step 4 includes: Based on the task dependencies established in step 1, determine the subsequent wavefront matrix to be updated; Read the decomposition factor data stored in step 3 from the disk; Create a numerical update memory buffer and load the read decomposition factor data into the numerical update memory buffer; A multi-wavefront numerical update algorithm is adopted, which uses the loaded decomposition factor data to perform numerical update calculations on the subsequent wavefront matrix.

7. The sparse direct solution method based on data kernel external storage according to claim 6, characterized in that, Step 2 also includes: The overall memory requirements for numerical decomposition are estimated based on the dimension of the coefficient matrix and the non-zero element distribution pattern. Based on the estimated total memory requirements and the detected available physical memory capacity, the adjustment value of the region size parameter is calculated through an optimization algorithm; Based on the detected disk input / output performance characteristics, an optimized value for the asynchronous storage memory buffer size is calculated using a performance optimization algorithm.

8. The sparse direct solution method based on data kernel external storage according to claim 7, characterized in that, Step 2 also includes: Before performing the left-look numerical decomposition calculation in step 2, the mathematical properties of the coefficient matrix are detected by a matrix analysis algorithm; When a symmetric positive definite matrix is ​​detected, the Cholesky decomposition algorithm is used to perform numerical decomposition. When an asymmetric matrix is ​​detected, the LU decomposition algorithm is used to perform numerical decomposition. During the numerical decomposition calculation process, both real and complex number data types are processed simultaneously.

9. The sparse direct solution method based on data kernel external storage according to claim 8, characterized in that, Step 4 involves using an input / output mechanism to monitor the loop execution status, including: The status monitoring mechanism collects disk storage operation status signals through input and output; Based on the monitored input and output operation status signals, a dynamic scheduling algorithm is used to dynamically schedule the execution tasks of the producer role computing threads. The scheduling of computational tasks and disk input / output tasks is carried out in parallel, completing the overlap of computation and input / output operations.

10. The sparse direct solution method based on data kernel-external storage according to claim 9, characterized in that, Step 5 involves performing numerical calculations on the factorization data using a previous-backward substitution algorithm in the back-substitution solution memory workspace, including: Based on the solution order determined by the task dependencies, read the complete stored decomposition factor data from the disk in reverse order; Load the read decomposition factor data into the back-substitution solution memory work area; According to the calculation order required by the previous generation and back generation algorithm, the previous generation calculation and back generation calculation steps of the previous generation and back generation algorithm are executed in sequence. By merging the results of previous and back-substitution calculations, a numerical solution to a high-dimensional sparse linear system of equations is generated.

Images