Patent Publication Number: US-8527569-B2

Title: Parallel processing method of tridiagonalization of real symmetric matrix for shared memory scalar parallel computer

Description:
CROSS-REFERENCE TO RELATED APPLICATION 
     This application is a continuation application of International PCT Application No. PCT/JP2007/000424, filed on Apr. 19, 2007, the entire contents of which are incorporated herein by reference. 
    
    
     FIELD 
     The embodiments discussed herein are related to a parallel processing method that enables high-speed processing of tridiagonalization of real symmetric matrices in shared memory scalar parallel computers. 
     BACKGROUND 
     As a method of solving eigenvalue problems of real symmetric matrices, there is a method in which tridiagonalization is performed so that the eigenvalue problems of real symmetric matrices are converted into eigenvalue problems of tridiagonal matrices. This method is costly in calculations for tridiagonalizing real symmetric matrices. Shared memory scalar parallel computers, which are less capable of accessing memory than vector computers, have improved algorithms in order to increase the amount of computation processed in units of memory accesses so that tridiagonalization is facilitated. 
     Patent Document 1 describes fundamental processing for block tridiagonalization. 
     In this method, however, the greatest part of the calculation cost consists of calculations of matrix vector products, and also matrix vector products greatly influence the memory access speed. Thus, this part of the calculation cost needs to be modified in order to improve the entire performance. 
     In a blocking method by which the computation amount is increased, it is important to calculate matrix vector products by making consecutive accesses and to perform those calculations parallelly and evenly by using respective CPUs. For this purpose, the updating of the entire matrix consists of the updating of the lower triangular portion and the copying of the upper triangular portion. In order to distribute the updating loads evenly, a matrix is divided into double of the number of CPUs in the column direction, and the computation amounts for updating are assigned to the respective CPUs in pairs including the i-th CPU and the “2×#CPU-(i−1)-th” CPU. However, this method is disadvantageous in that the memory areas to be referred to by the calculations of matrix vector products are distant so that the data in such areas interfere with data in cache, resulting in difficulty in storing the data in cache. 
     Finding a method that can harmonize the updating based on matrix vector products/matrix products with the load distribution and that can enable high-speed operation of matrix vector products is expected to lead to great improvement in performance. 
     Non Patent Document 1 and Non Patent Document 2 disclose a fundamental algorithm for tridiagonalization and parallel processing of tridiagonalization, respectively.
     Patent Document 1:   Japanese Laid-open Patent Publication No. 2004-5528   Non Patent Document 1   G. H. Golub, C. F. van Loan, Matrix Computation Second Edition, Johns Hopkins University Press 1989   Non Patent Document 2   J. Choi, J. J. Dongarra, and D. W. Walker, “THE DESIGN OF A PARALLEL DENSE LINEAR ALGEBRA SOFTWARE LIBRARY: REDUCTION TO HESSENBERG, TRIDIAGONAL, AND BIDIAGONAL FORM”, Engineering Physics and Mathematics Division, Mathematical Sciences Section, prepared by the Oak Ridge National Laboratory managed by Martin Marietta Energy System, Inc., for the U.S. DEPARTMENT OF ENERGY under Contract No. DE-AC05-84OR21400, ORNL/TM-12472.   

     SUMMARY 
     A parallel processing method according to the present invention is a parallel processing method for performing high-speed processing of tridiagonalization of a real symmetric matrix in a shared memory scalar parallel computer including a plurality of processors, the method including dividing a lower triangular matrix of the real symmetric matrix into as many column block matrices as there are processors, and loading the column block matrix onto a storage area in each of the processors; and making each of the processors perform computation in a longitudinal direction and a transverse direction of the real symmetric matrix for an element of a column block matrix read in a single loading operation. 
     The object and advantages of the invention will be realized and attained by means of the elements and combinations particularly pointed out in the claims. 
     It is to be understood that both the foregoing general description and the following detailed description are exemplary and explanatory and are not restrictive of the invention, as claimed 
    
    
     
       BRIEF DESCRIPTION OF DRAWINGS 
         FIG. 1  is a block diagram illustrating a hardware configuration of a shared memory parallel computer on which the present invention is based; 
         FIG. 2  illustrates processing of block tridiagonalization (first); 
         FIG. 3  illustrates processing of block tridiagonalization (second); 
         FIG. 4A  illustrates processing of block tridiagonalization (third); 
         FIG. 4B  illustrates processing of block tridiagonalization (fourth); 
         FIG. 4C  illustrates processing of block tridiagonalization (fifth); 
         FIG. 5  illustrates a parallel processing method according to an embodiment of the present invention (first); 
         FIG. 6  illustrates a parallel processing method according to an embodiment of the present invention (second); 
         FIG. 7  illustrates a parallel processing method according to an embodiment of the present invention (third); 
         FIG. 8  illustrates, in more detail, a parallel processing method according to an embodiment of the present invention (first); 
         FIG. 9  illustrates, in more detail, a parallel processing method according to an embodiment of the present invention (second); 
         FIG. 10  illustrates, in more detail, a parallel processing method according to an embodiment of the present invention (third); 
         FIG. 11  illustrates, in more detail, a parallel processing method according to an embodiment of the present invention (fourth); 
         FIG. 12  illustrates a flowchart for a process according to an embodiment of the present invention; 
         FIG. 13  illustrates a flowchart for a process according to an embodiment of the present invention; 
         FIG. 14  illustrates a flowchart for a process according to an embodiment of the present invention; 
         FIG. 15  illustrates a flowchart for a process according to an embodiment of the present invention; 
         FIG. 16  illustrates a flowchart for a process according to an embodiment of the present invention; 
         FIG. 17  illustrates a flowchart for a process according to an embodiment of the present invention; 
         FIG. 18  illustrates a flowchart for a process according to an embodiment of the present invention; 
         FIG. 19  illustrates a flowchart for a process according to an embodiment of the present invention; and 
         FIG. 20  illustrates a flowchart for a process according to an embodiment of the present invention. 
     
    
    
     DESCRIPTION OF EMBODIMENTS 
     The present invention is based on a result of modifying the method of tridiagonalization in Patent Document 1. 
     An embodiment of the present invention uses a blocked algorithm so that a portion to be updated can be calculated highly densely when performing tridiagonalization. The symmetry in the matrix vector product portion is utilized in order to parallelly perform the operation of matrix vector products only of the lower triangular portion. In this calculation, the matrix vector product appears twice; however, reference is only made to memory once, as a shared access. Because only the lower triangular portion is used, the lower triangular matrix is divided dynamically so that the computation amount is distributed evenly. This eliminates the need for redundant copying, and also makes it easier to store, in cache, the content of memory referred to by the matrix vector products. 
     A matrix is divided in the column direction so that not only is the computation amount for updating by the matrix product of the lower triangular portion (rank-p update) uniform but the computation amount of the matrix vector products that refers to the matrix is also uniform (the areas of the results of dividing the lower triangular portion are uniform). In this state, each of the areas obtained by the division can be treated as an area consisting of two portions, a lower triangular portion and a rectangular portion. In the updating of a lower triangular portion, this portion is divided appropriately, and the resultant pieces are updated using several matrix products. When divided lower triangular portions including diagonal elements are updated using normal matrix products, the upper triangular portion will also be updated partially. However, the computation amount is reduced by reducing the computation of extra portions to the minimum possible amount so that the updating is performed with a computation amount close to that in the updating of only the lower triangular portion. 
       FIG. 1  is a block diagram illustrating a hardware configuration of a shared memory parallel computer on which the present invention is based. 
     As illustrated in  FIG. 1 , the parallel computer on which the present invention is based has CPUs  12 - 1  through  12 - m  and memory modules  10 - 1  through  10 - n  connected to each other via an interconnection network  11  that includes a bus. The CPUs  12 - 1  through  12 - m  each include an L2 cache-and-bus interface module, and a plurality of L1 cache-and-CPU core modules connected to the L2 cache-and-bus interface module. This configuration is called a multi-core CPU. A multi-core CPU is constituted of a plurality of CPU cores (including L1 caches) embedded in a CPU, and each of the cores shares the L2 cache-and-bus interface. Each CPU core operates logically as if it were a CPU. In the example in  FIG. 1 , two CPU cores are embedded in one CPU; however, four CPU cores are embedded in one CPU in some cases. 
       FIGS. 2 through 4C  explain block tridiagonalization computation. 
     An explanation will be given for the mathematical framework of tridiagonalization for obtaining an eigenvalue/eigenvector by referring to  FIGS. 2 through 4C . 
     a) Mathematical Algorithm for Block Tridiagonalization 
     A matrix is tridiagonalized in units of the block width. A matrix is divided into blocks, and tridiagonalization is performed using the following block algorithm in units of the blocks.  FIG. 2  illustrates the processing of the n-th block. 
     The process for the last block is terminated by applying the algorithm to the leftmost block, whose block width is 2. 
     Step 1 Generate Householder vector u from the row vector of n+1 of A n    
     Step 2 calculate v i =A n+1 u, w i =v i −u (u t  v)/2 
     Step 3 update U i  and W i  so that U i =(U i−1 , u i ) W i =(W i−1 , w i ) 
     Step 4 if (i&lt;blks) then update “n+i+1”-th row of A n ; 
     A n (*,n+i+1)=A n (*,n+i+1)−U i W i (n+i+1,*) t −W i U i  (n+i+1,*) t    
     endif 
     enddo 
     step 5: A n+blks =A n −U blks  W blks   t −W blks  U blks   t    
     Next, an explanation is given for tridiagonalization based on the blocked Householder transformation. Blocked householder transformation is performed as described below.
 
 v =( v   1   ,v   2   , . . . ,v   n )
 
| v|   2   =v·v=h   2 ,
 
When U v =(h, 0, . . . , 0) t , U v =v−(v 1 −h, v 2 , . . . , vn) is satisfied. U=(I−uu t /|u| 2 )=(I−αuu t ) is satisfied, where u=(v 1 −h, v 2 , . . . , v n )
 
In the following calculation, α is omitted.
 
                     A     n   +   1       =       ⁢       U   t     ⁢     A   n     ⁢   U                 =       ⁢       (     I   -     uu   t       )     ⁢     A   ⁡     (     I   -     uu   t       )                     =       ⁢       A   n     -       uu   t     ⁢     A   n       -       A   n     ⁢     uu   t       +       uu   t     ⁢     A   n     ⁢     uu   t                     =       ⁢       A   n     -     uv   t     -     vu   t     +       uv   t     ⁢     uu   t                     =       ⁢       A   n     -     uw   t     -       uu   t     ⁢     u   t     ⁢     v   /   2       -     wu   t     -       uu   t     ⁢     u   t     ⁢     v   /   2       +       uu   t     ⁢     u   t     ⁢   v                   =       ⁢       A   n     -     uw   t     -       wu   t     ⁢           ⁢   …   ⁢           ⁢     (*   )                     
where w=v−u(u t  v)/2, v=A n  u
 
When this is repeatedly used,
 
A n+k =A n −U k  W k   t −W k  U k   t  . . . (**) is satisfied.
 
In the calculation in the k-th step, V n  can be calculated from (*) and (**) as follows:
 
 v   k   =A   n   u   k   −U   k−1   W   k−1   t   u   k   −W   k−1   t   u   k  
 
 w   k   =v   k   −u   k   u   k   t   v   k /2
 
 U   k =( U   k−1   ,u   k ), W   k =( W   k−1   ,w   k )
 
 A   n+k   =A   n   =U   k   W   k   t   −W   k   U   k   t  
 
b) Storing Information That Constitutes A Householder Transformation
 
     When an eigenvector is calculated, a Householder transformation used for the tridiagonalization is needed. Thus, Un and α are stored at the positions of the vector that constitute the Householder transformation. α is stored at the position of the corresponding diagonal element. 
     c) Method of Obtaining U i  Effectively 
     In order to perform Householder transformation for the tridiagonalization of the blocked part, the next vector has to be updated. In order to perform this calculation as locally as possible, the working area of the block width part is copied to perform tridiagonalization, and the result of the tridiagonalization is stored in the original area. The calculation is performed with increased computation density in the form of a matrix product, and the next column vector is not updated. For this reason, the blocked part is tridiagonalized using a recursive program. 
     Specifically, the block area of matrix A, at the left side of  FIG. 3 , is copied onto a working area, and the block area is divided into the first update portion and the last update portion to process them recursively.  FIGS. 4A through 4C  illustrate the recursive process. When it is assumed that the shaded area is updated in  FIG. 4A , the shaded area in  FIG. 4B  is updated next, and thereafter the shaded area in  FIG. 4C  is updated. In this manner, the first and last portions are updated in a nested manner. Specifically, when a recursive program is called to reach depth  2 , the shaded area in  FIG. 4A  is updated as B in the first process, then the shaded area in  FIG. 4B  is updated, and thereafter the shaded area in  FIG. 4C  is updated. In the parallelization for the updating, the block matrix of the update part is divided uniformly in the vector direction, and the areas obtained by uniformly dividing each part by the processor are updated parallelly. 
       FIGS. 5 through 7  illustrate a parallel processing method according to an embodiment of the present invention. 
     1) A method in which computation is performed twice on a matrix vector product using only the lower triangular portion and through one access (loading) 
     As illustrated in  FIG. 5 , a matrix vector product is calculated using only the lower triangular portion. The size of a matrix used for the matrix vector product is assumed to be N n ×N n . When this is used in a 4-parallel configuration, the positions at which the lower triangular portion is divided into four portions are obtained in the following way.
 
 i   m   =N   n ×(1−sqrt(1− m/#p ))
 
     The results obtained under m=1, 2, and 3 correspond to the boundary positions. In the above expression, #p represents the total number used for the division, i.e., the parallel number. 
     By utilizing the symmetry in A n , i.e., A n =A n   t , a matrix product is calculated from the information in the area to be processed by each processer. 
     Each portion includes lower triangular matrix portion D m  of the diagonal block matrix portion and rectangular matrix R m  below lower triangular matrix portion D m . 
     The m-th portion, which is a result of dividing the lower triangular matrix into #P, is transposed, and the resultant value corresponds to the result of dividing the upper triangular portion in the row direction. According to this relationship, the calculations below are satisfied.
 
 v   xm1 ( i   m +1: N   n )= R   m   ×u   1 ( i   m−1 +1 :i   m )  (expression 1)
 
 v   xm2 ( i   m−1 1: i   m )= R   m   t   ×u   1 ( i   m +1 :N   n )  (expression 2)
 
 v   xm3 ( i   m−1 +1 :i   m )= D   m   ×u   1 ( i   m−1 +1 :i   m )  (expression 3)
 
 v   xm4 ( i   m−1 +1 :i   m )=notdiag( D   m ) t   ×u   1 ( i   m−1 +2 :i   m )  (expression 4)
 
(notdiag(X) represents the lower triangular matrix obtained by removing the diagonal portion from lower triangular matrix X)
 
     
       
         
           
             
               
                 
                   
                     v 
                     1 
                   
                   = 
                   
                     
                       ∑ 
                       
                         m 
                         = 
                         1 
                       
                       
                         # 
                         ⁢ 
                         p 
                       
                     
                     ⁢ 
                     
                       ( 
                       
                         
                           v 
                           
                             xm 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             1 
                           
                         
                         + 
                         
                           v 
                           
                             xm 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             2 
                           
                         
                         + 
                         
                           v 
                           
                             xm 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             3 
                           
                         
                         + 
                         
                           v 
                           
                             xm 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             4 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   [ 
                   
                     Expression 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     1 
                   
                   ] 
                 
               
             
           
         
       
     
     In this calculation, the calculations of v xm1  and v xm2  are performed simultaneously in an access to one element while accessing the elements in the matrix successively, and thereby it is possible to double the performance of the matrix vector product. 
     An example of a program executed simultaneously using expressions 1 and 2 is described below. 
     For the sake of simplicity, it is assumed that R m  is a matrix n×m, x 1  is a vector of length of n, x 2  is a vector of length m, y 1  is a vector of length m, and y 2  is a vector of length n. 
     x 1 =R×y 1  and x 2 =R t ×y 2  are calculated commonly so that R is referred to.
     x1(1:n)=0   do j=1,m   tmp1=y1(j)   sum=0   do i=1,n   tmp2=y2(i)   tmp=R(i,j)   x1(i)=x1(i)+tmp*tmp1   sum=sum+tmp*tmp2   enddo   x2(j)=sum   enddo   

     In this case, x1 and y2 are small in data amount, and are accessed repeatedly, and accordingly they are cached so that they can be accessed at high speed. In comparison to accesses to R(i.j), the computation of x1 and x2 is so fast that it can be ignored, and thus the performance becomes twice the performance that would be obtained when performing the computation of matrix vector products separately. 
     The equation for calculating the boundary of the division, i m =N n ×(1−qsrt(1−m/#p)) can be obtained in the following manner. 
     Each of the CPUs is assigned the same number of computations. Thus, it is assumed that the areas corresponding to the matrices that are proportional to the number of elements (or the number of elements itself) are identical. 
     The calculation is performed on the assumption that N n   2 −(N n −i m ) 2 =N n   2 ×m/#p and r=i m /N n . 
     r 2 −2r=m/#p and 0&lt;=r&lt;=1, where r is the square root (r−1) 2 =m/#p
     r=1−sqrt(1−m/#p) (the smaller is employed in the range of the square root)
 
2) Method Of Parallelizing Updated Portion and Calculating Diagonal Block Portion
   

     The lower portion is calculated using the symmetry of the updated portion in step 5. A matrix is divided so that each CPU (thread) is assigned a uniform amount of the computation for the update. The division is performed in the same manner as in the division using a matrix vector product. 
     The size of a matrix used for the matrix vector product is assumed to be N n ×N n . When a four-parallel configuration is adopted, the positions at which the lower triangular portion is divided into four portions are obtained in the following way:
 
 i   m   =N   n ×(1−sqrt(1 −m/#p ))
 
     The results with m=1, 2, and 3 correspond to the boundary positions. In the above expression, #p represents the total number used for the division, i.e., the parallel number. 
     As has already been described, each of the areas obtained by the division in units of CPUs includes the diagonal block matrix portion and the rectangular matrix portion.  FIG. 6  illustrates matrix A n+I  divided in units of CPUs. The portion corresponding to the rectangular matrix portion is calculated using a normal matrix product, as in the expression below.
 
 A   n+k   =A   n   −U   k   W   k   t   −W   k   U   k   t  
 
     The updating of the diagonal block matrix portion is performed as illustrated in  FIG. 7 . Specifically, a point is determined that divides the diagonal line into two lines of equal length for the lower triangular portion of the diagonal block. The lower triangular matrix is divided in the row and column directions from this point. The rectangular portion is updated by a matrix product. Two half-sized lower triangular matrices are generated, and a point dividing the diagonal line into two lines of equal length is similarly determined for each of these portions, and the division in the row and column directions are performed similarly. When the target lower triangular matrix has become sufficiently small, the upper triangular matrix is calculated together by the matrix product. Thereby, a simple recursive program can be produced. 
     The above parallel processing can implement tridiagonalization of real symmetric matrices at high speed and with scalability. The above method has proven to be twice the conventional methods in performance. 
       FIGS. 8 through 11  illustrate the parallel processing method according to an embodiment of the present invention in more detail. 
     In block tridiagonalization computation, the following processes are performed. 
     1) Tridiagonalization is performed in units of blocks. The upper block is tridiagonalized, and matrices U and W resulting from the tridiagonalization are used to update the square matrix that has been reduced by the distance of the block width (the lower triangular matrix is illustrated by the thick dotted lines in  FIG. 8 ). However, this square matrix is a symmetric matrix, and therefore only the lower triangular matrix part (the triangle enclosed by the thick dotted lines) is updated to be used as the data of the upper triangular portion.
 
2) Tridiagonalization in units of blocks requires the calculation of a matrix vector product for a square matrix (An+1: the square matrix enclosed by the thick lines) represented by the thick lines in the block. This calculation continuously reduces the size of the square matrix by 1. This portion is also symmetric, and the lower triangular portion is used for the calculation.
 
     The above calculation requires a great computation cost to perform the entire calculation. In order to cause a high effect in parallel processing, the following points have to be taken into consideration. 
     # Most of the time taken for calculation of matrix vector products in a scalar computer allowing only slow memory accesses is used for memory accesses. Thus, it is important to reduce the cost of memory accesses. 
     # When a lower triangular matrix portion is used using the symmetry of matrices, the loads are desirably distributed evenly in view of parallel processing. 
     Thus, a program in a configuration as below is prepared. 
       FIG. 9  is a functional block diagram of a program according to an embodiment of the present invention. 
     A block tridiagonalization control unit  15  includes a block portion tridiagonalization unit  16  and an updating unit  19 . The block portion tridiagonalization unit  16  includes a load control unit  17  for a matrix vector product that determines the interval of the second dimension of the lower triangular matrix treated by each thread. The load control unit  17  includes a matrix vector product parallel calculation unit  18  for performing calculation for the lower triangular portion and the rectangular matrix of a block diagonal matrix. The updating unit  19  includes a load control unit for the calculation performed by the matrix product that determines the interval of the second dimension of the lower triangular matrix treated by each thread. The load control unit  20  further includes a block updating unit  21  for updating the portions other than the block diagonal portion, and an update control unit  22  for the lower triangular matrix portion of the block diagonal unit, and the update control unit  22  controls the division of the blocks. The update control unit  22  includes a calculation unit  23  for calculating the lower triangular matrix using the approximation using fine square matrices obtained by dividing a block diagonal matrix. 
       FIG. 10  illustrates a process of calculating a matrix vector product from a lower triangular matrix (that becomes smaller by 1 each time the process is performed). 
     In this example, the calculation of V(nptr:n)=A(nptr:n, nptr:n)*U(nptr:n) is performed, where nptr represents the leftmost element in the lower triangular matrix, n represents the rightmost element in the matrix, ms represents the point of origin treated by each thread, and me represents the value indicating the point of termination treated by each thread. 
     V(*, i), U(*,i), and vsum (1:n, numthrd), and entire matrix A having the attribute “shared”, which are arrays are secured, and the following processes are performed. 
     1) vsum is cleared to zero in each thread so that vsum(nptr:n,nothrd)=0. 
     2) The above two calculations are simultaneously performed in the column direction of L2 by making reference as the arrow inside L2 in  FIG. 10  shows. See the flowchart of matvec1, which will be explained later.
 
 v sum( me+ 1 :n )=L2 *U ( ms:me )
 
 v sum( ms:me )= L 2 t   *U ( me+ 1 :n )
 
3) Calculate block diagonal portion using L1 vsum(ms:me)=vsum(ms:me)+L1*U(ms:me)+ND(L1) t *U (ms:me,i) (ND (L1) is the lower triangular matrix obtained by setting 0 in the diagonal elements of L1)
 
4) Generate an interval in which nptr:n is divided evenly in units of threads. The point of origin is “is”, and the point of termination is “ie”. The termination of the calculation is confirmed by the barrier synchronization, and parallel processing is performed in each thread.
 
V(is:ie,i)=Σsum(is:ie,nothrd)
 
nothrd=1˜numthrd
 
(nothrd ranges between thread number 1 and numthrd, numthrd being the total number of threads)
 
       FIG. 11  illustrates processes performed by the updating unit. 
     The lower triangular matrix portion of A(nptr:n, nptr:n) is updated using the following expression.
 
 A ( nptr:n,nptr:n )= A ( nptr:n,nptr:n )− U ( nptr:n, 1:blk)* W ( nptr:n ,1:bl k) t   −W ( nptr:n, 1:blk)* U ( nptr:n, 1:blk) t  
 
     In the above expression, nptr represents the leftmost element in the lower triangular matrix, blk represents the block width, n represents the rightmost element in the matrix, ms represents the point of origin treated by each thread, and me represents the point of termination treated by each thread. 
     L1 is divided into L11, L12, and L22 using the midpoints set on the perpendicular lines and the base. L11 and L22 are similarly divided using the midpoints to divide the lower triangular matrix into a group of square matrices, and this process is repeated. 
     When the block diagonal matrix has become sufficiently small, the calculation is performed using matrix products of normal matrices. This is implemented by a recursive program (see the flowchart for subroutine ltgemmtrt for detailed information). 
       FIGS. 12 through 20  illustrate the flowcharts for processes according to an embodiment of the present invention. 
       FIG. 12  illustrates a flowchart for subroutine trid for tridiagonalizing a real symmetric matrix. In this subroutine, shared arrays A(k, n), diad(n), and sdiag(n) are input. diag and sdiag are arrays used for returning the diagonal element and the subdiagonal element of the calculated tridiagonalization matrix. 
     In step S 10 , the working areas U(n+1, iblk) and v (n+1, iblk) are secured in a routine to be used with the attribute “shared”. In step S 11 , a thread is generated. In each thread, the total number of threads is set in the local area, and the thread number (numerical order) assigned to each thread is set in nothrd, and the block width is set in iblk. Also, setting is performed so that nb=(n−2+iblk−1)/iblk, nbase=0, and i=1. In step S 12 , it is determined whether or not i&gt;nb−1 is satisfied. When the determination result in step S 12  is Yes, the process proceeds to step S 18 . When the determination result in step S 12  is No, setting is performed so that nbase=(i−1)*iblk, istart=1, and nwidth=iblk in step S 13 . In step S 14 , the target area of block tridiagonalization is copied onto working area U, i.e., U(nbase+1:n, 1:iblk)←A(nbase+1:n,nbase+1:nbase+iblk). In step S 15 , subroutine blktrid is called, and the portion copied onto U is tridiagonalized (istart=1 and iblk are given as the block width). In step S 16 , the result of the tridiagonalization is returned to array A, i.e., A(nbase+1:n, nbase+1:nbase+iblk)←U(nbase+1:n, In step S 17 , subroutine “update” is called in order to the lower triangular portion of A(nbase+iblk:n, nbase+iblk:n). Thereafter, the process returns to step S 12 . 
     In step S 18 , setting is performed so that nbase=(nb−1)*iblk, istart=1, and iblk2=n−nbase. In S 19 , the target area of tridiagonalization is copied onto work area U, i.e., U(nbase+1:n, 1:nwidth)←A(nbase+1:n, nbase+1:n). In step S 20 , subroutine blktrid is called, and the portion copied onto U is tridiagonalized (istart=1 and iblk2 is given as the block width). In step S 21 , the result of the tridiagonalization is returned to array A, i.e., A(nbase+1:n, nbase+1:n)←U(nbase+1:n, 1:width). In step S 22 , the threads generated for the parallel processing are deleted, and the process exits from this subroutine. 
       FIG. 13  illustrates a flowchart of subroutine blktrid. 
     This subroutine is a recursive program. The program statement that is to be called is subroutine blktrid (A, k, n, diag, sdiag, nbase, istart, nwidth, U, V, nothrd, numthrd). In this statement, nbase represents the offset indicating the position of the block, istart represents the offset in the block in the divided block that is the target in the recursive calling, and the initial value of istart is 1, and indicates the position of the target block. nwidth represents the block width. 
     In step S 25 , it is determined whether or not nwidth&lt;10 is satisfied. When the determination result in step s 25  is Yes, subroutine btunit is called to perform tridiagonalization, and the process exits from this subroutine. When the determination result in step S 25  is No, the subroutine is called with the target update position and block width changed so as to perform the recursive calling. Setting is performed so that istart2=istart and nwidth2=nwidth/2. The starting point and the block width of the divided block are given. In step S 28 , subroutine blktrid is called recursively. In step S 29 , barrier synchronization is achieved between the threads. In step S 30 , point of origin (is2, is3) and point of termination (ie2, ie3) treated in each thread for updating are calculated. Specifically, those two points are calculated using the following expressions.
 
 i start3 =i start+ n width/2 ,n width3 =n width− n width/2, is 2= i start2 ,ie 2 =i start+ n width2−1 ,is 3 =i start3 ,ie 3 =i start3 +n width3−1, iptr=n base+ i start3,len( n−iptr+ numthrd−1)/numthrd, is=iptr+ (nothrd−1)*len+1, ie =min( n,iptr +nothrd*len)
 
     In step S 31 , setting is performed so that U(is:ie, is3:ie3)=U(is:ie, is3:ie3)-U(is:ie, is2:ie2)*W(is3:ie3, is2:ie2) t −W(is:ie, is2:ie2)*U(is3:ie3, is2:ie2) t  is satisfied. 
     In step S 32 , barrier synchronization is established between threads, and subroutine blktrid is called recursively, and thereafter, the process exits from this subroutine. 
       FIG. 14  illustrates a flowchart for subroutine btunit. 
     The calling of this subroutine is performed by using btunit (A, k, n, diag, sdiag, nbase, istart, nwidth, U, V, nothrd, numthrd). In the explanation below, array tmp is a temporary array for calculation. In step S 35 , tmp(numthrd), sigma, and alpha are assigned with the attribute “shared”. In step S 36 , it is determined whether or not nbase+istart&gt;n−2 is satisfied. When the determination result in step S 36  is Yes, the process exits from this subroutine. When the determination result in step S 36  is No, setting is performed in step S 37  so that i=istart, and it is determined whether or not i&lt;=istart-1+nwidth is satisfied in step S 38 . When the determination result in step S 38  is NO, the process exits from this subroutine. When the determination result in step S 38  is Yes, point of origin (is) and point of termination (ie) treated by each thread are calculated in step S 39 . Specifically, the following expressions are calculated.
 
 iptr 2= n base+ i ,len=( n−iptr 2+numthrd−1)/numthrd, is=iptr 2+(nothrd−1)*len+1, ie =min( n,iptr 2+nothrd*len)
 
     In step S 40 , barrier synchronization is established. In step S 41 , setting is performed so that tmp (nothrd)=U (is:ie, i) t *U (is:ie, i), and in step S 42 , barrier synchronization is established. In step S 43 , it is determined whether or not nothrd=1 is satisfied. When the determination result in step S 43  is No, the process proceeds to step S 45 . When the determination result in step S 43  is Yes, a calculation is performed to obtain the square root of the sum of the data as the results of the partial calculations in each thread, and calculation is performed for the tridiagonalization (generation of Householder vector). Specifically, the following expressions are calculated, where sum represents the result of the addition, and sqrt represents the square root.
 
sigma=sqrt(sum(tmp(1:numthrd)))diag( iptr 2)= u ( iptr 2 ,i ), s diag( iptr 2)=−sigma, U ( n base+ i+ 1, i )= U ( n base+1 ,i )+sign( u ( n base+i+1 ,i )*sigma alpha=1.0/(sigma* u ( n base+ i+ 1, i ), U ( iptr 2 ,i )=alpha
 
     In step S 45 , barrier synchronization is established. In step S 46 , setting is performed so that iptr3=iptr2+1. In step S 47 , subroutine matvec is called, receiving A, V, U, n, k, I, iptr3, is, ie, nothrd, and numthrd. In other words, V=A*U is calculated from the lower triangular matrix of A(iptr3:n, iptr3:n). In step S 48 , barrier synchronization is established. In step s 49 , the following expression is calculated.
 
 V ( is:ie,i )=alpha*( V ( is:ie,i )− V ( is:ie, 1 :i− 1)*( U ( iptr 3 :n, 1 :i− 1) t   *U ( iptr 3 :n,i ))− U ( is:ie, 1 :i− 1)*( V ( iptr 3 :n, 1 :i− 1) t   *U ( iptr 3 :n,i )))
 
     In step S 50 , barrier synchronization is established. In step S 51 , setting is performed so that tmp(nothrd)=V(is:ie, i) t *U(is:ie, i), and in step S 52 , barrier synchronization is established. In step S 53 , it is determined whether or not nothrd=1 is satisfied. When the determination result in step S 53  is No, the process proceeds to step s 55 . When the determination result in step S 53  is Yes, setting is performed in step S 54  so that beta=0.5*alpha*sum(tmp(1:numthrd)), where sum represents the result of the addition of the vectors. In step S 55 , barrier synchronization is established, and in step S 56 , setting is performed so that V(is:ie, i)=V(is:ie, i)-beta*U(is:ie, i), and thereafter barrier synchronization is established in step S 57 . In step S 59 , it is determined whether or not ptr2&lt;n−2 is satisfied. When the determination result in step S 59  is Yes, U(is:ie, i+1)=U(is:ie, i+1)−U(is:ie, istart:i)*V(i+1, istart:i) t −V(is:ie, istart:i)*V(n+1, istart:i) t  is calculated in step S 60 , and the process returns to step S 38 . When the determination result in step S 59  is No, U(is:ie, i+1:i+2)=U(is:ie, i+1:i+2)−U(is:ie, istart:i)*V(i+1:n, istart:i) t −V(is:ie, istart:i)*U(n+1:n, istart:i) t  is calculated in step S 61 , and the process exits from the subroutine. 
       FIG. 16  illustrates a flowchart for subroutine matvec. 
     A, V, U, n, k, I, is, ie, iptr3, nothrd, numthrd are used as parameters. In step S 65 , vsum(n, numthrd) is secured with the attribute “shared”. In step S 66 , initialization is performed so that vsum(iptr3:n, nothrd)=0 in each thread. The control unit calculates the interval in which lower triangular matrix A (iptr3:n, iptr3:n) is assigned to each thread so that each thread can calculate a matrix vector product. Specifically, the following expressions are calculated.
 
 nn=n−iptr 3+1 ,x left=1.0−dble(nothrd−1)/dble(numthrd), x right=1.0−dble(nothrd)/dble(numthrd), ms=n base2 +nn *(1.0 −d qsrt( x left))+1, me=n base2 +nn*( 1.0 −d sqrt( x right))
 
     In step S 68 , matvec1 is called. In each thread, vsum(ms:me, nothrd)=A(me+1:n, ms:me) t *U(me+1:n,i) vsum(ne+1:n,nothrd)=A(me+1:n,ms:me)*U(ms:me) 
     are calculated from A(me+1:n, ms:me) in the diagonal matrix portion whose lower triangular matrix has the second dimensional scope with point of origin ms and point of termination me. The parameters are A, vsum, U, I, ms, me, is, ie, nothrd, and numthrd. 
     In step S 69 , matvec2 is called. In each thread, vsum(ms:me)=vsum(ms:me)+LOW(A(ms:me, ms:me))*U(ms:me,i)+NDLOW(A(ms:me, ms:me)) t *U(ms:me,i) are calculated from the diagonal matrix portion whose lower triangular matrix has a second dimensional scope including point of origin ms and point of termination me. 
     LOW represents the lower triangular matrix portion including an diagonal element, and NDLOW represents the lower triangular matrix portion whose diagonal element is zero. The parameters are A, vsum, U, I, ms, me, is, ie, nothrd, and numthrd. 
     In step S 70 , barrier synchronization is established, and in step S 71 , expression 2 below is calculated, and thereafter the process exits from this subroutine. 
     
       
         
           
             
               
                 
                   
                     V 
                     ⁡ 
                     
                       ( 
                       
                         
                           is 
                           ⁢ 
                           
                             : 
                           
                           ⁢ 
                           ie 
                         
                         , 
                         i 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       numthrd 
                     
                     ⁢ 
                     
                       vsum 
                       ⁡ 
                       
                         ( 
                         
                           
                             is 
                             ⁢ 
                             
                               : 
                             
                             ⁢ 
                             ie 
                           
                           , 
                           nothrd 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   [ 
                   
                     Expression 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     2 
                   
                   ] 
                 
               
             
           
         
       
     
       FIG. 17  illustrates a flowchart of subroutine matvec1. 
     In step S 75 , setting is performed so that k2=ms. In step S 76 , it is determined whether or not k2&lt;=me is satisfied. When the determination result in step S 76  is No, the process exits from this subroutine. When the determination result in step S 76  is Yes, setting is performed in step S 77  so that k1=me+1, sum=0.0, and tmp1=U(k2,i). In step S 78 , it is determined whether or not k1&lt;=n is satisfied. When the determination result in step S 78  is Yes, setting is performed in step S 79  so that tmp2=U(k1, i), tmp=A(k1, k2) (calculated twice with the element of A being loaded only once), vsum(k1)=vsum(k1)+tmp*tmp1, sum=sum+tmp*tmp2, and k1=k1+1, and thereafter, the process returns to step S 78 . When the determination result in step S 78  is No, setting is performed so that vsum(k2)=sum and k2=k2+1 in step S 80 , and the process returns to step S 76 . 
       FIG. 18  illustrates a flowchart for subroutine matvec2. 
     In step S 85 , setting is performed so that k2=ms. In step S 86 , it is determined whether or not k2&lt;=me is satisfied. When the determination result in step S 86  is No, the process exits from this subroutine. When the determination result in step S 86  is Yes, setting is performed in step S 87  so that k1=k2+1, tmp1=U(k2, i), tmpz=A(k2, k2), and sum=tmp1*tmpz. In step S 88 , it is determined whether or not k1&lt;=me is satisfied. When the determination result in step S 88  is Yes, setting is performed in step S 89  so that tmp2=U(k1, i), tmp=A(k1, k2) (calculated twice with the element of A being loaded only once), vsum(k1)=vsum(k1)+tmp*tmp1, sum=sum+tmp*tmp2, and k1=k1+1, and thereafter the process returns to step S 88 . When the determination result in step S 88  is No, setting is performed in step S 90  so that vsum(k2)=vsum(k2)+sum and k2=k2+1, and the process returns to step S 86 . 
       FIG. 19  illustrates a flowchart for a subroutine “update”. 
     In step S 95 , barrier synchronization is established. Step S 96  is executed by the load control unit. In step S 96 , pairs are formed of individual threads to determine the point of origin and the point of termination for assigning the processes of updating. Specifically, the following expressions are calculated.
 
 n base2 =n base+ i blk,
 
 nn=n−n base2,
 
 x left=1.0−dble(nothrd−1)/dble(numthrd),
 
 x right=1.0−dble(nothrd)/dble(numthrd),
 
 is=n base2 +nn *(1.0 −d qsrt( x left))+1,
 
 ie=n base2 +nn *(1.0 −d sqrt( x right))
 
     In step S 97 , the following expression is calculated.
 
 A ( ie+ 1: n,is;ie 1)= A ( ie+ 1: n,is:ie )− W ( ie+ 1: n, 1:blk)* U ( is:ie, 1:blk) t   −U ( ie+ 1: n, 1:blk)* W ( is:ie, 1:blk) t  
 
     In step S 98 , the lower triangular matrix portion is updated so that (A=A−W*U t ). Subroutine ltgemmtrt is called, and the lower triangular matrix portion of the diagonal block matrix is updated. W(is:ie, 1:blk), U(is:ie, 1:blk), A(is:ie, is:ie) and len=ie−is+1 (size of the diagonal block matrix) are given as the parameters. 
     In step S 99 , the lower triangular matrix is updated so that (A=A−U*W t ). Subroutine ltgemmtrt is called, and the lower triangular matrix portion of the diagonal block matrix is updated. U(is:ie, 1:blk), W(is:ie, 1:blk) and len=ie−is+1 (size of the diagonal block matrix) are given as the parameters. In other words, W is replaced with U to call the subroutine. 
     In step S 100 , barrier synchronization is established between the threads, and the process exits from this subroutine. 
       FIG. 20  illustrates a flowchart for subroutine ltgemmtrt. 
     This subroutine updates the lower triangular matrix portion of a diagonal matrix, and is called a recursive program. W(is:ie, 1:blk), U(is:ie, 1:blk), A(is:ie, is:ie) are received as W(len,blk), U(len, blk), and A(len, len) as parameters. len represents the size of the matrix. The block is divided repeatedly, and when the size of the block has reached the prescribed minimum value, the process of dividing a block is terminated. For example, the minimum value of the block is set so that nwc=8. In step S 105 , it is determined whether or not len&lt;=nwc is satisfied. When the determination result in step S 105  is Yes, A(len, len)=A(len, len)−W(len, blk)*U(len, blk) t  is calculated in step S 106 . When the determination result in step S 105  is No, setting is performed so that len1=len/2, len2=len−len1, and nptr=len1+1 in step S 107 , and subroutine ltgemmtrt is called to update the lower triangular matrix portion of the diagonal block matrix in step S 108 . For this update, U(1:len1, 1:blk), W(1:len, 1:blk), A(1:len, 1:len1), and len1 (the size of the diagonal block matrix) are used as the parameters. 
     In step S 109 , A(nptr:len, 1:len2)=A(nptr:len, 1:len2)−W(nptr:len, blk)*U(1:len2, blk) t  is calculated. In step S 110 , it is determined whether or not len1=len2. When the determination result in step S 110  is Yes, setting is performed in step S 111  so that len3=len2 and nptr2=nptr, and the process proceeds to step S 113 . When the determination result in step S 110  is No, setting is performed in step S 112  so that len3=len1 and nptr2=nptr+1, and thereafter the process proceeds to step S 113 . In step S 113 , subroutine ltgemmtrt is called, and the lower triangular matrix portion of the diagonal block matrix is updated. U(nptr2:len, 1:blk), W(nptr2:len, 1:blk), A(nptr2:len, nptr2:len), and len3 (size of the diagonal block) are used as the parameters. After step S 113 , the process exits from this subroutine. 
     All examples and conditional language recited herein are intended for pedagogical purposes to aid the reader in understanding the invention and the concepts contributed by the inventor to furthering the art, and are to be construed as being without limitation to such specifically recited examples and conditions, nor does the organization of such examples in the specification relate to a showing of the superiority and inferiority of the invention. Although the embodiment(s) of the present invention has(have) been described in detail, it should be understood that the various changes, substitutions, and alterations could be made hereto without departing from the spirit and scope of the invention.