Information processing program, information processing device, and information processing method

Abstract

To improve processing efficiency of a parallel process using a matrix.SOLUTION: A computer determines, by using each of a plurality of processes included in a matrix process as a first process and using a process next to the first process as a second process, a synchronization method for one or a plurality of processing units that process elements of a first portion of a matrix in parallel in the first process. At this time, the computer determines the synchronization method based on the number of the one or the plurality of processing units that process the elements of the first portion of the matrix in parallel and the number of one or a plurality of processing units that process elements of a second portion of the matrix in the second process. The computer executes the first process by using the one or the plurality of processing units that process the elements of the first portion of the matrix in parallel, and executes a synchronization process on the one or the plurality of processing units that process the elements of the first portion of the matrix in parallel by using the determined synchronization method. The computer executes the second process by using the one or the plurality of processing units that process the elements of the second portion of the matrix in parallel.SELECTED DRAWING: Figure 31

Description

[Technical Field]

[0001] The present invention relates to information processing technology. [Background technology]

[0002] In scientific and technical calculations, huge sparse matrices are sometimes used. Storing a sparse matrix in memory using a dense matrix format increases memory consumption, so a sparse matrix format is used that stores only the non-zero elements of the sparse matrix in memory. Known sparse matrix formats include the COO (Coordinate) format and the CSR (Compressed Sparse Row) format (see, for example, Non-Patent Document 1). Using a sparse matrix format can reduce memory consumption and the size of read data.

[0003] Figure 1 shows an example of memory access to matrix data stored in memory. Figure 1(a) shows an example of memory access in GEMV (General Matrix-Vector Multiplication) when the coefficient matrix A is stored in memory using a dense matrix format.

[0004] In this example, the matrix vector product is calculated as follows:

[0005] Ax=y (1)

[0006] x represents a vector, and y represents the vector resulting from the multiplication of matrix A and vector x. Matrix A is a sparse matrix. Thread 1 calculates element 121 of vector y, thread 2 calculates element 122 of vector y, thread 3 calculates element 123 of vector y, and thread 4 calculates element 124 of vector y.

[0007] In this case, thread 1 accesses element 101 of matrix A and elements 111 to 115 of vector x, thread 2 accesses element 102 of matrix A and elements 111 to 115 of vector x, thread 3 accesses element 103 of matrix A and elements 111 to 115 of vector x, and thread 4 accesses element 104 of matrix A and elements 111 to 115 of vector x.

[0008] In FIG. 1(a), memory access from each thread to matrix A and vector x is consecutive access.

[0009] Figure 1(b) shows an example of memory access in sparse matrix-vector multiplication (SpMV) when matrix A is stored in memory using the CSR format. In this case, the matrix data of matrix A is stored in memory using arrays csrValA, csrRowPtrA, and csrColIndA.

[0010] Thread 1 accesses element 131 of array csrValA, element 141 of array csrRowPtrA, element 151 of array csrColIndA, and elements 111 and 112 of vector x. Thread 2 accesses element 132 of array csrValA, element 142 of array csrRowPtrA, element 152 of array csrColIndA, and elements 112 and 113 of vector x.

[0011] Thread 3 accesses element 133 of array csrValA, element 143 of array csrRowPtrA, element 153 of array csrColIndA, and elements 111, 114, and 115 of vector x. Thread 4 accesses element 134 of array csrValA, element 144 of array csrRowPtrA, element 154 of array csrColIndA, and elements 113 and 115 of vector x.

[0012] In Figure 1(b), when threads 1 to 4 access matrix A, a load imbalance occurs, and each thread accesses vector x randomly. As such, memory access to matrix data in sparse matrix format is irregular, so it is desirable to devise a formatting algorithm.

[0013] In relation to matrix calculations using sparse matrices, the preconditioned conjugate gradient method (PCG method) and the biconjugate gradient method are known (see, for example, Non-Patent Document 2 and Non-Patent Document 3). A method for solving large-scale asymmetric simultaneous linear equations is also known (see, for example, Non-Patent Document 4).

[0014] A program for DILU (Diagonal-based Incomplete LU) preprocessing is also known (see, for example, Non-Patent Document 5). A method for reducing the bandwidth of a symmetric sparse matrix is ​​also known (see, for example, Non-Patent Document 6). A method for parallelizing DIC (Diagonal-based Incomplete Cholesky) preprocessing and DILU preprocessing by coloring is also known (see, for example, Non-Patent Document 7 and Non-Patent Document 8).

[0015] Methods for synchronizing threads in a GPU (Graphics Processing Unit) environment are also known (see, for example, Non-Patent Documents 9 and 10).

[0016] A computer-implemented system for efficient sparse matrix representation and processing is also known (see, for example, US Pat. No. 6,233,629). A parallel multi-color incomplete LU decomposition preprocessor is also known (see, for example, US Pat. No. 6,233,629). [Prior art documents] [Patent documents]

[0017] [Patent Document 1] Special Announcement No. 2016-119084 [License 2] U.S. Patent and Trademark Publication No. 2015 / 0042672 [Non-licensed literature]

[0018] [Non-licensed Document 1] "cuSPARSE Library", NVIDIA, September 2021, [online], [Reiwa September 10, 2021], インターネット<URL:https: / / docs.nvidia.com / cuda / pdf / CUSPARSE_Library.pdf> [Non-licensed Document 2] EF Kaasschieter, "Preconditioned conjugate gradients for solving singular systems", Journal of Computational and Applied mathematics 24, 1988, pages 265-275. [Non-licensed Document 3] "Biconjugate Gradient Method", Wolfram MathWorld, September 19, 2021, [online], [Reiwa October 1, 2020], インターネット<URL:https: / / mathworld.wolfram.com / BiconjugateGradientMethod.html> [Non-licensed Document 4] C. Pommerell, "Solution of large unsymmetric systems of linear equations", 1992, pages 1-207. [Non-licensed Document 5] "OpenFOAM-dev / src / OpenFOAM / matrices / lduMatrix / preconditioners / DILUPreconditioner / DILUPreconditioner.C", GitHub, [online], [Retrieved October 1, 2021], Internet <URL:https: / / github.com / OpenFOAM / OpenFOAM-dev / blob / master / src / OpenFOAM / matrices / lduMatrix / preconditioners / DILUPreconditioner / DILUPreconditioner.C> [Non-patent document 6] E. Cuthill and J. McKee, "REDUCING THE BANDWIDTH OF SPARSE SYMMETRIC MATRICES", In Proceedings of the 1969 24th national conference, Association for Computing Machinery, August 1969, pages 157-172. [Non-Patent Document 7] "Challenges to Parallelize OpenFOAM Threads," Naoki Yoshifuji, November 5, 2019, [online], [Retrieved September 10, 2021], Internet<URL:https: / / proc-cpuinfo.fixstars.com / 2019 / 11 / openfoam-dic-pcg / > [Non-patent document 8] "AMGX REFERENCE MANUAL", NVIDIA, October 2017, [online], [Retrieved September 10, 2021], Internet<URL:https: / / github.com / NVIDIA / AMGX / blob / main / doc / AMGX_Reference.pdf> [Non-Patent Document 9] "CUDA 9 AND MORE", Akira Naruse, December 12, 2017, [online], [Retrieved September 10, 2021], Internet<URL:https: / / www.nvidia.com / content / apac / gtc / ja / pdf / 2017 / 1041.pdf> [Non-Patent Document 10] "AMGX / core / src / solvers / multicolor_dilu_solver.cu", GitHub, [online], [Retrieved September 10, 2021], Internet <URL:https: / / github.com / NVIDIA / AMGX / blob / main / core / src / solvers / multicolor_dilu_solver.cu> Summary of the Invention [Problem to be solved by the invention]

[0019] When DIC or DILU preprocessing is parallelized by coloring, synchronization occurs between multiple threads. Depending on the synchronization method used, the cost of synchronization may increase and processing efficiency may decrease.

[0020] Note that this problem is not limited to DIC preprocessing and DILU preprocessing using threads, but occurs in various matrix calculations using various processing units.

[0021] In one aspect, the present invention aims to improve the processing efficiency of parallel processing using matrices. [Means for solving the problem]

[0022] In one proposal, the information processing program causes a computer to perform the following processes.

[0023] The computer uses each of the multiple processes included in the matrix processing as a first process and the process next to the first process as a second process, and determines a synchronization method for one or more processing units that process elements of a first portion of the matrix in parallel in the first process. At this time, the computer determines the synchronization method based on the number of one or more processing units that process elements of the first portion of the matrix in parallel and the number of one or more processing units that process elements of a second portion of the matrix in parallel in the second process.

[0024] The computer performs a first process using one or more processing units that process elements of a first portion of the matrix in parallel, performs synchronization for the one or more processing units that process elements of the first portion of the matrix in parallel using the determined synchronization method, and performs a second process using one or more processing units that process elements of a second portion of the matrix in parallel. [Effects of the Invention]

[0025] According to one aspect, it is possible to improve the processing efficiency of parallel processing using matrices. [Brief explanation of the drawings]

[0026] [Figure 1] FIG. 10 is a diagram illustrating memory access to matrix data. [Figure 2] FIG. 1 is a diagram illustrating the algorithm of the PCG method. [Figure 3] FIG. 10 is a diagram illustrating a coefficient matrix. [Figure 4] FIG. 1 is a diagram illustrating the algorithm of the PBiCG method. [Figure 5A] FIG. 10 is a diagram illustrating object initialization processing in a DILU preprocessing program. [Figure 5B] FIG. 10 is a diagram showing preprocessing of a DILU preprocessing program. [Figure 6] FIG. 10 is a diagram illustrating a lower triangular matrix. [Figure 7] FIG. 10 is a diagram showing the repeat direction of a face. [Figure 8]FIG. 10 is a diagram illustrating four substitution processes included in DILU preprocessing. [Figure 9] FIG. 10 is a diagram showing the iteration directions in four substitution processes. [Figure 10] FIG. 10 is a diagram illustrating substitution processing included in DILU preprocessing and DIC preprocessing. [Figure 11] FIG. 10 is a diagram illustrating an upper triangular matrix. [Figure 12] FIG. 1 illustrates column coloring in an upper triangular matrix. [Figure 13] FIG. 10 is a diagram illustrating parallelization of forward substitution in DIC preprocessing. [Figure 14] FIG. 1 illustrates parallelization based on column coloring. [Figure 15] FIG. 1 illustrates row coloring in an upper triangular matrix. [Figure 16] FIG. 10 is a diagram illustrating parallelization of back substitution in DIC preprocessing. [Figure 17] FIG. 1 illustrates parallelization based on row coloring. [Figure 18] FIG. 10 is a diagram illustrating loop processing of DIC preprocessing and DILU preprocessing. [Figure 19] FIG. 10 is a diagram illustrating parallelization of forward substitution in DILU preprocessing. [Figure 20] FIG. 10 is a diagram illustrating a calculation process based on graph coloring. [Figure 21] FIG. 10 is a diagram showing the number of columns for each color in a fluid analysis simulation. [Figure 22] FIG. 1 illustrates a method for synchronization between threads. [Figure 23] FIG. 10 is a diagram illustrating a calculation process to which a synchronization method is applied. [Figure 24] FIG. 2 is a functional configuration diagram of an information processing device. [Figure 25] 10 is a flowchart of a calculation process. [Figure 26] FIG. 1 is a functional configuration diagram showing a specific example of an information processing device. [Figure 27] FIG. 3 is a diagram illustrating information stored in a storage unit. [Figure 28] FIG. 2 is a diagram illustrating information stored in a storage unit within a GPU. [Figure 29] FIG. 1 is a diagram illustrating synchronization processing in a calculation process using threads. [Figure 30] FIG. 10 is a diagram illustrating an upper triangular matrix. [Figure 31] FIG. 10 is a diagram illustrating a calculation process using an upper triangular matrix. [Figure 32] FIG. 1 is a diagram illustrating a synchronization method in a fluid analysis simulation. [Figure 33] FIG. 10 is a diagram illustrating a comparison result of parallelization methods. [Figure 34] 10 is a flowchart of an analysis process. [Figure 35] 10 is a flowchart of a column coloring information generation process. [Figure 36] FIG. 10 is a diagram illustrating a first color ID update process. [Figure 37] FIG. 10 is a diagram illustrating a second color ID update process. [Figure 38] 10 is a flowchart of a process for calculating a preprocessing residual of DIC preprocessing. [Figure 39] 10 is a flowchart of a process for calculating a preprocessing residual of DILU preprocessing. [Figure 40] 10 is a flowchart of an enqueue process. [Figure 41] 10 is a flowchart of thread processing. [Figure 42] 10 is a flowchart of a DIC forward assignment thread process. [Figure 43] FIG. 10 is a diagram illustrating sf for an upper triangular matrix. [Figure 44] 10 is a flowchart of a DILU forward assignment thread process. [Figure 45] 10 is a flowchart of a DILU backward substitution thread process. [Figure 46] FIG. 10 illustrates a timeline of the first DIC forward assignment thread processing. [Figure 47] FIG. 10 is a diagram showing a timeline of the second DIC forward assignment thread processing. [Figure 48] FIG. 1 is a hardware configuration diagram of an information processing device including a plurality of nodes. [Figure 49] FIG. 10 is a diagram showing performance in a fluid analysis simulation. [Figure 50] FIG. 10 is a diagram illustrating preprocessing and the number of smoothers. [Figure 51] FIG. 2 is a hardware configuration diagram of an information processing device. DETAILED DESCRIPTION OF THE INVENTION

[0027] Hereinafter, embodiments will be described in detail with reference to the drawings.

[0028] First, as an example of matrix calculation using a sparse matrix, a method for solving simultaneous linear equations such as the following will be described.

[0029] Ax=b (2)

[0030] A represents a coefficient matrix, b represents a constant vector, and x represents an unknown vector. A is a sparse matrix.

[0031] When matrix A is a symmetric matrix, equation (2) becomes a symmetric simultaneous linear equation, and when matrix A is an asymmetric matrix, equation (2) becomes an asymmetric simultaneous linear equation. As an iterative method for a symmetric simultaneous linear equation, for example, the PCG method is used, and as an iterative method for an asymmetric simultaneous linear equation, for example, the preconditioned biconjugate gradient method (PBiCG method) is used.

[0032] FIG. 2 shows an example of the algorithm of the PCG method described in Non-Patent Document 2. M in Equation 201 represents a preconditioning matrix, and M -1 represents the inverse matrix of M. The vector r i represents the residual. Equation 201 is the PCG method. -1 and r i This represents a preprocessing step to calculate the sparse matrix-vector product of M -1 and r i The preprocessing of computing the sparse matrix-vector product of is an example of matrix processing.

[0033] The if statement 202 indicates the termination condition of the for loop of the PCG method. i The calculation is terminated when the magnitude of becomes smaller than a suitable threshold ε. Equation 203 represents the sparse matrix-vector multiplication using matrix A. According to the algorithm in Figure 2, preprocessing improves the convergence of the solution.

[0034] Figure 3 shows an example of a coefficient matrix. The coefficient matrix in Figure 3 is a square matrix with 27 rows and 27 columns, and is used by a PCG solver to calculate pressure on a 3x3x3 cubic grid in a fluid analysis application. In this cubic grid, each grid point is adjacent to up to six grid points located above, below, left, right, front and rear. Small squares in the coefficient matrix represent non-zero elements.

[0035] Figure 4 shows an example of the algorithm of the PBiCG method. Equation 401 corresponds to equation 201 in Figure 2. M in equation 402 -T is M -1 represents the transpose of M -T is the transpose matrix M T Equation 401 and Equation 402 represent preprocessing in the PBiCG method. The preprocessing in the PBiCG method is an example of matrix processing. When matrix A is a symmetric matrix, the algorithm of the PBiCG method is the same as the algorithm of the PCG method.

[0036] Incomplete decomposition preprocessing is a method for reducing the amount of calculation required for preprocessing by simplifying the LDU decomposition of matrix M. The LDU decomposition of matrix M is described by the following equation:

[0037] M=LDU (3)

[0038] L represents a lower triangular matrix, D represents a diagonal matrix, and U represents an upper triangular matrix. Non-Patent Document 4 describes Jacobi preprocessing, SSOR (Symmetric Successive Over-Relaxation) preprocessing, DILU preprocessing, etc. as incomplete decomposition preprocessing. Depending on the type of incomplete decomposition preprocessing, there is a trade-off between the calculation accuracy and the calculation amount of the incomplete decomposition preprocessing.

[0039] The matrix A can be written as follows:

[0040] A=L A +D A +U A (4)

[0041] L A is a lower triangular matrix representing the lower triangular elements of matrix A, and D A is a diagonal matrix representing the diagonal elements of matrix A, and U A is an upper triangular matrix that represents the upper triangular elements of matrix A. However, the lower triangular matrix L A and the upper triangular matrix U A does not include the diagonal elements of matrix A. Hereinafter, this decomposition method will be referred to as the LDU format.

[0042] The LDU decomposition of the matrix M in the DILU preprocessing is described by the following equation using the matrix D in equation (3).

[0043] M=(L A +D)D -1 (D+U A ) (5)

[0044] In this case, the matrix M in Figure 4 -1 is described by the following equation:

[0045] M -1 =((L A +D)D -1 (D+U A )) -1 =(D+U A ) -1 D(L A +D) -1 (6)

[0046] FIG. 5A shows the object initialization process included in the DILU preprocessing program described in Non-Patent Document 5. The argument rD in statement 501 is the matrix D A and statements 502 and 503 represent the inverse matrix D of matrix D. -1Calculating and storing (D -1 rDPtr[i] (i=cell) in statement 503 represents the element d in the i-th row and i-th column of the matrix D. ii The reciprocal of 1 / d ii is stored.

[0047] In the FDIC (Faster version of DIC), rDu=a ij / d jj and rDl=a ij / d ii By pre-calculating a, the pre-processing is accelerated. ij represents the element in the i-th row and j-th column of matrix A.

[0048] Figure 5B shows the preprocessing steps included in the DILU preprocessing program described in Non-Patent Document 5. For simplicity, the names of sequences included in the program will sometimes be referred to without the Ptr. For example, rAPtr will be referred to as rA, and wAPtr will be referred to as wA.

[0049] The argument rA in the formula 505 of the for loop 504 is the residual b-Ax i The return value wA in the expression 509 of the for loop 508 is the preprocessing residual M -1 (b-Ax i ) by the for loop 504, the for loop 506, and the for loop 508. -1 (b-Ax i ) is calculated.

[0050] The expression 505 in the for loop 504 is D -1 This represents the process of calculating rA and storing it in wA. Expression 507 in for loop 506 and expression 509 in for loop 508 represent subtraction assignment to wA. Expression 507 is (L A +D) -1 This represents the forward substitution process of calculating rA and storing it in wA, and equation 509 is (D+U A ) -1 This represents the back-substitution process that calculates DwA and stores it in wA.

[0051] The for loop 506 and the for loop 508 simultaneously perform the iterative inverse matrix calculation of the lower triangular matrix and the upper triangular matrix and the matrix vector product calculation. When the processing of the for loop 508 is completed, wA contains the preprocessing residual M -1 (b-Ax i ) is stored.

[0052] The DILU preprocessing program in Non-Patent Document 5 also includes a transposed version of the preprocessing in Figure 5B. The preprocessing in Figure 5B and the transposed version are each called once for each step of the PBiCG method. Note that the DIC preprocessing program does not include the transposed version.

[0053] Next, the processing of the for loop 504 and the for loop 506 in FIG. 5B and wA←(L A +D) -1 For simplicity, we will write wA = w and rA = r. Matrix L A and matrix D is a square matrix with n rows and n columns, and w and r are n-dimensional vectors, then matrix D, matrix L A , the elements of vector w and vector r can be written as follows:

[0054]

number

[0055] w=(L A +D) -1 From r, the following equation is obtained:

[0056]

number

[0057] In the first line of both sides of equation (15), 1 / d 1,1 Multiplying by gives the following equation:

[0058]

number

[0059] In the first line of both sides of equation (16), 2,1 Multiply by and subtract the result from the second row, then add 1 / d to the second row 2,2 Multiplying by gives the following equation:

[0060]

number

[0061] By inductively considering both sides of equation (17), the following equation holds:

[0062]

number

[0063] v j represents the j-th element of the vector on the right-hand side of equation (18). The first term on the right-hand side of equation (18) corresponds to equation 505 in for loop 504, and the second term corresponds to equation 507 in for loop 506.

[0064] Figure 6 shows the lower triangular matrix L in the LDU format sparse matrix. A In this example, the matrix L A is a 4-row, 4-column matrix, and the row and column positions are represented by 0 to 3. A The non-zero elements of are a, b, c, and d.

[0065] In the lower part of Equation 507, the elements in the column direction are consecutive in column-major order. A Only non-zero elements of lower are stored. The column position corresponding to each element of lower is stored in l, and the row position is stored in u. losort uses the lower triangular matrix L A The index of lower corresponding to each non-zero element when the non-zero elements of are sorted in row-major order is stored. Therefore, lower, l, u, and losort are as follows:

[0066] lower=[a,b,c,d] l=[0,0,1,1] u=[2,3,2,3] losort=[0,2,1,3]

[0067] In this case, each array has four indices, 0 to 3. When the non-zero elements are sorted by column-major, they become a, b, c, and d, and when sorted by row-major, they become a, c, b, and d. The lower indices corresponding to a, b, c, and d are 0, 1, 2, and 3, respectively, so losort stores these indices in the order 0, 2, 1, and 3.

[0068] The index face in the for loop 506 represents the index of lower, and changes in ascending order from 0 to nFaces-1. nFaces is the value of the matrix L A represents the number of non-zero elements in the matrix L. If face ranges from 0 to nFaces-1, then sface=losort[face] is A represents the row-major index to access the non-zero elements of

[0069] If i=u[sface] and j=l[sface], then wA[u[sface]] on the left side of Eq. 507 is the w i The right-hand side of Equation 507, rD[u[sface]], corresponds to the second term d i,i -1 and lower[sface] corresponds to a i,j and wA[l[sface]] corresponds to v j Therefore, it can be seen that equation 507 represents the subtraction of the second term on the right-hand side of equation (18).

[0070] Next, the process of the for loop 508 in FIG. 5B and wA←(D+U A ) -1 For simplicity, we will write wA = z for the storage destination and wA = w for the read source.A and the elements of vector z can be written as follows:

[0071]

number

[0072] z=(D+U A ) -1 From Dw, the following equation is obtained:

[0073]

number

[0074] In the nth row of both sides of equation (21), 1 / d n,n Multiplying by gives the following equation:

[0075]

number

[0076] In the nth row of both sides of equation (22), n-1,n Multiply by and subtract the result from the n-1th row, then add 1 / d to the n-1th row n-1,n-1 Multiplying by gives the following equation:

[0077]

number

[0078] By inductively considering both sides of equation (23), the following equation holds:

[0079]

number

[0080] FIG. 7 shows an example of the direction of iteration of the index face of the for loop 508. The upper part of the expression 509 of the for loop 508 contains the row-major matrix U A Only the non-zero elements of are stored.

[0081] In the for loop 508, the index face varies in descending order from nFaces-1 to 0. nFaces is the upper triangular matrix U A In this case, the index face represents the number of non-zero elements in the matrix U A represents the row-major index to access the non-zero elements of

[0082] For sparse matrices in LDU format, the matrix L A and the location of non-zero elements in matrix U A It is assumed that the positions of non-zero elements in are symmetric about the diagonal elements. A The element in the i-th row and j-th column of the matrix U A is equivalent to the element in the jth row and ith column of

[0083] Therefore, in the same way as in the case of the for loop 504 and the for loop 506, if i=l[face] and j=u[face], then wA[l[face]] on the left side of the equation 509 of the for loop 508 is z on the left side of the equation (24). i The rD[l[face]] on the right side of Equation 509 corresponds to the d i,i -1 and upper[face] corresponds to a i,j corresponds to wA[u[face]] and q j Therefore, it can be seen that equation 509 represents the subtraction of the second term on the right-hand side of equation (24).

[0084] In the transposed version of the preprocessing in Figure 5B, M -1 Instead of M -T is used. M -T can be written as follows:

[0085] M-T =((L A +D)D -1 (D+U A )) -T =((D+U A ) -1 D(L A +D) -1 ) T =(L A +D) -T D T (D+U A ) -T =(D+U AT ) -1 D(L AT +D) -1 (31)

[0086] The AT in Equation (31) represents the transposed matrix A T T of matrix A. From Equation (31), in the preprocessing of the transposed version, by swapping the application and non-application of losort and swapping upper and lower between forward substitution and backward substitution, the matrix A T T is accessed.

[0087] In Equations (16), (17), (22), and (23), substitution is performed row by row. However, even if all rows with j > i (forward substitution) or j < i (backward substitution) are subtracted after the calculation of the i-th row, the calculation order is maintained. Therefore, in the preprocessing of the transposed version, access to non-zero elements in column-major order is performed in both forward substitution and backward substitution.

[0088] Figure 8 shows an example of four substitution processes included in the DILU preprocessing. Transpose indicates whether it is a process of the transposed version. The transpose checkmark indicates that it is a process of the transposed version, and no checkmark indicates that it is not a process of the transposed version. Step indicates whether it is forward substitution or backward substitution.

[0089] The face descending order indicates whether the faces change in ascending or descending order in the for loop. A check mark next to face descending order indicates that the faces change in descending order, and no check mark indicates that the faces change in ascending order.

[0090] The lower triangle is the lower triangular matrix L A or upper triangular matrix U A The check mark in the lower triangle indicates whether the lower triangular matrix L A indicates that the upper triangular matrix U is used. A Indicates that the following is to be used.

[0091] The loop body indicates the subtraction assignment formula included in the for loop. For simplicity, face is written as f and sface as sf. Similar notations may be used in the following explanations. wT is the preprocessing residual M -T (bA T x i )

[0092] Figure 9 shows examples of iteration directions in the four substitution processes shown in Figure 8. DILU forward substitution represents a non-transposed version of forward substitution, DILU backward substitution represents a non-transposed version of backward substitution, transposed DILU forward substitution represents a transposed version of forward substitution, and transposed DILU backward substitution represents a transposed version of backward substitution.

[0093] The triangle 901 is the lower triangular matrix L A , and triangle 902 represents the lower triangular part of the upper triangular matrix U A Triangle 903 represents the triangle obtained by transposing triangle 902 (transposed upper triangle), and triangle 904 represents the triangle obtained by transposing triangle 901 (transposed lower triangle).

[0094] In DILU forward substitution, the lower triangular matrix L Ais used, and the iteration direction is the direction indicated by arrow 911. In this case, the non-zero elements within triangle 901 are accessed in row-major order. The sface ascending order indicates that the face ascending order and sf = losort[f] are applied.

[0095] In the DILU backward substitution, the upper triangular matrix U A is used, and the iteration direction is the direction indicated by arrow 912. In this case, the face descending order is applied, and the non-zero elements within triangle 902 are accessed in row-major order.

[0096] In the transposed DILU forward substitution, the upper triangular matrix U A is used, and the iteration direction is the direction indicated by arrow 913. In this case, the face ascending order is applied, and the non-zero elements within triangle 903 are accessed in column-major order.

[0097] In the transposed DILU backward substitution, the lower triangular matrix L A is used, and the iteration direction is the direction indicated by arrow 914. In this case, the non-zero elements within triangle 904 are accessed in column-major order. The sface descending order indicates that the face descending order and sf = losort[f] are applied.

[0098] Next, the DIC preprocessing will be described. In the DIC preprocessing, since the matrix A is a symmetric matrix, A T = A. In this case, since lower[sf] = upper[sf] holds, in the forward substitution of the DIC preprocessing, Equation 507 in Figure 5B can be transformed as follows.

[0099] wA[u[sf]] -= rD[u[sf]] * upper[sf] * wA[l[sf]] (32)

[0100] After the calculation of the i-th row is completed, even if all the j-th rows where j > i (forward substitution) or j < i (backward substitution) are subtracted, the calculation order is maintained. Therefore, Equation (32) can be transformed as follows.

[0101] wA[u[f]] -=rD[u[f]]*upper[f]*wA[l[f]] (33)

[0102] By using equation (33), random access by losort is eliminated, and memory access to upper, l, and u is made more efficient. Furthermore, by not using lower, memory usage is reduced.

[0103] 10 shows an example of the substitution process included in the DILU preprocessing and DIC preprocessing. The forward substitution in the DIC preprocessing uses Equation (33), and the backward substitution in the DIC preprocessing uses the same calculation formula as the backward substitution in the DILU preprocessing.

[0104] In the calculation formulas of the loop body shown in Figures 8 and 10, a value is read from wA or wT on the right side and a value is written to wA or wT on the left side, so there is a data dependency within the loop. Graph coloring is used as a method for parallelizing loop processing that has such data dependency.

[0105] For example, in the forward substitution of the DIC preprocessing shown in Figure 10, the loop body contains upper. In this case, the upper triangular matrix U A By coloring the non-zero elements of using graph coloring, the loop processing can be parallelized without compromising the data dependency.

[0106] Figure 11 shows the upper triangular matrix U A The matrix U in Figure 11 is an example. A is a 4-by-4 matrix, and the row and column positions are represented by 0 to 3. A The non-zero elements of are a, b, c, and d. In this case, nFaces=4, and upper, u, and l are as follows:

[0107] upper=[a,b,c,d] u=[1,3,2,3] l=[0,0,1,1]

[0108] Graph coloring of a matrix can be performed using either column coloring or row coloring, where each vertex of the colored graph is associated with a column, or row coloring, where each vertex of the colored graph is associated with a row.

[0109] For example, in graph coloring based on "A Nodal Numbering Scheme" described in Non-Patent Document 6, vertices that do not have edge endpoints and vertices adjacent only to vertices that have already been reached are colored one color at a time. Column coloring is used in the parallelization of forward substitution in DIC preprocessing, and row coloring is used in the parallelization of backward substitution in DIC preprocessing.

[0110] Figure 12 shows the upper triangular matrix U A First, we consider the upper triangular matrix U A is regarded as an adjacency matrix, and the upper triangular matrix U A is converted into a directed graph 1201. Vertices 1211 to 1214 of directed graph 1201 correspond to column positions 0 to 3, respectively.

[0111] A non-zero element in row i, column j is associated with an edge extending from the vertex indicating column position i to the vertex indicating column position j. For example, non-zero element a in row 0, column 1 is associated with an edge extending from vertex 1211 indicating column position 0 to vertex 1212 indicating column position 1. Similarly, b is associated with an edge extending from vertex 1211 to vertex 1214, c is associated with an edge extending from vertex 1212 to vertex 1213, and d is associated with an edge extending from vertex 1212 to vertex 1214.

[0112] The initial value of the color ID is -1, and each time coloring is performed, the color ID is incremented by 1. Only columns that do not contain non-zero elements correspond to the color ID "-1", and since such columns are not used in preprocessing, the color IDs of columns used in preprocessing start from 0.

[0113] Because vertex 1211 does not have an end point of an edge (end point of an arrow), a color ID of "-1" is assigned to vertex 1211. Next, a color ID of "0" is assigned to vertex 1212, which can be reached from vertex 1211. Next, a color ID of "1" is assigned to vertex 1213, which can be reached from vertex 1212, and to vertex 1214, which can be reached from vertex 1211 and vertex 1212.

[0114] Next, the color ID assigned to each vertex is assigned to the column position indicated by the vertex. For example, color ID "0" of vertex 1212 is assigned to column position 1 indicated by vertex 1212. Similarly, color ID "1" of vertex 1213 is assigned to column position 2, and color ID "1" of vertex 1214 is assigned to column position 3. As a result, a is colored in the color indicated by color ID "0," and b, c, and d are colored in the color indicated by color ID "1."

[0115] Column coloring allows calculations using non-zero elements in one or more columns colored with the same color to be performed using only the results of calculations that have already been completed, so there are no data dependencies between calculations using non-zero elements of the same color, and these calculations can be performed in parallel.

[0116] 13 shows an example of parallelization of forward substitution in DIC preprocessing. Loop processing 1301 represents forward substitution in DIC preprocessing, and loop processing 1302 represents processing in which loop processing 1301 is parallelized by using column coloring.

[0117] In the loop process 1302, nColors_col represents the number of colors used in the column coloring (number of column colors), and color represents a color ID. color_cols[color] represents the column position to which the color indicated by color is assigned. nColors_col and color_cols[color] are calculated in advance.

[0118] A for loop 1311 included in the loop processing 1302 is parallelized. On the other hand, the innermost for loop 1312 is processed by sequential processing because parallelization using atomic add is expected to result in a decrease in performance.

[0119] Figure 14 shows an example of parallelization based on the column coloring of Figure 12. In the column coloring of Figure 12, the color IDs of the columns used in preprocessing are 0 and 1, so nColors_col=2. In this case, color_cols is as follows:

[0120] color_cols[0]=[1] color_cols[1]=[2,3]

[0121] When loop process 1301 in Figure 13 is expanded, it becomes process 1401. Process 1401 includes equations 1411 to 1414. Of these, equation 1411 is a calculation formula using non-zero elements assigned color ID "0", and equations 1412 to 1414 are calculation formulas using non-zero elements assigned color ID "1". Therefore, the calculations of equations 1412 to 1414 can be parallelized.

[0122] For example, when process 1401 is parallelized using thread 0 and thread 1, process 1402 can be assigned to thread 0, and process 1403 can be assigned to thread 1. Process 1402 includes expressions 1411 and 1413, and process 1403 includes expressions 1412 and 1414.

[0123] In process 1402, thread 0 executes the calculation of formula 1413 after completing the calculation of formula 1411. In process 1403, thread 1 executes the calculation of formulas 1412 and 1414 after completing the calculation of formula 1411 by thread 0.

[0124] In this case, the calculations of formulas 1412 and 1414 are performed in parallel with the calculation of formula 1413. Therefore, the calculation formula is reduced from four lines to three lines, and the speed of process 1401 is increased. In this way, by assigning calculations using non-zero elements in multiple columns colored with the same color to different threads for each column, forward substitution in DIC preprocessing can be easily parallelized.

[0125] Figure 15 shows the upper triangular matrix U A Here is an example of row coloring in the upper triangular matrix U A is regarded as an adjacency matrix, and the upper triangular matrix U A is converted into a directed graph 1501. Vertices 1511 to 1514 of directed graph 1501 correspond to row positions 0 to 3, respectively.

[0126] In row coloring, unlike column coloring, a non-zero element in row i, column j is associated with an edge from the vertex indicating row position j to the vertex indicating row position i. For example, non-zero element a in row 0, column 2 is associated with an edge from vertex 1513 indicating row position 2 to vertex 1511 indicating row position 0. Similarly, b is associated with an edge from vertex 1514 to vertex 1511, c is associated with an edge from vertex 1513 to vertex 1512, and d is associated with an edge from vertex 1514 to vertex 1513.

[0127] The initial value of the color ID is -1, and each time coloring is performed, the color ID is incremented by 1. Only rows that do not contain non-zero elements correspond to the color ID "-1", and since such rows are not used in preprocessing, the color IDs of rows used in preprocessing start from 0.

[0128] Because vertex 1514 does not have an end point of an edge, a color ID of "-1" is assigned to vertex 1514. Next, a color ID of "0" is assigned to vertex 1513, which can be reached only from vertex 1514. Next, a color ID of "1" is assigned to vertex 1512, which can be reached from vertex 1513, and to vertex 1511, which can be reached from vertex 1513 and vertex 1514.

[0129] Next, the color ID assigned to each vertex is assigned to the row position indicated by the vertex. For example, color ID "0" of vertex 1513 is assigned to row position 2 indicated by vertex 1513. Similarly, color ID "1" of vertex 1512 is assigned to row position 1, and color ID "1" of vertex 1511 is assigned to row position 0. As a result, d is colored in the color indicated by color ID "0," and a, b, and c are colored in the color indicated by color ID "1."

[0130] Row coloring allows computations using non-zero elements in one or more rows colored with the same color to be performed using only the results of computations that have already been completed, so there are no data dependencies between computations using non-zero elements of the same color, and these computations can be performed in parallel.

[0131] 16 shows an example of parallelization of back substitution in DIC preprocessing. Loop processing 1601 represents back substitution in DIC preprocessing, and loop processing 1602 represents processing in which loop processing 1601 is parallelized by row coloring.

[0132] In the loop process 1602, nColors_row represents the number of colors used in row coloring (number of row colors), and color represents a color ID. color_rows[color] represents the row position to which the color indicated by color is assigned. nColors_row and color_rows[color] are calculated in advance.

[0133] The for loop 1611 included in the loop processing 1602 is parallelized. In the case of row coloring, the start and end points of the edges of the graph are always in descending order of row position, and there is no data dependency in the innermost for loop 1612. Therefore, in the for loop 1611 and the for loop 1612, calculations are executed in ascending order.

[0134] Figure 17 shows an example of parallelization based on the row coloring of Figure 15. In the row coloring of Figure 15, the matrix U AThe non-zero elements of are a, b, c, and d, so nFaces=4. Also, the color IDs of the rows used in preprocessing are 0 and 1, so nColors_row=2. In this case, upper, u, l, and color_rows are as follows:

[0135] upper=[a,b,c,d] u=[2,3,2,3] l=[0,0,1,2] color_rows[0]=[2] color_rows[1]=[1,0]

[0136] When loop process 1601 in Figure 16 is expanded, it becomes process 1701. Process 1701 includes equations 1711 to 1714. Of these, equation 1711 is a calculation formula using non-zero elements assigned color ID "0", and equations 1712 to 1714 are calculation formulas using non-zero elements assigned color ID "1". Therefore, the calculations of equations 1712 to 1714 can be parallelized.

[0137] For example, when process 1701 is parallelized using thread 0 and thread 1, process 1702 can be assigned to thread 0, and process 1703 can be assigned to thread 1. Process 1702 includes expressions 1711 and 1712, and process 1703 includes expressions 1713 and 1714.

[0138] In process 1702, thread 0 executes the calculation of formula 1712 after completing the calculation of formula 1711. In process 1703, thread 1 executes the calculation of formulas 1713 and 1714 after completing the calculation of formula 1711 by thread 0.

[0139] In this case, the calculations of formulas 1713 and 1714 are executed in parallel with the calculation of formula 1712. Therefore, the calculation formula is reduced from four lines to three lines, and the speed of process 1701 is increased. In this way, by assigning calculations using non-zero elements in multiple lines colored with the same color to different threads for each line, it is possible to easily parallelize the backsubstitution in the DIC preprocessing.

[0140] Fig. 18 shows an example of loop processing of DIC preprocessing and DILU preprocessing. The forward substitution and backward substitution of the DILU preprocessing in Fig. 18 are the same as those in Fig. 8. The forward substitution and backward substitution of the DIC preprocessing in Fig. 18 are the same as those in Fig. 10.

[0141] DIC pretreatment -1 In the calculation, the following formula is used:

[0142] rD[u[f]] -=upper[f]*upper[f] / rD[l[f]] (34)

[0143] DILU pretreatment -1 In the calculation, the following formula is used:

[0144] rD[u[f]] -=upper[f]*lower[f] / rD[l[f]] (35)

[0145] DIC pretreatment and DILU pretreatment -1 The array to be calculated in the calculation is rD. The array to be calculated in the forward substitution and backward substitution of DIC preprocessing and DILU preprocessing is wA. The array to be calculated in the forward substitution of the transposed version and backward substitution of the transposed version of DILU preprocessing is wT.

[0146] All the loops in Figure 18 have similar data dependencies on the arrays they are calculating, so parallelization based on graph coloring can be applied to these loops.

[0147] Figure 19 shows an example of parallelizing the forward substitution of the DILU preprocessing. Loop process 1901 represents the forward substitution of the DILU preprocessing, and loop process 1902 represents the process of parallelizing loop process 1901 by column coloring. For the forward substitution of the DILU preprocessing as well, similar to the case of the forward substitution of the DIC preprocessing, it can be parallelized by column coloring.

[0148] Here, the reason why the same column coloring as that of the forward substitution of the DIC preprocessing can be used in the forward substitution of the DILU preprocessing will be explained.

[0149] As shown in Figure 18, in the forward substitution of the DILU preprocessing, different from the forward substitution of the DIC preprocessing, for the index indicating the non-zero elements of the upper triangular matrix U A sf = losort[f] is used.

[0150] First, for the index f1 and index f2 indicating the non-zero elements of the upper triangular matrix U A assume that f1 < f2 and u[f1] = l[f2] hold. In this case, there is a data dependency in the calculation order of f1 and f2.

[0151] Furthermore, assume that losort[f1] > losort[f2] holds. In this case, the application of losort reverses the calculation order of f1 and f2, and the data dependency disappears. Since losort corresponds to the lexicographical sorting of the upper u and lower l, from losort[f1] > losort[f2], any of the following relationships holds.

[0152] u[f1] > u[f2] (36) u[f1] = u[f2] and l[f1] > l[f2] (37)

[0153] From Equation (36) and u[f1] = l[f2], l[f2] > u[f2] holds, and from Equation (37) and u[f1] = l[f2], l[f1] > u[f1] holds. However, l[f2] > u[f2] means that f2 is the upper triangular matrix U AThis contradicts the fact that l[f1]>u[f1] indicates that f1 is an upper triangular matrix U A Therefore, the order of calculation of f1 and f2 does not change when losort is applied, and data dependency still exists after applying losort.

[0154] The same conclusion can be drawn in the reverse case, since the inverse of losort corresponds to a lexicographical sort of the top l and bottom u. Thus, the existence of a data dependency in a computation using two nonzero elements before applying losort is necessary and sufficient for the existence of a data dependency in a computation using two nonzero elements after applying losort.

[0155] From the above explanation, it can be seen that the same column coloring can be used in the forward substitution of DILU preprocessing after applying losort as in the forward substitution of DIC preprocessing before applying losort.

[0156] Similarly, the backward substitution of the DILU preprocessing, the forward substitution of the transposed version of the DILU preprocessing, and the backward substitution of the transposed version of the DILU preprocessing can also be parallelized by graph coloring.

[0157] Figure 20 shows an example of a calculation process based on graph coloring. The upper triangular matrix in Figure 20 is a matrix with 6 rows and 6 columns, and the row and column positions are represented by 0 to 5. The non-zero elements of this upper triangular matrix are a, b, c, d, e, f, g, and h.

[0158] In this example, column position 1 is assigned a color ID of "0," column positions 2 and 4 are assigned a color ID of "1," and column positions 3 and 5 are assigned a color ID of "2." As a result, a is colored with the color indicated by color ID "0," c, d, and f are colored with the color indicated by color ID "1," and b, e, g, and h are colored with the color indicated by color ID "2."

[0159] Calculation 2011 is a calculation using a in column position 1, calculation 2012-1 is a calculation using d in column position 2, and calculation 2012-2 is a calculation using c and f in column position 4. Calculation 2013-1 is a calculation using b, e, and g in column position 3, and calculation 2013-2 is a calculation using h in column position 5. Thus, calculations 2012-1 and 2012-2 are calculations using non-zero elements of the same color, and calculations 2013-1 and 2013-2 are also calculations using non-zero elements of the same color.

[0160] When these calculations are performed using only thread 0, thread 0 executes calculation 2011, calculation 2012-1, calculation 2012-2, calculation 2013-1, and calculation 2013-2 in chronological order.

[0161] When these calculations are performed using thread 0 and thread 1, thread 0 performs calculation 2011, calculation 2012-1, and calculation 2013-1 in chronological order. Thread 1 performs calculation 2012-2 in parallel with calculation 2012-1, and performs calculation 2013-2 in parallel with calculation 2013-1.

[0162] When multiple threads are used, synchronization operations are inserted between calculations using non-zero elements of different colors to ensure the order of calculations where data dependencies exist. In this example, after calculation 2011 is completed, synchronization operation 2021 is performed between thread 0 and thread 1, and after calculations 2012-1 and 2012-2 are completed, synchronization operation 2022 is performed between thread 0 and thread 1.

[0163] When thread 0 and thread 1 are used, the calculation time is shorter than when only thread 0 is used, but overhead occurs due to synchronization processes 2021 and 2022. Thus, in calculation processes based on graph coloring, the fewer the number of threads, the longer the calculation time, and the greater the number of threads, the greater the overhead of synchronization processes.

[0164] Figure 21 shows an example of the number of columns per color in a fluid analysis simulation. The horizontal axis represents the color ID, and the vertical axis represents the number of columns assigned the same color ID (number of columns). The maximum number of threads per block of the GPU used for calculation is 1024. In this case, calculations using up to 1024 columns can be performed in one block.

[0165] A block represents a logical grouping of threads. The maximum number of threads per block is sometimes called the block size. A block is an example of a group of threads, and the block size is an example of a group size.

[0166] Figure 21(a) shows an example of the change in the number of columns in a fluid analysis simulation that calculates pressure on a 100x100x100 cubic grid. Figure 21(b) shows an example of the change in the number of columns in a fluid analysis simulation that calculates pressure on a 290,000 unstructured grid.

[0167] 21(a) and 21(b), the peak in the number of columns is steep, and the number of columns for each color fluctuates significantly. Therefore, it is desirable to schedule an appropriate number of threads in order to set an optimal number of threads according to the number of columns for each color.

[0168] Fig. 22 shows an example of a synchronization method between threads in a GPU environment. Fig. 22(a) shows an example of stream synchronization. Kernel 2201 and kernel 2202 each include block 0 and block 1. A kernel represents a computation unit executed in time series in a GPU stream. A kernel is sometimes called a kernel function. Block 0 includes threads 0 to 2, and block 1 includes threads 0 to 2.

[0169] In stream synchronization, synchronization processing is performed for all threads between consecutively executed kernels 2201 and 2202. In this case, the preceding kernel 2201 is temporarily terminated before the subsequent kernel 2202 is started, which increases the startup cost of the kernel 2202.

[0170] 22(b) shows an example of intra-block synchronization. Kernel 2203 includes block 0 and block 1. In intra-block synchronization, synchronization process 2211 for all threads in block 0 and synchronization process 2212 for all threads in block 1 are executed separately. In this case, the number of synchronized threads is limited to the block size. Intra-block synchronization is an example of intra-group synchronization.

[0171] 22(c) shows an example of inter-block synchronization. Kernel 2204 includes block 0 and block 1. In inter-block synchronization, synchronization processing 2221 is executed for all threads included in block 0 and block 1 in kernel 2204. In this case, the cost of the synchronization processing is greater than that of intra-block synchronization. Inter-block synchronization is an example of inter-group synchronization.

[0172] Figure 23 shows an example of a calculation process applying the synchronization method shown in Figure 22. In this example, five color IDs are used.

[0173] Calculation 2311 is a calculation using the non-zero elements of the column assigned color ID "0". Calculations 2312-1 to 2312-3 are calculations using the non-zero elements of the column assigned color ID "1". Calculations 2313-1 to 2313-3 are calculations using the non-zero elements of the column assigned color ID "2". Calculation 2314 is a calculation using the non-zero elements of the column assigned color ID "3". Calculation 2315 is a calculation using the non-zero elements of the column assigned color ID "4".

[0174] 23(a) shows an example of computation processing to which stream synchronization is applied. Kernel 2301 executes computation 2311. Kernel 2302 executes computations 2312-1 to 2312-3 in parallel. Kernel 2303 executes computations 2313-1 to 2313-3 in parallel. Kernel 2304 executes computation 2314. Kernel 2305 executes computation 2315.

[0175] When stream synchronization is applied, a new kernel is started every time the color is switched, which increases the cost of starting the kernel.

[0176] 23(b) shows an example of a calculation process to which intra-block synchronization is applied. A kernel 2306 includes one block, and the block size is 2.

[0177] Of the two threads in the block, one thread executes in chronological order calculation 2311, calculation 2312-1, calculation 2312-3, calculation 2313-1, calculation 2313-3, calculation 2314, and calculation 2315. The other thread executes calculation 2312-2 in parallel with calculation 2312-1, and executes calculation 2313-2 in parallel with calculation 2313-1.

[0178] After the calculation 2311 is completed, a synchronization process 2321 is performed between the two threads, and after the calculations 2312-1 to 2312-3 are completed, a synchronization process 2322 is performed between the two threads. After the calculations 2313-1 to 2313-3 are completed, a synchronization process 2323 is performed between the two threads, and after the calculation 2314 is completed, a synchronization process 2324 is performed between the two threads.

[0179] When intra-block synchronization is applied, if the number of columns of the same color exceeds the block size, the calculation is expanded in the time direction, which increases the calculation time.

[0180] 23(c) shows an example of a calculation process in which inter-block synchronization is applied. A kernel 2307 includes two blocks, and the block size is 2.

[0181] One thread included in one of the two blocks chronologically executes calculation 2311, calculation 2312-1, calculation 2313-1, calculation 2314, and calculation 2315. The other thread executes calculation 2312-2 in parallel with calculation 2312-1, and executes calculation 2313-2 in parallel with calculation 2313-1.

[0182] One thread included in the other block executes calculation 2312-3 in parallel with calculation 2312-1 and calculation 2312-2, and executes calculation 2313-3 in parallel with calculation 2313-1 and calculation 2313-2.

[0183] After the calculation 2311 is completed, a synchronization process 2331 is performed among the three threads, and after the calculations 2312-1 to 2312-3 are completed, a synchronization process 2332 is performed among the three threads. After the calculations 2313-1 to 2313-3 are completed, a synchronization process 2333 is performed among the three threads, and after the calculation 2314 is completed, a synchronization process 2334 is performed among the three threads.

[0184] When inter-block synchronization is applied, the number of threads becomes redundant in color calculations with a small number of columns, and the cost of synchronization processing increases.

[0185] As described above, in a GPU environment, the application conditions and costs of synchronization processing differ depending on the synchronization method between threads, so it is desirable to use different synchronization methods so as to improve the efficiency of the entire computational processing.

[0186] 24 shows an example of the functional configuration of an information processing device (computer) according to the embodiment. An information processing device 2401 in FIG.

[0187] FIG. 25 is a flowchart showing an example of the calculation process performed by the information processing device 2401 in FIG.

[0188] The determination unit 2411 uses each of the multiple processes included in the matrix processing as a first process, and uses the process next to the first process as a second process. First, the determination unit 2411 determines a synchronization method for one or more processing units that process the elements of the first part of the matrix in parallel in the first process (step 2501). At this time, the determination unit 2411 determines the synchronization method based on the number of one or more processing units that process the elements of the first part of the matrix in parallel and the number of one or more processing units that process the elements of the second part of the matrix in parallel in the second process.

[0189] Next, the arithmetic processing unit 2412 executes a first process using one or more processing units that process the elements of the first part of the matrix in parallel (step 2502). Next, the arithmetic processing unit 2412 executes synchronization processing for the one or more processing units that process the elements of the first part of the matrix in parallel using the determined synchronization method (step 2503). Next, the arithmetic processing unit 2412 executes a second process using one or more processing units that process the elements of the second part of the matrix in parallel (step 2504).

[0190] According to the information processing device 2401 in FIG. 24, it is possible to improve the processing efficiency of parallel processing using matrices.

[0191] Fig. 26 shows a specific example of the information processing device 2401 in Fig. 24. The information processing device 2601 in Fig. 26 includes a CPU (Central Processing Unit) 2611, a storage unit 2612, a GPU 2613, and an output unit 2614. The GPU 2613 includes a calculation processing unit 2621 and a storage unit 2622. The CPU 2611 and the calculation processing unit 2621 correspond to the determination unit 2411 and the calculation processing unit 2412 in Fig. 24, respectively.

[0192] The information processing device 2601 uses an iterative method to find a solution to the simultaneous linear equations of equation (2) in various numerical calculations such as fluid analysis simulation, climate simulation, and molecular dynamics simulation in the fields of materials science and biochemistry. For example, the fluid analysis simulation may be a simulation that performs steady-state analysis of a fluid using the PISO (Pressure-Implicit with Splitting of Operators) method. In this case, the vector x in equation (2) represents a physical quantity. x may be a vector representing pressure or velocity.

[0193] As the iterative method, for example, the PCG method or the PBiCG method is used. The PCG method uses DIC preprocessing, and the PBiCG method uses DILU preprocessing. Either the PCG method or the PBiCG method can be selected depending on the physical quantity represented by x.

[0194] 27 shows an example of information stored in the storage unit 2612 of FIG. 26. The storage unit 2612 of FIG. 27 stores an upper triangular matrix U A and the lower triangular matrix L A Each CPU stores its own CPU coloring information 2711. The CPU coloring information 2711 includes CPU column coloring information and CPU row coloring information. The CPU column coloring information includes nColors, maxThreads, threads, and kernels.

[0195] nColors represents the number of colors used in column coloring. threads is an array representing the number of columns per color. maxThreads represents the maximum number of columns per color. kernels is an array representing the launch information for each kernel used in the iteration method.

[0196] The kernel launch information includes a start offset and a flag. The start offset indicates the index of threads that corresponds to the color of the non-zero element used for calculation by the thread when the kernel is launched. The flag indicates whether the number of threads in the kernel is greater than the block size. For implementation simplicity, the start offset of the last element of kernels is set to nColors, and the flag is set to N / A (not applicable).

[0197] The CPU row coloring information includes the same information as the CPU column coloring information. In the CPU row coloring information, nColors represents the number of colors used in row coloring, threads represents the number of rows per color, and maxThreads represents the maximum number of rows per color.

[0198] 28 shows an example of information stored in the storage unit 2622 in the GPU 2613 in FIG. 26. The storage unit 2622 in FIG. 28 stores sparse matrix information 2811, an upper triangular matrix U A and the lower triangular matrix L A Each GPU stores coloring information 2812 and an iteration method object 2813 .

[0199] The sparse matrix information 2811 is information about a sparse matrix in LDU format, and includes diag, upper, lower, u, l, ownerStart, losortStart, and losort.

[0200] diag is an array representing the diagonal elements of the coefficient matrix A, and upper is a row-major, upper triangular matrix U A is an array representing the non-zero elements of the matrix L, and lower is a column-major, lower triangular matrix L. A is an array representing the non-zero elements of

[0201] u is the upper triangular matrix U A is an array representing the column positions of non-zero elements in the upper triangular matrix U A In the LDU format sparse matrix, the upper triangular matrix U Aand the location of non-zero elements in the lower triangular matrix L A It is assumed that the positions of non-zero elements in are symmetric with respect to the diagonal elements. Therefore, u is a lower triangular matrix L A represents the row position of the non-zero elements of the lower triangular matrix L A represents the column position of a non-zero element in

[0202] ownerStart is the upper triangular matrix U A This is an array that represents the index of upper corresponding to the first non-zero element of each row when the non-zero elements of are sorted by row-major. Each element of ownerStart corresponds to the number of non-zero elements in the row above that position.

[0203] losortStart is the upper triangular matrix U A losortStart is an array that represents the index of lower corresponding to the first non-zero element in each column when the non-zero elements of are sorted by column-major. Each element of losortStart corresponds to the number of non-zero elements in the column to the left of that column position.

[0204] losort is an upper triangular matrix U A is an array that represents the index of upper corresponding to each non-zero element when the non-zero elements of are sorted by column-major.

[0205] By using such sparse matrix information 2811, the upper triangular matrix U A or lower triangular matrix L A It is possible to read out information about non-zero elements contained in a specific column or row of in the order of a constant or the order of the number of non-zero elements in that column or row. The information about non-zero elements that is read out is the number, positions, and values ​​of the non-zero elements.

[0206] The GPU coloring information 2812 includes GPU column coloring information and GPU row coloring information, which includes sorted, starts, and threads.

[0207] "sorted" is an array that represents the column positions sorted by color, and "starts" is an array that represents the starting offset for each color in "sorted." The starting offset represents the index in "sorted" that corresponds to the first column position assigned to each color. "threads" is an array similar to "threads" in the CPU column coloring information.

[0208] The GPU row coloring information contains the same information as the GPU column coloring information. The "sorted" in the GPU row coloring information represents the row position sorted by color, and "starts" represents the index of "sorted" corresponding to the first row position where each color is assigned. The "threads" is an array similar to the "threads" in the CPU row coloring information.

[0209] The iteration method object 2813 includes rD, rA, wA, rT, wT, and iteration method data. rD is an array representing the reciprocal of the diagonal elements. rA is the residual b-Ax i is an array representing the preprocessing residual M -1 (b-Ax i ) is an array representing the transposed residual bA T x i is the array representing the transposed preprocessing residual M -T (bA T x i ) is an array representing the solution x at step i of the iteration method. i Includes data such as:

[0210] For example, in the case of the upper triangular matrix shown in FIG. 20, the sparse matrix information 2811 is as follows:

[0211] upper=[a,b,c,d,e,f,g,h] u=[1,3,4,2,3,4,3,5] l=[0,0,0,1,1,1,2,2] ownerStart=[0,3,6,8,8,8,8] losortStart=[0,0,1,2,5,7,8] losort=[0,3,1,4,6,2,5,7]

[0212] In this case, upper, u, l, and losort have eight indices ranging from 0 to 7. diag and lower are omitted. Sorting the non-zero elements by row-major gives a, b, c, d, e, f, g, and h, and sorting by column-major gives a, d, b, e, g, c, f, and h.

[0213] The first non-zero element in row position 0 is a, and the index of a in upper is 0, so ownerStart[0] = 0. The first non-zero element in row position 1 is d, and the index of d in upper is 3, so ownerStart[1] = 3. The first non-zero element in row position 2 is g, and the index of g in upper is 6, so ownerStart[2] = 6.

[0214] Since there are no non-zero elements in row positions 3 to 5, ownerStart[3] = ownerStart[4] = ownerStart[5] = 8. Element 8 indicates that there is no corresponding index in upper. ownerStart[6] represents the number of non-zero elements in all rows, so ownerStart[6] = 8.

[0215] The upper indices corresponding to column-major a, d, b, e, g, c, f, and h are 0, 3, 1, 4, 6, 2, 5, and 7, respectively, so losort = [0, 3, 1, 4, 6, 2, 5, 7].

[0216] In the case of the upper triangular matrix shown in FIG. 20, the CPU column coloring information of the CPU coloring information 2711 and the GPU column coloring information of the GPU coloring information 2812 are as follows:

[0217] nColors=3 maxThreads=2 sorted=[1,2,4,3,5] starts=[0,1,3,5] threads=[1,2,2] kernels=[(0,false),(1,true),(3,N / A)]

[0218] In this example, the block size is 1. There are five indices for sorted, from 0 to 4. Element 1 at index 0 represents column position 1, which is assigned color ID "0". Element 2 at index 1 and element 4 at index 2 represent column positions 2 and 4, which are assigned color ID "1". Element 3 at index 3 and element 5 at index 4 represent column positions 3 and 5, which are assigned color ID "2".

[0219] The first column position with color ID "0" is 1, and the index of column position 1 is 0, so starts[0] = 0. The first column position with color ID "1" is 2, and the index of column position 2 is 1, so starts[1] = 1.

[0220] The first column position with color ID "2" is 3, and the index of column position 3 is 3, so starts[2]=3. Color ID "3" is unused, so starts[3]=5. The element 5 indicates that the corresponding index does not exist in sorted.

[0221] The CPU 2611 regards one or more columns assigned the same color as a portion of an upper triangular matrix and determines a synchronization method for a first process using the non-zero elements of the first portion and a second process using the non-zero elements of the second portion. The synchronization method is determined using the number N1 of one or more threads that process the non-zero elements of the first portion in the first process and the number N2 of one or more threads that process the non-zero elements of the second portion in the second process. A thread is an example of a processing unit.

[0222] If either N1 or N2 is equal to or smaller than the block size and the other is greater than the block size, the CPU 2611 determines that the synchronization method is stream synchronization. Therefore, stream synchronization is applied where the number of columns with the same color changes across block sizes.

[0223] If both N1 and N2 are greater than the block size, the CPU 2611 determines the synchronization method to be inter-block synchronization, and if both N1 and N2 are equal to or less than the block size, the CPU 2611 determines the synchronization method to be intra-block synchronization.

[0224] 20, in the case of calculation processing using threads 0 and 1, calculation 2011 is regarded as the first processing, calculation 2012-1 and calculation 2012-2 are regarded as the second processing, and the synchronization method of synchronization processing 2021 is determined. In this case, since N1=1=block size and N2=2>block size, the synchronization method of synchronization processing 2021 is determined to be stream synchronization.

[0225] Next, regarding calculation 2012-1 and calculation 2012-2 as the first process, and calculation 2013-1 and calculation 2013-2 as the second process, a synchronization method is determined for synchronization process 2022. In this case, since N1=2>block size and N2=2>block size, the synchronization method for synchronization process 2022 is determined to be inter-block synchronization.

[0226] Fig. 29 shows an example of synchronization processing in a computation process using threads 0 and 1 in Fig. 20. Because the synchronization method of synchronization process 2021 is determined to be stream synchronization, kernels 2901 and 2902 are launched in chronological order. Kernel 2901 includes block 0, and block 0 includes thread 0. Because the block size is 1, kernel 2902 includes block 0 and block 1, and block 0 and block 1 each include thread 0.

[0227] Thread 0 of block 0 of kernel 2901 executes calculation 2011. Thread 0 of block 0 of kernel 2902 executes calculation 2012-1 and calculation 2013-1 in chronological order. Thread 0 of block 1 executes calculation 2012-2 in parallel with calculation 2012-1, and executes calculation 2013-2 in parallel with calculation 2013-1. After calculation 2012-1 and calculation 2012-2 are completed, inter-block synchronization 2911 is executed between thread 0 of block 0 and thread 0 of block 1.

[0228] kernels[0]=(0, false) represents the kernel 2901. The starting offset 0 of the kernel 2901 represents the index of element 1 of threads, and the flag false represents that the number of threads in the kernel 2901 is equal to or less than the block size.

[0229] kernels[1]=(1, true) represents the kernel 2902. The starting offset 1 of the kernel 2902 represents the index of the first element 2 of threads, and the flag true represents that the number of threads in the kernel 2902 is greater than the block size.

[0230] kernels[2]=(3,N / A) is the last element of kernels, and the starting offset 3 represents nColors.

[0231] Figure 30 shows an example of an upper triangular matrix with 10 rows and 10 columns. The row and column positions of the upper triangular matrix in Figure 30 are represented by 0 to 9. The non-zero elements of this upper triangular matrix are a, b, c, d, e, f, g, h, i, j, k, l, and m.

[0232] In this example, column position 1 is assigned a color ID of "0," column positions 2, 7, and 8 are assigned a color ID of "1," and column positions 3, 5, and 6 are assigned a color ID of "2." Furthermore, column position 4 is assigned a color ID of "3," and column position 9 is assigned a color ID of "4."

[0233] As a result, a is colored in the color indicated by color ID "0", c, d, g, and h are colored in the color indicated by color ID "1", and b, e, f, i, j, and k are colored in the color indicated by color ID "2". Furthermore, l is colored in the color indicated by color ID "3", and m is colored in the color indicated by color ID "4".

[0234] In this case, the CPU column coloring information of the CPU coloring information 2711 and the GPU column coloring information of the GPU coloring information 2812 are as follows:

[0235] nColors=5 maxThreads=3 sorted=[1,2,7,8,3,5,6,4,9] starts=[0,1,4,7,8,9] threads=[1,3,3,1,1] kernels=[(0,false),(1,true),(3,false),(5,N / A)]

[0236] In this example, the block size is 2. There are nine indices for sorted, from 0 to 8. Element 1 at index 0 represents column position 1, which is assigned color ID "0". Element 2 at index 1, element 7 at index 2, and element 8 at index 3 represent column positions 2, 7, and 8, which are assigned color ID "1".

[0237] Element 3 at index 4, element 5 at index 5, and element 6 at index 6 represent column positions 3, 5, and 6, which are assigned color ID "2." Element 4 at index 7 represents column position 4, which is assigned color ID "3." Element 9 at index 8 represents column position 9, which is assigned color ID "4."

[0238] The first column position with color ID "0" is 1, and the index of column position 1 is 0, so starts[0]=0. The first column position with color ID "1" is 2, and the index of column position 2 is 1, so starts[1]=1. The first column position with color ID "2" is 3, and the index of column position 3 is 4, so starts[2]=4.

[0239] The first column position with color ID "3" is 4, and the index of column position 4 is 7, so starts[3]=7. The first column position with color ID "4" is 9, and the index of column position 9 is 8, so starts[4]=8. Color ID "5" is unused, so starts[5]=9. The element 9 indicates that there is no corresponding index in sorted.

[0240] Fig. 31 shows an example of calculation processing using the upper triangular matrix shown in Fig. 30. Calculation 3111 is a calculation using a in column position 1, calculation 3112-1 is a calculation using d in column position 2, calculation 3112-2 is a calculation using g in column position 7, and calculation 3112-3 is a calculation using c and h in column position 8.

[0241] Calculation 3113-1 is a calculation using b, e, and i in column position 3, calculation 3113-2 is a calculation using f and j in column position 5, and calculation 3113-3 is a calculation using k in column position 6. Calculation 3114 is a calculation using l in column position 4, and calculation 3115 is a calculation using m in column position 9.

[0242] Therefore, calculations 3112-1 to 3112-3 are calculations using non-zero elements of the same color, and calculations 3113-1 to 3113-3 are also calculations using non-zero elements of the same color.

[0243] First, the calculation 3111 is regarded as the first process, and the calculations 3112-1 to 3112-3 are regarded as the second process, and the synchronization method after the first process is completed is determined. In this case, since N1=1<block size and N2=3>block size, the synchronization method is determined to be stream synchronization.

[0244] Next, calculations 3112-1 to 3112-3 are regarded as the first process, and calculations 3113-1 to 3113-3 are regarded as the second process, and the synchronization method after the first process is completed is determined. In this case, since N1=3>block size and N2=3>block size, the synchronization method is determined to be inter-block synchronization.

[0245] Next, calculations 3113-1 to 3113-3 are regarded as the first process, and calculation 3114 is regarded as the second process, and the synchronization method after the first process is completed is determined. In this case, since N1=3>block size and N2=1<block size, the synchronization method is determined to be stream synchronization.

[0246] Next, regarding calculation 3114 as the first process and calculation 3115 as the second process, the synchronization method after the first process is completed is determined. In this case, since N1=1<block size and N2=1<block size, the synchronization method is determined to be intra-block synchronization.

[0247] Since the synchronization method after the calculation 3111 is completed and the synchronization method after the calculations 3113-1 to 3113-3 are completed are determined to be stream synchronization, the kernels 3101 to 3103 are activated in chronological order.

[0248] Kernel 3101 executes calculation 3111. Kernel 3102 executes calculations 3112-1 to 3112-3 in parallel, and then executes inter-block synchronization 3121. Next, kernel 3102 executes calculations 3113-1 to 3113-3 in parallel. Kernel 3103 executes calculation 3114, and then executes intra-block synchronization 3122. Next, kernel 3103 executes calculation 3115.

[0249] kernels[0]=(0,false) represents the kernel 3101. The starting offset 0 of the kernel 3101 represents the index of the first element 1 of threads, and the flag false represents that the number of threads in the kernel 3101 is equal to or less than the block size.

[0250] kernels[1]=(1, true) represents the kernel 3102. The starting offset 1 of the kernel 3102 represents the index of the first element 3 of threads, and the flag true represents that the number of threads in the kernel 3102 is greater than the block size.

[0251] kernels[2]=(3,false) represents the kernel 3103. The starting offset 3 of the kernel 3103 represents the index of the second element 1 of threads, and the flag false represents that the number of threads in the kernel 3103 is equal to or less than the block size.

[0252] kernels[3]=(5,N / A) is the last element of kernels, and the starting offset 5 represents nColors.

[0253] The CPU 2611 enqueues kernels corresponding to each element of the kernels into a stream, causing the GPU 2613 to execute preprocessing. The calculation processing unit 2621 of the GPU 2613 executes preprocessing by activating the enqueued kernels in chronological order. The calculation processing unit 2621 executes stream synchronization processing by terminating the preceding kernel and activating the next kernel, and executes inter-block synchronization or intra-block synchronization processing depending on the number of threads included in each kernel.

[0254] FIG. 32 shows an example of a synchronization method for the fluid analysis simulation of FIG. 21(a). One kernel is launched for each of color ID ranges R1 to R3. Therefore, the number of launched kernels is three. In ranges R1 and R3, 1024 threads corresponding to the block size are used, and in range R2, more than 1024 threads are used.

[0255] In range R1, intra-block synchronization is performed after the calculation using each color column except the last color is completed, and stream synchronization is performed after the calculation using the last color column is completed. In range R2, inter-block synchronization is performed after the calculation using each color column except the last color is completed, and stream synchronization is performed after the calculation using the last color column is completed. In range R3, intra-block synchronization is performed after the calculation using each color column except the last color is completed.

[0256] The area of ​​rectangle 3201 represents the product of the number of colors and the number of threads in range R1, the area of ​​rectangle 3202 represents the product of the number of colors and the number of threads in range R2, and the area of ​​rectangle 3203 represents the product of the number of colors and the number of threads in range R3. The area of ​​each rectangle corresponds to the computational complexity of a calculation using the sequence of colors in the corresponding range.

[0257] According to this parallelization method, stream synchronization is applied only where the number of columns with the same color changes across block sizes. This reduces the kernel startup cost. Furthermore, when the number of columns is greater than the block size, inter-block synchronization is applied with an emphasis on the number of parallel columns, thereby promoting speedup through parallelization. When the number of columns is equal to or less than the block size, intra-block synchronization, which has a low synchronization cost, is applied, thereby reducing the synchronization cost. This improves the processing efficiency of parallel processing.

[0258] Figure 33 shows an example of the comparison results of parallelization methods based on graph coloring. The labels indicate the names of the parallelization methods, the number of launched kernels indicates the number of launched kernels, and the number of blocks indicates the number of blocks included in a kernel. The synchronization method indicates the synchronization method after the calculation using each colored column is completed, and the problems indicate the problems with the parallelization methods.

[0259] THRUST is a parallelization method that uses Thrust, a template library of CUDA (Compute Unified Device Architecture), and its synchronization method is stream synchronization. STREAM is a parallelization method that uses CUDA, and its synchronization method is also stream synchronization.

[0260] The synchronization method of BLOCK is intra-block synchronization, and the synchronization method of GRID is inter-block synchronization. ADAPTIVE is a parallelization method of an embodiment, and the synchronization method of ADAPTIVE uses intra-block synchronization, inter-block synchronization, and stream synchronization in combination.

[0261] In the case of THRUST and STREAM, the number of kernels launched is equal to the number of colors, which means the overhead of synchronization is high. In the case of THRUST, the overhead of Thrust is further added. In the case of BLOCK, the number of parallel processes is limited to the block size, and in the case of GRID, the synchronization cost between blocks is high.

[0262] On the other hand, in the case of ADAPTIVE, the number of kernels is the number of times the number of columns crosses the block size + 1. This reduces the kernel startup cost compared to THRUST and STREAM. Furthermore, by adaptively selecting inter-block or intra-block synchronization depending on whether the number of columns exceeds the block size, the number of parallel processes and synchronization cost are optimized.

[0263] When a lower triangular matrix is ​​used, the synchronization method is determined in the same manner as when an upper triangular matrix is ​​used. When row coloring information is used to determine the synchronization method, CPU 2611 regards one or more rows assigned the same color as part of an upper triangular matrix or a lower triangular matrix, and determines the synchronization method in the same manner as when column coloring information is used.

[0264] Fig. 34 is a flowchart showing an example of analysis processing performed by the information processing device 2601 in Fig. 26. The CPU 2611 executes an application program to perform the analysis processing in Fig. 34. The analysis processing corresponds to a fluid analysis simulation, a climate simulation, a molecular dynamics simulation, etc.

[0265] First, the CPU 2611 acquires sparse matrix information 2811 of each sparse matrix in the LDU format used in the analysis process from the GPU 2613. Then, the CPU 2611 uses the sparse matrix information 2811 to calculate the upper triangular matrix U A and the lower triangular matrix L A For each, column coloring information is generated (step 3401). The column coloring information includes nColors, maxThreads, threads, kernels, sorted, and starts.

[0266] Next, the CPU 2611 calculates the upper triangular matrix U A and the lower triangular matrix L A For each, row coloring information is generated (step 3402). The row coloring information includes nColors, maxThreads, threads, kernels, sorted, and starts.

[0267] Next, the CPU 2611 uses the GPU 2613 to calculate the D -1 Calculation is performed (step 3403). Then, CPU 2611 repeats the time step loop process of steps 3404 to 3407. A time step represents a processing step at a predetermined time interval to which an iterative method is applied in the analysis process.

[0268] In the processing of the time step loop, the CPU 2611 instructs the GPU 2613 to calculate the first residual (step 3404). The GPU 2613 calculates the first residual. The first residual is b-Ax0 in FIG. 2 and FIG. 4 and bAx0 in FIG. 4. T Next, the CPU 2611 repeats the solver loop processing of steps 3405 and 3406.

[0269] A solver represents a numerical calculation that finds a solution to the simultaneous linear equations in equation (2), and one or more solvers are used depending on the type of physical quantity represented by x. For example, in the case of a fluid analysis simulation, a pressure solver that finds the pressure at each grid point and a velocity solver that finds the velocity at each grid point may be used. The solver loop processing is repeated for each solver.

[0270] In the solver loop processing, the CPU 2611 repeats the iteration loop processing of steps 3405 and 3406. The iteration method used is the PCG method or the PBiCG method. In the iteration loop processing, the CPU 2611 calculates a preprocessing residual using the GPU 2613 (step 3405). In the PCG method, a preprocessing residual in DIC preprocessing is calculated, and in the PBiCG method, a preprocessing residual in DILU preprocessing is calculated.

[0271] Next, the CPU 2611 instructs the GPU 2613 to perform iterative calculations following the calculation of the preprocessing residual (step 3406). The GPU 2613 performs iterative calculations following the calculation of the preprocessing residual. The iterative calculations following the calculation of the preprocessing residual include updating the solution and residual using the preprocessing residual. The iterative loop processing is repeated until the residual, etc., satisfies the termination condition.

[0272] After the solver loop processing iterations are completed, the CPU 2611 instructs the GPU 2613 to perform other processing (step 3407). The GPU 2613 performs the other processing, and the output unit 2614 outputs the processing results of the other processing. The other processing includes calculations using the solutions obtained by the solver loop processing. The time step loop processing is repeated until the final time step in the analysis processing is reached.

[0273] Fig. 35 is a flowchart showing an example of the column coloring information generation process in step 3401 of Fig. 34. In the column coloring information generation process of Fig. 35, the upper triangular matrix U A and the lower triangular matrix L A Each of these is used as a matrix to be processed.

[0274] First, the CPU 2611 initializes an array cc that indicates the color ID of each column of the matrix to be processed (step 3501), thereby setting each element of cc to the initial value −1.

[0275] Next, the CPU 2611 repeats the processing of the cell loop in step 3502. Cell is a variable indicating a non-zero element included in the matrix to be processed. In the processing of the cell loop, the CPU 2611 updates cc[u[cell]] using the following equation (step 3502).

[0276] cc[u[cell]]= max(cc[u[cell]],cc[l[cell]]+1) (41)

[0277] In equation (41), max(cc[u[cell]],cc[l[cell]]+1) represents the maximum value of cc[u[cell]] and cc[l[cell]]+1.

[0278] Equation (41) indicates that if there is an edge from vertex l[cell] to vertex u[cell] in the directed graph, the color ID next to the color ID of vertex l[cell] is assigned to vertex u[cell]. By repeating the cell loop process for each non-zero element, vertices that do not have an end point of an edge and vertices adjacent only to vertices that have already been reached are colored one by one, as shown in Figure 12, and the color of each vertex is associated with each column.

[0279] Figure 36 shows an example of the first color ID update process for the upper triangular matrix of Figure 20. u[cell] represents the column position of the non-zero element indicated by cell, and l[cell] represents the row position of the non-zero element indicated by cell. The numbers 0 to 5 to the right of l[cell] represent cc[0] to cc[5], and cc[j] (j = 0 to 5) represents the color ID at column position j.

[0280] In this case, cell varies from 0 to 7, and cc[u[cell]] is updated in ascending order of cell. For example, when cell=0, cc[u[0]]=cc[1] is calculated by the following formula.

[0281] cc[1]=max(cc[1],cc[l[0]]+1) =max(-1,cc[0]+1) =max(-1,-1+1) =max(-1,0) =0 (42)

[0282] Therefore, cc[1] is updated from -1 to 0 when cell=0. Similarly, cc[2] is updated from -1 to 1 when cell=3. cc[3] is updated from -1 to 0 when cell=1, from 0 to 1 when cell=4, and from 1 to 2 when cell=6. cc[4] is updated from -1 to 0 when cell=2, and from 0 to 1 when cell=5. cc[5] is updated from -1 to 0 when cell=7.

[0283] Finally, column position 0 is assigned a color ID of "-1", column position 1 is assigned a color ID of "0", column positions 2 and 4 are assigned a color ID of "1", and column positions 3 and 5 are assigned a color ID of "2". Thus, cc=[-1,0,1,2,1,2].

[0284] Figure 37 shows an example of the second color ID update process for the upper triangular matrix in Figure 30. The 0 to 9 on the right side of l[cell] represent cc[0] to cc[9], and cc[j] (j = 0 to 9) represents the color ID at column position j. In this case, cell ranges from 0 to 12, and cc[u[cell]] is updated in ascending order of cell.

[0285] cc[1] is updated from -1 to 0 when cell=0. cc[2] is updated from -1 to 1 when cell=3. cc[3] is updated from -1 to 0 when cell=1, from 0 to 1 when cell=4, and from 1 to 2 when cell=8. cc[4] is updated from -1 to 3 when cell=11.

[0286] cc[5] is updated from -1 to 1 when cell=5, and from 1 to 2 when cell=9. cc[6] is updated from -1 to 2 when cell=10. cc[7] is updated from -1 to 1 when cell=6. cc[8] is updated from -1 to 0 when cell=2, and from 0 to 1 when cell=7. cc[9] is updated from -1 to 4 when cell=12.

[0287] Finally, column position 0 is assigned a color ID of "-1", column position 1 is assigned a color ID of "0", column positions 2, 7, and 8 are assigned a color ID of "1", column positions 3, 5, and 6 are assigned a color ID of "2", column position 4 is assigned a color ID of "3", and column position 9 is assigned a color ID of "4". Thus, cc=[-1,0,1,2,3,2,2,1,1,4].

[0288] After the cell loop is completed, the CPU 2611 calculates nColors by adding 1 to the maximum value of the elements of cc (step 3503). Next, the CPU 2611 generates sorted by sorting the column positions having color IDs equal to or greater than 0 in ascending order of color ID, and generates starts from sorted (step 3504).

[0289] Next, the CPU 2611 generates threads by obtaining the number of columns for each color ID (step 3505), and obtains the maximum value of the elements of threads as maxThreads (step 3506).

[0290] Next, CPU 2611 regards columns assigned the same color ID as a part of the matrix to be processed and determines a synchronization method for the synchronization process between a first process using the non-zero elements of the first part and a second process using the non-zero elements of the second part. In determining the synchronization method, the number N1 of threads that processes the non-zero elements of the first part in the first process and the number N2 of threads that processes the non-zero elements of the second part in the second process are used.

[0291] If either N1 or N2 is equal to or less than the block size and the other is greater than the block size, the CPU 2611 determines the synchronization method to be stream synchronization. If both N1 and N2 are greater than the block size, the CPU 2611 determines the synchronization method to be inter-block synchronization, and if both N1 and N2 are equal to or less than the block size, the CPU 2611 determines the synchronization method to be intra-block synchronization.

[0292] Next, the CPU 2611 divides the elements of sorted at the point determined for stream synchronization, and assigns different kernels to one or more elements before the division point and one or more elements after the division point, respectively.The CPU 2611 then generates kernels whose elements are the startup information of each kernel (step 3507).

[0293] The CPU 2611 stores nColors, maxThreads, threads, and kernels in the storage unit 2612 as CPU column coloring information of the CPU coloring information 2711. Then, the CPU 2611 transfers sorted, starts, and threads to the GPU 2613.

[0294] The arithmetic processing unit 2621 of the GPU 2613 stores the sorted, starts, and threads received from the CPU 2611 in the storage unit 2622 as GPU column coloring information of the GPU coloring information 2812.

[0295] The row coloring information generation process in step 3402 of Figure 34 is similar to the column coloring information generation process of Figure 35. In the row coloring information generation process, an array cr representing the color ID of each row is used instead of the array cc representing the color ID of each column, and the following equation is used instead of equation (41):

[0296] cr[l[cell]]= max(cr[l[cell]],cr[u[cell]]+1) (43)

[0297] FIG. 38 is a flowchart showing an example of processing for calculating the preprocessing residual of DIC preprocessing in step 3405 of FIG. 34 when the PCG method is selected as the iterative method.

[0298] First, the CPU 2611 instructs the GPU 2613 to multiply rD[i] and rA[i] (step 3801). The GPU 2613 calculates the product of rD[i] and rA[i] and stores the result in wA[i].

[0299] Next, the CPU 2611 performs wA forward substitution using the GPU 2613 (step 3802), and performs wA backward substitution using the GPU 2613 (step 3803). The wA forward substitution and wA backward substitution correspond to the forward substitution and backward substitution of the DIC preprocessing, respectively. In the wA forward substitution and wA backward substitution, the GPU 2613 executes a series of kernels included in the stream to update wA.

[0300] Next, the CPU 2611 waits until the execution of all kernels included in the stream is completed (step 3804).

[0301] In a typical GPU application, the update process for wA is often enqueued to the same stream in the next iteration loop, ensuring the order of calculations. In this case, the wait process in step 3804 can be omitted.

[0302] FIG. 39 is a flowchart showing an example of processing for calculating the preprocessing residual of DILU preprocessing in step 3405 of FIG. 34 when the PBiCG method is selected as the iterative method.

[0303] First, the CPU 2611 instructs the GPU 2613 to multiply rD[i] and rA[i] (step 3901). The GPU 2613 calculates the product of rD[i] and rA[i] and stores the result in wA[i].

[0304] Next, the CPU 2611 performs wA forward substitution using the GPU 2613 (step 3902), and performs wA backward substitution using the GPU 2613 (step 3903). The wA forward substitution and wA backward substitution correspond to the forward substitution and backward substitution of the DILU preprocessing, respectively. In the wA forward substitution and wA backward substitution, the GPU 2613 executes a series of kernels included in the stream to update wA.

[0305] Next, the CPU 2611 instructs the GPU 2613 to multiply rD[i] and rT[i] (step 3904). The GPU 2613 calculates the product of rD[i] and rT[i] and stores the result in wT[i].

[0306] Next, the CPU 2611 performs wT forward substitution using the GPU 2613 (step 3905), and performs wT backward substitution using the GPU 2613 (step 3906). The wT forward substitution and wT backward substitution correspond to the transposed forward substitution and transposed backward substitution of the DILU preprocessing, respectively. In the wT forward substitution and wT backward substitution, the GPU 2613 executes a series of kernels included in the stream to update wT.

[0307] Next, the CPU 2611 waits until the execution of all kernels included in the stream is completed (step 3907).

[0308] In step 3403 in Fig. 34, step 3802 and step 3803 in Fig. 38, and step 3902, step 3903, step 3905, and step 3906 in Fig. 39, kernels are enqueued into streams using CPU coloring information 2711. Hereinafter, the processing performed in these steps may be referred to as enqueue processing.

[0309] 40 is a flowchart showing an example of enqueue processing. First, the CPU 2611 checks whether the faces change in descending order based on the preprocessing type, kernel type, and transposition type (step 4001).

[0310] The preprocessing type indicates whether it is DIC preprocessing or DILU preprocessing, and the kernel type is D -1 The transpose type indicates whether the forward and backward substitutions in the DILU preprocessing are transposed versions.

[0311] Whether or not the faces change in descending order can be determined from the information in Fig. 18. In this case, the faces change in descending order in the backward substitution of DIC preprocessing, the backward substitution of DILU preprocessing, and the backward substitution of the transposed version of DILU preprocessing, and the faces change in ascending order in other processes.

[0312] Therefore, in the enqueue processing of steps 3803, 3903, and 3906, the faces change in descending order, and in the enqueue processing of steps 3403, 3802, 3902, and 3905, the faces change in ascending order.

[0313] If the faces change in descending order (step 4001, YES), the CPU 2611 determines to use the kernels of the row coloring information (step 4002). On the other hand, if the faces change in ascending order (step 4001, NO), the CPU 2611 determines to use the kernels of the column coloring information (step 4003).

[0314] Next, the CPU 2611 repeats the processing of the kernels loop in steps 4004 to 4008. The processing of the kernels loop uses index k of kernels. kernels[k].first represents the start offset included in the element of index k of kernels. kernels[k].second represents the flag included in the element of index k of kernels.

[0315] Since the last element of kernels is not enqueued to the stream, if the number of elements of kernels is K, the processing of the kernels loop is repeated for k=0 to K-2.

[0316] In processing the kernels loop, the CPU 2611 determines the start color ID and end color ID using the following formula (step 4004).

[0317] Starting color ID = kernels[k].first (44) End color ID=kernels[k+1].first(45)

[0318] Next, the CPU 2611 checks kernels[k].second (step 4005). If kernels[k].second=true (step 4005, YES), the CPU 2611 sets the number of active threads to maxThreads (step 4006). On the other hand, if kernels[k].second=false (step 4005, NO), the CPU 2611 sets the number of active threads to the block size (step 4007).

[0319] Next, the CPU 2611 enqueues the kernel corresponding to the element with index k of kernels to the stream of the GPU 2613 (step 4008). The arguments of the kernel include a start color ID, an end color ID, the number of activated threads, a preprocessing type, a kernel type, and a transposition type.

[0320] The arithmetic processing unit 2621 of the GPU 2613 activates the kernels enqueued in the stream in chronological order, and executes the processing indicated by the pre-processing type, kernel type, and transposition type using the number of activated threads. If the number of activated threads is greater than the block size, multiple blocks are used, and if the number of activated threads is equal to the block size, a single block is used.

[0321] 41 is a flowchart showing an example of thread processing performed using the nth thread (n=0 to the number of activated threads−1) among the number of activated threads. First, the arithmetic processing unit 2621 of the GPU 2613 checks whether the face changes in descending order based on the preprocessing type, kernel type, and transposition type (step 4101).

[0322] If face changes in descending order (step 4101, YES), the calculation processing unit 2621 decides to use threads, sorted, and starts from the row coloring information (step 4102). On the other hand, if face changes in ascending order (step 4101, NO), the calculation processing unit 2621 decides to use threads, sorted, and starts from the column coloring information (step 4103).

[0323] Next, the calculation processing unit 2621 repeats the color loop processing of steps 4104 to 4114. In the color loop processing, a variable c indicating a color ID is used. The color loop processing is repeated for c=start color ID to end color ID−1.

[0324] In the processing of the color loop, the calculation processing unit 2621 compares n with threads[c] (step 4104). If n is smaller than threads[c] (step 4104, YES), the calculation processing unit 2621 calculates cell using the following formula (step 4105).

[0325] cell=sorted[starts[c]+n] (46)

[0326] Next, the calculation processing unit 2621 checks the kernel type and the matrix to be used (step 4106). The kernel type and the matrix to be used can be determined from the information in FIG.

[0327] The kernel type is forward or backward substitution, and the matrix used is an upper triangular matrix U A If so (step 4106, YES), the processing unit 2621 determines the start face and the end face according to the following equations (step 4107).

[0328] Start face=losortStart[cell] (47) End face=losortStart[cell+1] (48)

[0329] On the other hand, if the kernel type is D -1 The calculation is performed or the matrix used is a lower triangular matrix L A If so (step 4106, NO), the processing unit 2621 determines the start face and the end face according to the following equations (step 4108).

[0330] Start face=ownerStart[cell] (49) End face=ownerStart[cell+1] (50)

[0331] Therefore, equations (49) and (50) are used in the enqueue processing of steps 3403, 3902, and 3906. Furthermore, equations (47) and (48) are used in the enqueue processing of steps 3802, 3803, 3903, and 3905.

[0332] Next, the calculation processing unit 2621 repeats the processing of the face loop in steps 4109 to 4111. In the processing of the face loop, a variable f indicating a face is used. The processing of the face loop is repeated for f=start face to end face-1. The face loop can be expanded to increase the speed.

[0333] In the processing of the face loop, the arithmetic processing unit 2621 checks the kernel type and the matrix used (step 4109). If the kernel type is forward substitution or backward substitution and the matrix used is an upper triangular matrix U A If so (step 4109, YES), the calculation processing unit 2621 calculates sf using the following equation (step 4110).

[0334] sf=losort[f] (51)

[0335] Next, the arithmetic processing unit 2621 performs the processing of the loop body shown in FIG. 18 (step 4111).

[0336] On the other hand, if the kernel type is D -1 The calculation is performed or the matrix used is a lower triangular matrix L A If so (step 4109, NO), the processing unit 2621 skips the processing of step 4110 and performs the processing of the loop body (step 4111).

[0337] Therefore, in the enqueue processes of steps 3403, 3902, and 3906, the process of step 4110 is skipped, and in the enqueue processes of steps 3802, 3803, 3903, and 3905, the process of step 4110 is performed.

[0338] In the kernel processing in Figure 41, the default is to parallelize by columns, so the upper triangular matrix U A In the case of forward or backward substitution using (51), losort in equation (51) is applied. Then, losortStart in equations (47) and (48) is used instead of ownerStart in equations (49) and (50).

[0339] On the other hand, D -1 Calculation, lower triangular matrix L A Forward substitution using, or lower triangular matrix L A In the case of back substitution using (51), the losort in equation (51) is not applied, and the ownerStart in equations (49) and (50) is used.

[0340] Therefore, if the kernel type is forward assignment or backward assignment, in step 4111, the arithmetic processing unit 2621 calculates wA or wT using the calculation formula in which f and sf included in the loop body before parallelization shown in FIG. 18 are swapped.

[0341] On the other hand, if the kernel type is D -1 If it is a calculation, the arithmetic processing unit 2621 calculates rD ​​using the loop body before parallelization shown in FIG. 18 in step 4111. In step 3403, the enqueue process in FIG. 40 is executed for each solver. In the case of a solver using the PCG method, the DIC preprocessing D-1 The calculation loop body is used, and for solvers using the PBiCG method, the DILU preconditioner D -1 A loop body of computation is used.

[0342] Next, the processing unit 2621 checks the number of blocks included in the kernel being executed (step 4112). If the kernel includes multiple blocks (step 4112, YES), the processing unit 2621 applies inter-block synchronization to execute synchronization processing (step 4113). On the other hand, if the kernel includes only a single block (step 4112, NO), the processing unit 2621 applies intra-block synchronization to execute synchronization processing (step 4114).

[0343] However, if c=end color ID-1, the processing of steps 4112 to 4114 is omitted. If n is equal to or greater than threads[c] (step 4104, NO), the calculation processing unit 2621 performs the processing of step 4112 and subsequent steps.

[0344] FIG. 42 is a flowchart showing an example of DIC forward assignment thread processing when the preprocessing type is DIC preprocessing and the kernel type is forward assignment in the thread processing of FIG.

[0345] In this case, since face changes in ascending order, the calculation processing unit 2621 decides to use threads, sorted, and starts of the column coloring information (step 4201).

[0346] Next, the calculation processing unit 2621 repeats the color loop processing of steps 4202 to 4209. The color loop processing is repeated for c=start color ID to end color ID−1.

[0347] The processing in steps 4202 and 4203 is the same as the processing in steps 4104 and 4105 in Fig. 41. In this case, the kernel type is forward substitution, and the matrix used is an upper triangular matrix U ATherefore, the calculation processing unit 2621 determines the start face and the end face according to equations (47) and (48) (step 4204).

[0348] Next, the calculation processing unit 2621 repeats the face loop processing of steps 4205 and 4206. The face loop processing is repeated for f=start face to end face-1 in ascending order of face.

[0349] The processing of step 4205 is the same as the processing of step 4110 in Fig. 41. In this case, since the preprocessing type is DIC preprocessing and the kernel type is forward substitution, the arithmetic processing unit 2621 calculates wA by the following equation (step 4206).

[0350] wA[u[sf]] -=rD[u[sf]]*upper[sf]*wA[l[sf]] (52)

[0351] Equation (52) represents a calculation formula in which f included in the body of the forward substitution loop in the DIC preprocessing shown in Fig. 18 is replaced with sf. The processing in steps 4207 to 4209 is the same as the processing in steps 4112 to 4114 in Fig. 41.

[0352] Fig. 43 shows an example of sf for the upper triangular matrix in Fig. 20. The sparse matrix information 2811 of the upper triangular matrix in Fig. 20 is as follows:

[0353] upper=[a,b,c,d,e,f,g,h] u=[1,3,4,2,3,4,3,5] l=[0,0,0,1,1,1,2,2] losortStart=[0,0,1,2,5,7,8] losort=[0,3,1,4,6,2,5,7]

[0354] Also, sorted and starts are as follows:

[0355] sorted=[1,2,4,3,5] starts=[0,1,3,5]

[0356] When n=0 and c=2, according to formula (46), cell=3. Therefore, according to formula (47) and formula (48), the start face=2 and the end face=5.

[0357] In this case, the face loop processing is repeated for f = 2 to 4. From equation (51), when f = 2, sf = 1, when f = 3, sf = 4, and when f = 4, sf = 6. Since upper[1] = b, upper[4] = e, and upper[6] = g, the processing of step 4206 in Figure 42 is performed for b, e, and g, which are non-zero elements at column position 3.

[0358] FIG. 44 is a flowchart showing an example of DILU forward assignment thread processing when the preprocessing type is DILU preprocessing and the kernel type is forward assignment in the thread processing of FIG.

[0359] In this case, since face changes in ascending order, the calculation processing unit 2621 decides to use threads, sorted, and starts of the column coloring information (step 4401).

[0360] Next, the calculation processing unit 2621 repeats the color loop processing of steps 4402 to 4408. The color loop processing is repeated for c=start color ID to end color ID−1.

[0361] The processing in steps 4402 and 4403 is the same as the processing in steps 4104 and 4105 in Fig. 41. In this case, the kernel type is forward substitution, and the matrix used is a lower triangular matrix L A Therefore, the calculation processing unit 2621 determines the start face and the end face according to equations (49) and (50) (step 4404).

[0362] Next, the calculation processing unit 2621 repeats the processing of the face loop in step 4405. The processing of the face loop is repeated for f=start face to end face-1 in ascending order of face.

[0363] In this case, since the preprocessing type is DILU preprocessing and the kernel type is forward substitution, the arithmetic processing unit 2621 calculates wA by the following equation (step 4405).

[0364] wA[u[f]] -=rD[u[f]]*lower[f]*wA[l[f]] (53)

[0365] Equation (53) represents a calculation formula in which sf included in the body of the forward substitution loop in the DILU preprocessing shown in Fig. 18 is replaced with f. The processing in steps 4406 to 4408 is the same as the processing in steps 4112 to 4114 in Fig. 41.

[0366] FIG. 45 is a flowchart showing an example of DILU backward substitution thread processing when the preprocessing type is DILU preprocessing and the kernel type is backward substitution in the thread processing of FIG.

[0367] In this case, since face changes in descending order, the calculation processing unit 2621 decides to use threads, sorted, and starts of the row coloring information (step 4501).

[0368] Next, the calculation processing unit 2621 repeats the color loop processing of steps 4502 to 4509. The color loop processing is repeated for c=start color ID to end color ID−1.

[0369] The processing in steps 4502 and 4503 is the same as the processing in steps 4104 and 4105 in Fig. 41. In this case, the kernel type is back-substitution and the matrix used is an upper triangular matrix U ATherefore, the calculation processing unit 2621 determines the start face and the end face according to equations (47) and (48) (step 4504).

[0370] Next, the calculation processing unit 2621 repeats the face loop processing of steps 4505 and 4506. The face loop processing is repeated for f=end face-1 to start face in descending order of face.

[0371] The processing of step 4505 is the same as the processing of step 4110 in Fig. 41. In this case, since the preprocessing type is DILU preprocessing and the kernel type is backward substitution, the calculation processing unit 2621 calculates wA by the following equation (step 4506).

[0372] wA[l[sf]] -=rD[l[sf]]*upper[sf]*wA[u[sf]] (54)

[0373] Equation (54) represents a calculation formula in which f included in the body of the backward substitution loop in the DILU preprocessing shown in Fig. 18 is replaced with sf. The processing in steps 4507 to 4509 is the same as the processing in steps 4112 to 4114 in Fig. 41.

[0374] Figure 46 shows an example timeline of the first DIC forward substitution thread processing using the upper triangular matrix of Figure 20. In this example, the block size is 1.

[0375] (Block, Thread) represents a combination of a block ID and a thread ID. (0,0) represents thread 0 in block 0 of each kernel, and (1,0) represents thread 0 in block 1 of each kernel. Thread 0 in block 0 corresponds to the 0th thread among the number of running threads. Thread 0 in block 1 corresponds to the 1st thread among the number of running threads.

[0376] R represents a read access to upper and wA, and W represents a write access to wA. a to h represent non-zero elements stored in upper, and [s] (s=0 to 5) represents wA[s].

[0377] Kernel K1 corresponds to kernels[0] and includes a single block. Kernel K2 corresponds to kernels[1] and includes two blocks. Read accesses to upper and wA in periods 4601 and 4602 are executed by unrolling the face loop.

[0378] At timing 4611, stream synchronization is performed, and at timing 4612, inter-block synchronization is performed.

[0379] Figure 47 shows an example timeline of the second DIC forward substitution thread processing using the upper triangular matrix of Figure 30. In this example, the block size is 2.

[0380] (0,0) represents thread 0 in block 0 of each kernel, and (0,1) represents thread 1 in block 0 of each kernel. (1,0) represents thread 0 in block 1 of each kernel, and (1,1) represents thread 1 in block 1 of each kernel.

[0381] Thread 0 included in block 0 corresponds to the 0th thread among the number of running threads. Thread 1 included in block 0 corresponds to the 1st thread among the number of running threads. Thread 0 included in block 1 corresponds to the 3rd thread among the number of running threads.

[0382] a to m represent the non-zero elements stored in upper, and [s] (s=0 to 9) represents wA[s].

[0383] Kernel K1 corresponds to kernels[0] and includes a single block. Kernel K2 corresponds to kernels[1] and includes two blocks. Kernel K3 corresponds to kernels[2] and includes a single block. Read accesses to upper and wA in periods 4701 and 4702 are executed by unrolling the face loop.

[0384] At timing 4711 and timing 4713, stream synchronization is performed, at timing 4712, inter-block synchronization is performed, and at timing 4714, intra-block synchronization is performed.

[0385] Fig. 48 shows an example of the hardware configuration of an information processing device including multiple nodes. The information processing device in Fig. 48 includes nodes 4801-1 to 4801-P (P is an integer equal to or greater than 2). Each node 4801-p (p=1 to P) includes a CPU 4811, a memory 4812, and a GPU 4813-1 to GPU 4813-3.

[0386] Each GPU 4813-q (q=1 to 3) includes a processing unit and a storage unit, similar to the GPU 2613 in Fig. 26. The CPU 4811 and the processing unit in the GPU 4813-q correspond to the determination unit 2411 and the processing unit 2412 in Fig. 24, respectively. The memory 4812 stores the same information as in Fig. 27, and the storage unit in the GPU 4813-q stores the same information as in Fig. 28.

[0387] The nodes 4801-1 to 4801-P can communicate with each other via a communication network 4802. Each node 4801-p may include four or more GPUs 4813-q.

[0388] The information processing device in FIG. 48 can divide a sparse matrix in the LDU format and perform analysis processing. In this case, 3P GPUs 4813-q divide the sparse matrix and perform DIC preprocessing or DILU preprocessing. -1The CPU 4811 executes thread processing for calculation, forward substitution, or backward substitution. The CPU 4811 starts a process for each GPU 4813-q to perform analysis processing.

[0389] Each GPU 4813-q does not communicate with other GPUs 4813-q during thread processing. When boundary communication of SpMV, global dot product calculation, etc. occurs, communication occurs between GPUs 4813-q or between nodes 4801-p.

[0390] Figure 49 shows an example of the performance in a fluid dynamics simulation that calculates pressure and velocity for a 200x200x200 cubic grid. In this example, the 3-D Lid Driven cavity flow from the OpenFOAM HPC Benchmark suite is used, with the PCG method used to calculate pressure and the PBiCG method used to calculate velocity.

[0391] The horizontal axis represents block size, and the vertical axis represents processing time (seconds) per time step. Point 4901 represents the performance of THRUST shown in Figure 33. In the case of THRUST, since the block size is unknown, point 4901 is plotted at the position where block size = 128.

[0392] A polygonal line 4902 represents the performance of STREAM, a polygonal line 4903 represents the performance of GRID, and a polygonal line 4904 represents the performance of GRID (unroll=3), a polygonal line 4905 represents the performance of ADAPTIVE, and a polygonal line 4906 represents the performance of ADAPTIVE (unroll=3).

[0393] unroll=3 means that the first three iterations in the face loop are unrolled and executed. For a three-dimensional cubic lattice, unroll=3 is optimal because the coefficient matrix A is a heptadian matrix.

[0394] It can be seen that the processing speed of ADAPTIVE (unroll=3) indicated by the polygonal line 4906 is approximately 1.5 times faster than THRUST indicated by the point 4901, and approximately 1.1 times faster than GRID (unroll=3) indicated by the polygonal line 4904.

[0395] Figure 50 shows examples of the number of preprocessing and smoothers used in the OpenFOAM v8 tutorial. Figure 50(a) shows an example of the number of preprocessing. DIC stands for DIC preprocessing, DILU stands for DILU preprocessing, and GAMG (generalized geometric-algebraic multi-grid) stands for GAMG preprocessing.

[0396] 50(b) shows an example of the number of smoothers, where DIC stands for DIC smoother, DIC GaussSeidel stands for DIC Gauss-Seidel smoother, GaussSeidel stands for Gauss-Seidel smoother, and symGaussSeidel stands for symmetric Gauss-Seidel smoother.

[0397] DIC and DILU preprocessing are frequently used in fields such as fluid dynamics. Therefore, a fast GPU implementation of DIC and DILU preprocessing would be a significant contribution to software solutions in these fields.

[0398] The configurations of the information processing device 2401 in Fig. 24, the information processing device 2601 in Fig. 26, and the information processing device in Fig. 48 are merely examples, and some of the components may be omitted or changed depending on the application or conditions of the information processing device. For example, instead of a GPU, it is also possible to execute parallel processing using another arithmetic processing device such as a CPU. In this case, instead of a thread, another processing unit such as a process may be used.

[0399] The information shown in FIGS. 27 and 28 is merely an example, and some of the information may be omitted or changed depending on the use or conditions of the information processing device.

[0400] The flowcharts in FIGS. 25, 34, 35, 38 to 42, 44, and 45 are merely examples, and some of the processes may be omitted or changed depending on the configuration or conditions of the information processing device.

[0401] The coefficient matrices shown in Figures 1 and 3 are merely examples, and the coefficient matrices change depending on the matrix processing included in the analysis processing. The lower triangular matrices shown in Figures 6 and 9 are merely examples, and the lower triangular matrices change depending on the coefficient matrices used in the analysis processing. The upper triangular matrices shown in Figures 7, 9, 11, 12, 15, 20, and 30 are merely examples, and the upper triangular matrices change depending on the coefficient matrices used in the analysis processing.

[0402] The algorithm of the PCG method shown in Fig. 2 and the algorithm of the PBiCG method shown in Fig. 4 are merely examples, and the PCG method and the PBiCG method can also be written in other formats. The programs of the DILU preprocessing shown in Fig. 5A and Fig. 5B are merely examples, and the programs of the DILU preprocessing can also be written in other formats.

[0403] The DILU preprocessing shown in Figure 8 and the DIC and DILU preprocessing shown in Figures 10 and 18 are merely examples, and the DIC and DILU preprocessing can also be written in other forms. The parallelization shown in Figures 13, 14, 16, 17, and 19 is merely an example, and the parallelization results will vary depending on whether an upper triangular matrix or a lower triangular matrix is ​​used.

[0404] The simulation results shown in Figures 21, 32, and 49 are merely examples, and the simulation results will vary depending on the coefficient matrix used in the simulation. The synchronization methods shown in Figures 22, 23, 29, 31, and 33 are merely examples, and other synchronization methods for threads may be used.

[0405] The color ID update process shown in Figures 36 and 37 is merely an example, and the color ID update process varies depending on the coefficient matrix used in the analysis process. The indexes shown in Figure 43 are merely an example, and the indexes of non-zero elements vary depending on the coefficient matrix used in the analysis process. The timeline shown in Figures 46 and 47 is merely an example, and the timeline of the DIC forward substitution thread process varies depending on the coefficient matrix used in the analysis process.

[0406] The number of preprocessings and smoothers shown in FIG. 50 is just an example, and DIC preprocessing and DILU preprocessing are used in various applications.

[0407] Equations (1) to (54) are merely examples, and matrix processing can also be described using other calculation formulas.

[0408] Fig. 51 shows an example of the hardware configuration of the information processing device 2401 in Fig. 24 and the information processing device 2601 in Fig. 26. The information processing device in Fig. 51 includes a CPU 5101, a memory 5102, an input device 5103, an output device 5104, an auxiliary storage device 5105, a media drive device 5106, a network connection device 5107, and a GPU 5108. These components are hardware and connected to each other via a bus 5109.

[0409] The memory 5102 is, for example, a semiconductor memory such as a ROM (Read Only Memory) or a RAM (Random Access Memory), and stores programs and data used in processing. The memory 5102 may operate as the storage unit 2612 in FIG.

[0410] The CPU 5101 operates as the determination unit 2411 in Fig. 24 by executing a program using the memory 5102. The CPU 5101 also operates as the CPU 2611 in Fig. 26. The CPU 5101 is sometimes called a processor.

[0411] The input device 5103 is, for example, a keyboard, a pointing device, etc., and is used to input instructions or information from a user or operator. The output device 5104 is, for example, a display device, a printer, etc., and is used to output inquiries or instructions to a user or operator and processing results. The processing results may be the results of calculations using solutions for each time step. The output device 5104 may operate as the output unit 2614 in FIG. 26.

[0412] The auxiliary storage device 5105 is, for example, a magnetic disk device, an optical disk device, a magneto-optical disk device, a tape device, etc. The auxiliary storage device 5105 may be a hard disk drive or an SSD (Solid State Drive). The information processing device stores programs and data in the auxiliary storage device 5105 and can use them by loading them into the memory 5102. The auxiliary storage device 5105 may operate as the storage unit 2612 in FIG. 26.

[0413] The medium drive device 5106 drives the portable recording medium 5110 and accesses the recorded contents thereof. The portable recording medium 5110 is a memory device, a flexible disk, an optical disk, a magneto-optical disk, etc. The portable recording medium 5110 may be a CD-ROM (Compact Disk Read Only Memory), a DVD (Digital Versatile Disk), a USB (Universal Serial Bus) memory, etc. A user or operator can store programs and data in the portable recording medium 5110 and load them into the memory 5102 for use.

[0414] In this way, the computer-readable recording medium that stores the program and data used in the processing is a physical (non-transitory) recording medium such as the memory 5102, the auxiliary storage device 5105, or the portable recording medium 5110.

[0415] The network connection device 5107 is a communication interface circuit connected to a communication network such as a LAN (Local Area Network) or a WAN (Wide Area Network) and performs data conversion associated with communication. The information processing device receives programs and data from an external device via the network connection device 5107 and loads them into the memory 5102 for use. The network connection device 5107 may operate as the output unit 2614 in FIG. 26 .

[0416] The GPU 5108 includes a processor 5111 and a memory 5112, performs processing instructed by the CPU 5101, and outputs the processing results to the CPU 5101. The GPU 5108 may operate as the GPU 2613 in Fig. 26. The processor 5111 may operate as the arithmetic processing unit 2412 in Fig. 24 or the arithmetic processing unit 2621 in Fig. 26, and the memory 5112 may operate as the storage unit 2622 in Fig. 26.

[0417] It is not necessary for the information processing device to include all of the components shown in Figure 51, and some components may be omitted depending on the purpose or conditions of the information processing device. For example, if an interface with a user or operator is not required, the input device 5103 and the output device 5104 may be omitted. If the portable recording medium 5110 or a communication network is not used, the medium drive device 5106 or the network connection device 5107 may be omitted.

[0418] Although the disclosed embodiments and their advantages have been described in detail, those skilled in the art may make various modifications, additions, and omissions without departing from the scope of the invention as clearly set forth in the claims.

[0419] The following notes are further provided regarding the embodiment described with reference to FIGS. (Appendix 1) using each of a plurality of processes included in the matrix processing as a first process and a process subsequent to the first process as a second process, and determining a synchronization method for one or more processing units that process the elements of the first part of the matrix in parallel based on the number of one or more processing units that process the elements of the first part of the matrix in parallel in the first process and the number of one or more processing units that process the elements of the second part of the matrix in parallel in the second process; performing the first operation using one or more processing units that process elements of the first portion of the matrix in parallel; performing a synchronization process for one or more processing units that process the elements of the first portion of the matrix in parallel using the synchronization method; performing the second operation using one or more processing units that process the elements of the second portion of the matrix in parallel; An information processing program that causes a computer to execute a process. (Appendix 2) The information processing program described in Appendix 1, characterized in that the process of determining the synchronization method includes a process of determining the synchronization method to stream synchronization when either the number of one or more processing units that process elements of the first part of the matrix in parallel or the number of one or more processing units that process elements of the second part of the matrix in parallel is less than or equal to the group size, and the other number is greater than the group size. (Appendix 3) The information processing program described in Appendix 2, characterized in that the process of executing the synchronization process using the synchronization method includes a process of applying inter-group synchronization as the synchronization method when the number of one or more processing units that process elements of the first part of the matrix in parallel is greater than the group size and the number of one or more processing units that process elements of the second part of the matrix in parallel is greater than the group size. (Appendix 4) The information processing program described in Appendix 2, characterized in that the process of performing the synchronization process using the synchronization method includes a process of applying intra-group synchronization as the synchronization method when the number of one or more processing units that process elements of a first part of the matrix in parallel is less than or equal to the group size and the number of one or more processing units that process elements of a second part of the matrix in parallel is less than or equal to the group size. (Appendix 5) An information processing program according to any one of appendices 1 to 4, characterized in that one or more processing units that process elements of a first part of the matrix in parallel process one or more columns or rows included in the first part of the matrix in parallel, and one or more processing units that process elements of a second part of the matrix in parallel process one or more columns or rows included in the second part of the matrix in parallel. (Appendix 6) The information processing program according to any one of appendices 1 to 5, wherein the matrix is ​​an upper triangular matrix or a lower triangular matrix representing coefficients of simultaneous linear equations, and the matrix processing is included in preprocessing in the solution of the simultaneous linear equations. (Appendix 7) a determination unit that uses each of a plurality of processes included in the matrix processing as a first process and a process subsequent to the first process as a second process, and determines a synchronization method for one or more processing units that process elements of a first portion of the matrix in parallel based on the number of one or more processing units that process elements of a first portion of the matrix in parallel in the first process and the number of one or more processing units that process elements of a second portion of the matrix in parallel in the second process; a processing unit that executes the first processing using one or more processing units that process elements of a first part of the matrix in parallel, executes synchronization processing for the one or more processing units that process elements of the first part of the matrix in parallel using the synchronization method, and executes the second processing using one or more processing units that process elements of a second part of the matrix in parallel; An information processing device comprising: (Appendix 8) The information processing device described in Appendix 7, characterized in that the determination unit determines the synchronization method to be stream synchronization when the number of one or more processing units that process elements of the first part of the matrix in parallel and the number of one or more processing units that process elements of the second part of the matrix in parallel is less than or equal to the group size, and the number of the other is greater than the group size. (Appendix 9) The information processing device described in Appendix 8, characterized in that the calculation processing unit applies inter-group synchronization as the synchronization method when the number of one or more processing units that process elements of the first part of the matrix in parallel is greater than the group size and the number of one or more processing units that process elements of the second part of the matrix in parallel is greater than the group size. (Appendix 10) The information processing device described in Appendix 8, characterized in that the calculation processing unit applies intra-group synchronization as the synchronization method when the number of one or more processing units that process elements of the first part of the matrix in parallel is less than or equal to the group size and the number of one or more processing units that process elements of the second part of the matrix in parallel is less than or equal to the group size. (Appendix 11) An information processing device as described in any one of Supplementary Notes 7 to 10, characterized in that one or more processing units that process elements of the first part of the matrix in parallel process one or more columns or rows included in the first part of the matrix in parallel, and one or more processing units that process elements of the second part of the matrix in parallel process one or more columns or rows included in the second part of the matrix in parallel. (Appendix 12) The information processing device described in any one of appendices 7 to 11, wherein the matrix is ​​an upper triangular matrix or a lower triangular matrix representing coefficients of a system of linear equations, and the matrix processing is included in preprocessing in the solution of the system of linear equations. (Appendix 13) using each of a plurality of processes included in the matrix processing as a first process and a process subsequent to the first process as a second process, and determining a synchronization method for one or more processing units that process the elements of the first part of the matrix in parallel based on the number of one or more processing units that process the elements of the first part of the matrix in parallel in the first process and the number of one or more processing units that process the elements of the second part of the matrix in parallel in the second process; performing the first operation using one or more processing units that process elements of the first portion of the matrix in parallel; performing a synchronization process for one or more processing units that process the elements of the first portion of the matrix in parallel using the synchronization method; performing the second operation using one or more processing units that process the elements of the second portion of the matrix in parallel; An information processing method characterized in that the processing is executed by a computer. (Appendix 14) The information processing method described in Appendix 13, characterized in that the process of determining the synchronization method includes a process of determining the synchronization method to stream synchronization when the number of one or more processing units that process elements of the first part of the matrix in parallel and the number of one or more processing units that process elements of the second part of the matrix in parallel is less than or equal to the group size, and the number of the other is greater than the group size. (Appendix 15) The information processing method of claim 14, wherein the process of performing the synchronization process using the synchronization method includes a process of applying inter-group synchronization as the synchronization method when the number of one or more processing units that process elements of a first part of the matrix in parallel is greater than the group size and the number of one or more processing units that process elements of a second part of the matrix in parallel is greater than the group size. (Appendix 16) The information processing method of claim 14, wherein the process of performing the synchronization process using the synchronization method includes a process of applying intra-group synchronization as the synchronization method when the number of one or more processing units that process elements of a first part of the matrix in parallel is less than or equal to the group size and the number of one or more processing units that process elements of a second part of the matrix in parallel is less than or equal to the group size. (Appendix 17) An information processing method according to any one of appendices 13 to 16, characterized in that one or more processing units that process elements of the first part of the matrix in parallel process one or more columns or rows included in the first part of the matrix in parallel, and one or more processing units that process elements of the second part of the matrix in parallel process one or more columns or rows included in the second part of the matrix in parallel. (Appendix 18) The information processing method described in any one of appendices 13 to 17, wherein the matrix is ​​an upper triangular matrix or a lower triangular matrix representing coefficients of simultaneous linear equations, and the matrix processing is included in preprocessing in the solution of the simultaneous linear equations. [Explanation of symbols]

[0420] 101-104, 111-115, 121-124, 131-134, 141-144, 151-154 elements 202 if statements 501~503 sentences 504, 506, 508, 1311, 1312, 1611, 1612 for loop 505, 507, 509, 1711~1714 formula 901~904 Triangle 911~914 Arrow 1201, 1501 directed graphs 12:11~12:14, 15:11~15:14 Vertex 1301, 1302, 1601, 1602, 1901, 1902 Loop processing 1401~1403, 1701~1703 Processing 2011, 2012-1, 2012-2, 2013-1, 2013-2, 2311, 2312-1~2312-3, 2313-1~ 2313-3, 2314, 2315, 3111, 3112-1~3112-3, 3113-1~3113-3, 3114, 3115 calculation 2021, 2022, 2211, 2212, 2221, 2321-2324, 2331-2334 Synchronous processing 2201~2204, 2301~2307, 2901, 2902, 3101~3103 kernels 2401, 2601 Information processing equipment 2411 Decision Section 2412, 2621 Processing unit 2611, 4811, 5101 CPUs 2612, 2622 Storage section 2613, 4813-1 to 4813-3 GPUs 2614 Output section 2711 CPU Coloring Information 2811 Sparse Matrix Information 2812 GPU coloring information 2813 Iteration Method Objects 2911, 3121 Inter-block synchronization 3122 Intra-block synchronization 3201~3203 Rectangle 4601, 4602, 4701, 4702 period 4611, 4612, 4711~4714 Timing Nodes 4801-1 to 4801-P 4802 Communication Networks 4901 points 4902~4906 Line 4812, 5102, 5112 memory 5103 Input Device 5104 Output device 5105 Auxiliary storage device 5106 Media drive unit 5107 Network connection device 5109 Bus 5110 Portable recording media 5111 processor

Claims

1. using each of a plurality of processes included in the matrix processing as a first process and a process subsequent to the first process as a second process, and determining a synchronization method for one or more processing units that process the elements of the first part of the matrix in parallel based on the number of one or more processing units that process the elements of the first part of the matrix in parallel in the first process and the number of one or more processing units that process the elements of the second part of the matrix in parallel in the second process; performing the first operation using one or more processing units that process elements of a first portion of the matrix in parallel; performing synchronization of one or more processing units that process the elements of the first portion of the matrix in parallel using the synchronization method; performing the second operation using one or more processing units that process elements of a second portion of the matrix in parallel; An information processing program that causes a computer to execute a process.

2. The information processing program according to claim 1, characterized in that the process of determining the synchronization method includes a process of determining the synchronization method to stream synchronization when either the number of one or more processing units that process elements of the first part of the matrix in parallel or the number of one or more processing units that process elements of the second part of the matrix in parallel is less than or equal to the block size, and the other number is greater than the block size.

3. The information processing program according to claim 2, characterized in that the process of executing the synchronization process using the synchronization method includes a process of applying inter-block synchronization as the synchronization method when the number of one or more processing units that process elements of a first part of the matrix in parallel is greater than the block size and the number of one or more processing units that process elements of a second part of the matrix in parallel is greater than the block size.

4. The information processing program according to claim 2, characterized in that the process of executing the synchronization process using the synchronization method includes a process of applying intra-block synchronization as the synchronization method when the number of one or more processing units that process elements of the first part of the matrix in parallel is less than or equal to the block size and the number of one or more processing units that process elements of the second part of the matrix in parallel is less than or equal to the block size.

5. An information processing program according to any one of claims 1 to 4, characterized in that one or more processing units that process elements of the first part of the matrix in parallel process one or more columns or rows included in the first part of the matrix in parallel, and one or more processing units that process elements of the second part of the matrix in parallel process one or more columns or rows included in the second part of the matrix in parallel.

6. 6. The information processing program according to claim 1, wherein the matrix is ​​an upper triangular matrix or a lower triangular matrix representing coefficients of simultaneous linear equations, and the matrix processing is included in preprocessing in a method for solving the simultaneous linear equations.

7. a determination unit that uses each of a plurality of processes included in the matrix processing as a first process and a process subsequent to the first process as a second process, and determines a synchronization method for one or more processing units that process elements of a first portion of the matrix in parallel based on the number of one or more processing units that process elements of a first portion of the matrix in parallel in the first process and the number of one or more processing units that process elements of a second portion of the matrix in parallel in the second process; performing the first operation using one or more processing units that process elements of the first portion of the matrix in parallel; a processor that executes a synchronization process using the synchronization method, the second process being performed using one or more processing units that process elements of the second part of the matrix in parallel; An information processing device comprising:

8. using each of a plurality of processes included in the matrix processing as a first process and a process subsequent to the first process as a second process, and determining a synchronization method for one or more processing units that process the elements of the first part of the matrix in parallel based on the number of one or more processing units that process the elements of the first part of the matrix in parallel in the first process and the number of one or more processing units that process the elements of the second part of the matrix in parallel in the second process; performing the first operation using one or more processing units that process elements of a first portion of the matrix in parallel; performing synchronization of one or more processing units that process the elements of the first portion of the matrix in parallel using the synchronization method; performing the second operation using one or more processing units that process elements of a second portion of the matrix in parallel; An information processing method characterized in that the processing is executed by a computer.

Images