CN102486727A - Parallel kraut decomposition method for ultra-large-scale matrix multi-core based on thread building blocks - Google Patents

Abstract

The invention relates to a multinuclear parallel crout decomposition method for an ultra-large scale matrix based on a TBB (Treading Building Block), which comprises the following steps: 1) rewriting a parallel computing part according to a traditional Crout method into a standard part meeting TBB requirement; 2) setting a question initial value; 3) entering into line circulation of the Crout method; 4) calling a parallel_reduce parallel module, storing the maximum pivot element in each line into a temporary vector, and meanwhile, calculating and storing a scale factor; 5) entering into row circulation of the Crout method; 6) confirming a new pivot element and the scale factor, modifying TINY, dividing each line by the pivot element; and 7) finishing the decomposition of the ultra-large scale matrix, thereby obtaining a row and line arranging sequence changed due to a local pivot element method. Compared with the prior art, the multinuclear parallel crout decomposition method provided by the invention has the advantages that the running efficiency of matrix LU decomposition is greatly increased, the method is performed on different platforms, the method is extendable, the application is wide, and the like.

Description

Parallel kraut decomposition method for ultra-large-scale matrix multi-core based on thread building blocks

technical field

The invention relates to an ultra-large-scale matrix decomposition method, in particular to a multi-core parallel Crout decomposition method of an ultra-large-scale matrix based on a thread building block.

Background technique

The Threading Building Blocks (TBB) just promoted by Intel is a multi-threaded parallel programming model based on C++, which is used to support parallel computing of multi-core processors. It has a mature data structure, supports scalable thread nesting parallelism, and supports Scalable memory allocation, support for multiple system platforms, and support for multiple compilers. TBB's programming model is to use templates as a parallel iteration model. This makes it easy for programmers to deal with issues such as synchronization, load balancing, and cache optimization, and to easily implement parallel programs that are automatically scheduled, making full use of the CPU's multi-core computing capabilities and other system resources. TBB improves the scalability and scalability of the program, and fully supports nested parallel programming. Programmers can create their own parallel components to build large parallel programs. At the same time, because TBB is developed by chip manufacturer Intel, it has better support for multi-core platforms.

TBB is an open source C++ templating library that abstracts threads into tasks and then creates reliable, portable, and scalable parallel applications. Writing task-based parallel applications using TBB helps improve the development efficiency of scalable software that runs across multi-core platforms. Compared with other threading methods such as native threads and thread wrappers, TBB is the most efficient way to execute parallel applications, and it can fully utilize the parallel performance of multi-core platforms. TBB has the following three significant advantages: 1) TBB abstracts threads and raises them to the level of tasks, which simplifies the development of parallel applications; 2) The performance of applications developed based on TBB can increase 3) TBB can reduce common thread errors such as deadlock and resource competition, and provides a cross-platform, scalable parallel solution.

At present, in theoretical research and practical application, matrix inversion has a wide range of applications, such as solving linear equations, dynamic programming, computer graphics processing, control engineering and so on. At present, there are mainly the following methods for matrix inversion: (1) adjoint matrix method; (2) elementary transformation method; (3) Gauss-Jordan elimination method; (4) triangular decomposition method. To be introduced below.

(1) Adjoint matrix method. Use the adjoint matrix of the matrix and the determinant of the matrix to find the inverse matrix of the matrix. This method is generally suitable for matrices with small dimensions

(2) Elementary transformation method. The original matrix and the identity matrix undergo elementary row (or column) transformation at the same time. When the original matrix becomes the identity matrix, the identity matrix becomes the inverse matrix of the original matrix.

(3) Gaussian-Jordan elimination method: also known as Gaussian elimination method and Gaussian elimination method, it is actually what we commonly call addition and subtraction method. Mathematically, the Gauss-Jordan elimination method, named after Gauss and Jordan, is an algorithm in linear algebra for solving the solution of a system of linear equations, solving the rank of a matrix, and solving the inverse of a reversible square matrix. When applied to a matrix, Gaussian elimination produces "row-eliminated trapezoidal form".

In the process of eliminating elements, it may happen that the pivot is zero, and the elimination method cannot be carried out at this time; even if the pivot is not zero but very small, if it is used as a divisor, the order of magnitude of other elements will increase seriously and be discarded. The diffusion of the input error also makes the calculation solution unreliable at last. In order to make the Gauss-Jordan elimination method have better numerical stability, the complete principal element elimination method can be used, but the complete principal element elimination method takes more machine time when selecting the pivot, so consider using the column pivot Elimination method, but because the pivot still needs to be selected, the amount of calculation is also relatively large. For well-conditioned problems, the Gauss-Jordan elimination method may also give very bad results, which shows that the algorithm of the Gauss-Jordan elimination method is very unstable. In fact, general matrices are ill-conditioned matrices, and the Gauss-Jordan elimination method cannot obtain satisfactory results.

(4) Triangular decomposition method

The triangular decomposition method is to decompose the square matrix into an upper triangular matrix and a lower triangular matrix, which is also called LU decomposition method. The purpose of this decomposition method is mainly to simplify the calculation of the determinant value of a large matrix, matrix inversion operation and solution of simultaneous equations. It should be noted that the upper and lower triangular matrix obtained by this decomposition method is not unique, and several pairs of different upper and lower triangular matrix pairs can be found, and their product will also obtain the original matrix. It is conditional that a matrix can be decomposed by LU, which requires it to be a square matrix and all its sequential principal subforms are not equal to zero. The main advantage of LU decomposition is that it can formulate the time-consuming elimination step, operating only on the coefficient matrix. One motivation for LU factorization is that it provides an efficient way to invert matrices and also provides a way to evaluate the state of a system of equations.

Among the above matrix inversion methods, the adjoint matrix method and the elementary transformation method are easy to understand and suitable for manual solution, but for the solution of large-scale matrices, these two methods are not suitable; Gauss-Jordan elimination method and LU decomposition method It is convenient for the implementation of the program, and under most conditions, it can adapt to the solution of large-scale matrices. However, for large-scale or ultra-large-scale matrices, the Gauss-Jordan elimination method and the traditional serial LU decomposition method often require too much calculation, which requires high performance of the computer, and the real-time performance of the algorithm is not high. Applications with strict requirements are not suitable.

Contents of the invention

The purpose of the present invention is to provide a super-large-scale matrix multi-core parallelism based on thread building blocks that can greatly improve the operating efficiency of matrix LU decomposition in order to overcome the defects in the above-mentioned prior art, and can be cross-platform, scalable, and widely used. Kraut decomposition method.

The object of the present invention can be realized through the following technical solutions: the ultra-large-scale matrix multi-core parallel kraut decomposition method based on thread building block, it is characterized in that, this decomposition method comprises the following steps: 1) installation and environment setting of TBB parallel computing platform ;2) Analyze and extract the part that can be calculated in parallel in the traditional Crout method, and use the TBB parallel template class to rewrite it as a specification class that meets the needs of TBB; 3) Set the initial value of the problem; 4) Enter the row loop of the Crout method; 5) Call parallel_reduce Parallel module class, save the maximum pivot of each row in a temporary vector, and calculate and save the scale factor at the same time; 6) Enter the column loop of the Crout method; 7) Determine the new pivot and scale factor, modify TINY, divide each row Use the pivot to judge whether the column cycle is over, if it is judged otherwise, go to step 6), judge whether the row cycle is over, if it is judged otherwise, go to step 4); The order of rows and columns changed by the pivot method.

The initial value of the problem in the step 3) includes the scale of the matrix, the number of rows and columns of the matrix, and the input matrix.

Compared with the prior art, the present invention has the advantages of high efficiency, cross-platform, scalability, and wide application, and the parallel Crout method of the present invention can greatly improve the operating efficiency of matrix LU decomposition; on this basis, the present invention can The efficiency of matrix inversion is greatly improved, and the effectiveness and high efficiency of the present invention are further verified. The applicable conditions of the parallel Crout method are: if all the subforms of the matrix are non-zero, the invention can be used for LU decomposition, and then the inverse of the matrix or the solution of the linear equation system can be obtained. When the scale of the matrix reaches one million levels, the invention improves about 13.7% than the serial Crout method and about 77% than the Gauss-Jordan elimination method. Therefore, the present invention has very large practical value and wide application prospect.

Description of drawings

Fig. 1 satisfies the TBB of the present invention and seeks the maximum pivot class structural diagram;

Fig. 2 is the parallel Crout method flowchart of the present invention;

Fig. 3 is to utilize parallel Crout method of the present invention to ask the flow chart of matrix inverse;

Fig. 4 is a schematic diagram of the result of solving the matrix inverse by using the present invention.

Detailed ways

The present invention will be described in detail below in conjunction with the accompanying drawings and specific embodiments.

Example

TBB defines the concept of tasks. When initializing TBB task scheduling, the task scheduler object task_scheduler_init implements multi-task allocation and parallel computing, and supports multi-thread division. When calling the template class for parallel computing, the value range of the loop processing and the task granularity parameters are specified by the template class parameters. The task granularity parameter determines the granularity of task division. If the granularity is too large, the operating efficiency cannot be fully improved; if the granularity is too small, the overhead caused by excessive parallel task allocation will reduce the operating efficiency. If the appropriate task granularity cannot be obtained, the automatic allocation function auto_partitioner() provided by TBB can be used to help users set appropriate task granularity parameters.

Taking the maximum principal element template class designed by parallel_reduce as an example, the implementation process of the parallel maximum principal element algorithm designed by TBB is described in detail for the Crout method to call. First, based on the parallel_reduce template, this process is encapsulated into a class MaxElement, whose function is to solve the largest principal element. The interfaces that must be included in class MaxElement are operator and constructor, the specific form is as follows:

class MaxElement{

...

void operator()(const blocked_range<size_t>&r){

void operator()(const blocked_range<size_t>&r){

for(size_t i=r.begin(); i!=r.end();++i){

...

}}}

//Create a branch, if you want to access x, you should guarantee atomic operation

MaxElement(MaxElement&x, tbb::split):

my_a(x.my_a), value_of_max(FLT_MIN), index_of_max(-1){...}

//merge branches

void join(const MaxElement&y){

...

}}

//Constructor

MaxElement(DOUBLE*a):

my_a(a), value_of_max(FLT_MIN), index_of_max(-1){}

};

Among them, the operator interface is the main part of parallel processing. Its main function is to optimize the parallel loop, modify the parameters of the loop to the blocked_range template class defined by TBB, and support the parallel division of tasks in the loop body. In the operator interface, the comparison size operation in each small block is performed in parallel, and the maximum value in the small block is returned. When TBB decides to split a blocked_range, it will call MaxElement(MaxElement&x, tbb::split). This tbb::split is just a point symbol, which is used to distinguish it from the copy constructor. When TBB merges blocked_range, it will call voidjoin(MaxElement&y). In this example, merging means operating multiple results, saving and returning the final result. The constructor mainly implements parameter initialization and separates and constructs subtasks from the entire task space. After completing the writing of the parallel template class, initialize the TBB task scheduler, call the parallel template written above, return the calculation result, and finally end the TBB task scheduling.

The basic steps of using TBB to design a parallel Crout method are as follows:

Step 1: Installation and environment setting of TBB parallel computing platform;

Step 2: Analyze and extract the part that can be calculated in parallel in the traditional Crout method, and use the TBB parallel template class to rewrite it into a standard class that meets the needs of TBB;

Step 3: Set the initial value of the problem, that is, the size of the matrix, the number of rows and columns of the matrix, the input matrix, etc.;

Step 4: Enter the line loop of the Crout method;

Step 5: Call the parallel_reduce parallel module class, save the maximum pivot of each row in a temporary vector, and calculate and save the scale factor at the same time, so using TBB for parallel design is very important to improve the efficiency of the method;

Step 6: Enter the column loop of the Crout method;

Step 7: Determine the new pivot and scale factor, modify TINY, and divide each row by the pivot. Go to step 6 if the column cycle is not over, judge whether the row cycle is over, if judged otherwise go to step 4); otherwise go to step 4;

Step 8: At the end, get the order of ranks and columns changed by the partial pivot method.

In order to verify the correctness, effectiveness and evaluation of the present invention's practicability of the parallel Crout method, the inverse of the solution matrix of the present invention is verified and evaluated below, and its basic steps are as follows:

Step 1: Installation and environment setting of TBB parallel computing platform;

Step 2: Read in the file, including the input matrix and related parameter settings;

Step 3: Set the initial value of the problem, that is, the number of rows and columns of the matrix, and input the matrix;

Step 4: Utilize the parallel Crout method of the present invention to carry out LU decomposition to the input matrix;

Step 5: Enter the row loop of the matrix;

Step 6: Enter the column loop of the matrix;

Step 7: Invert the matrix column by column through the back substitution function;

Step 8: Increase the column count, if the column cycle is not over, go to step 6; increase the row count, if the row cycle is not over, go to step 5;

Step 9: End, get the inverse matrix of the input matrix and write it to the file.

As shown in Figures 1 and 2, the present invention analyzes the traditional Crout decomposition method, parallelizes the key steps, and utilizes the thread building blocks of Intel Corporation to realize parallelization. The present invention makes full use of the flexibility, cross-platform and stability of the thread building block, realizes the parallel Crout decomposition method based on the thread building block, and applies the method to matrix inversion. A large number of simulation experiments show that the present invention can greatly The efficiency of the traditional Crout decomposition method is greatly improved, and then the efficiency of the matrix inversion method using the present invention is greatly improved.

The present invention is composed of a maximum pivot class, a maximum pivot call class and a Crout decomposition method using the maximum pivot call class that conform to the thread building block specification, and the maximum pivot class that meets the thread construction block specification provides specific Part of parallel processing, task thread branch division and merging and other functions; the maximum main element class call class provides the call encapsulation of the largest main element class, and provides the granularity control parameters for thread division; parallelized Crout The decomposition method takes advantage of thread parallelization, makes full use of computer resources, and improves the efficiency of the traditional Crout decomposition method. The present invention can run on multiple system platforms, such as Windows, Linux, Unix and other systems.

Fig. 1 seeks the structural diagram of maximum pivot class for the present invention, the detailed description of each part in the figure below:

In part 101, the operator interface is implemented, which is the main part of parallel processing, and its main function is to perform parallel loop optimization, and modify the parameters of the loop to the blocked_range template class defined by TBB, which can support the parallel division of tasks in the loop body;

In part 102, the branch creation function interface with tbb::split is implemented. When TBB decides to split a blocked_range, it will call MaxElement(MaxElement&x, tbb::split). This tbb::split is just a point symbol. Use It is different from the copy constructor;

In part 103, implement the join merge function interface. When TBB merges blocked_range, it will call void_join(MaxElement&y). In this example, merge means to operate multiple results, save and return the final result;

In part 104, the constructor interface is implemented, and the constructor mainly implements parameter initialization and separates and constructs subtasks from the entire task space.

Fig. 2 is the parallel Crout method flowchart of the present invention, each step in the figure is described in detail below:

In step 201, read in the matrix to be decomposed, initialize the row value i, the column value j, the vector vv[i] containing the scale factor, the very small number TINY, temporary variables and so on. Then execute step 202;

In step 202, loop by row, increment row, and find the contained scale factor. Then execute step 203;

In step 203, a parallel algorithm for finding the maximum principal element is invoked: the principal element and the scale factor of each row are calculated. Then execute step 204;

In step 204, loop by column, incrementing column. Then execute step 205;

In step 205, new candidate principal elements and exchange scale factors are calculated to replace the original candidate principal elements and exchange scale factors. Then execute step 206;

In step 206, it is judged whether the matrix is a singular matrix, if yes, TINY is assigned zero, otherwise TINY remains unchanged. Then execute step 207;

In step 207, each row in the matrix is divided by the principal element. If the circulation by column is not finished, then perform step 204; if the circulation by column is finished but the circulation by row is not finished, then perform step 202; if the circulation by column and the circulation by row have all ended, then perform step 208;

In step 208, the end function: output the arrangement order of the rows changed by the partial pivot method, and the LU decomposition matrix of the matrix can be obtained according to the arrangement order.

Fig. 3 is the flow chart of using the parallel Crout method to find the inverse matrix of the matrix, and each step in the figure is described in detail below:

In step 301, the required input file is selected through the dialog box, the matrix file selected is in TXT format, and the data structure format is: the row number (row) of the first row matrix; the column number (col) of the second row matrix, in In this application, the number of rows is equal to the number of columns (row=col=n), that is, the matrix here is an n-dimensional square matrix; the third row starts with the value of the matrix. Then execute step 302;

In step 302, the required values are read according to the dimensions of the read-in matrix, and related variables and data structures are initialized. Then execute step 303;

In step 303, the read-in matrix is decomposed by using the parallel Crout method of the present invention. Then execute step 304;

In step 304, the row loop i of the control matrix is executed in step 305;

In step 305, the column cycle j of the control matrix, and then perform step 306;

In step 306, the calculated matrix is processed through a back substitution function. If i<n and j<n, then return to step 305; if i<n and j>=n, then return to step 304; if i>=n and j>=n, execute step 307;

In step 307, the matrix inversion calculation is ended, and the inverse of the matrix and related calculation parameters are output: the name of the method, the operation time, the number of method threads, the size of the method sub-matrix, and the like. Execute step 308;

In step 308, the result is written into a TXT file for further analysis by the user.

Figure 4 compares the efficiency of the parallel crout method, serial crout method, and Gauss-Jordan elimination method in solving the matrix in detail. It can be seen from the reference line that when the matrix dimension is greater than or equal to 256, the running time of the parallel crout method is better than that of the serial crout method, and is significantly better than that of the Gauss-Jordan elimination method; and as the matrix dimension increases increases, this trend becomes more pronounced. When the scale of the matrix reaches one million levels, the invention improves about 13.7% than the serial Crout method and about 77% than the Gauss-Jordan elimination method. That is, for large-scale or very large-scale matrices, the parallel crout method is obviously better than the Gauss-Jordan elimination method, and it is also better than the serial crout method. Therefore, the present invention has high practical value and wide application prospect.

Claims (2)

1. based on the ultra-large matrix multi-core parallel concurrent Ke Laote decomposition method of thread constructing piece, it is characterized in that this decomposition method may further comprise the steps:

1) installation and the environment setting of thread constructing piece (TBB) parallel computing platform;

2) analyze and extract in traditional Ke Laote (Crout) method can parallel computation part, utilize TBB parallel templates class to be rewritten as the standard class that meets the TBB needs;

3) the problem initial value is set;

4) getting into going of Crout method circulates;

5) call the parallel module class of parallel_reduce, each maximum pivot of going is kept in the interim vector, calculate and preserve scale factor simultaneously;

6) the row circulation of entering Crout method;

7) confirm new pivot and scale factor, revise TINY, each row judges divided by pivot whether the row circulation finishes, if be judged as otherwise forward step 6) to, judges whether the row circulation finishes, if be judged as otherwise forward step 4) to;

8) end of run is accomplished ultra-large matrix decomposition, obtains because of the altered ranks ordering of partial pivot method.

2. the ultra-large matrix multi-core parallel concurrent Ke Laote decomposition method based on the thread constructing piece according to claim 1 is characterized in that the problem initial value in the described step 3) comprises scale, the line number of matrix, the columns of matrix, input matrix.

Images