JP2011145999A - Method and apparatus for calculating simultaneous linear equation - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a preconditioned iterative solution that has a highly accelerated effect of preconditioning, and a practical preconditioned iterative solution that does not leave parameter determination to a user by providing the step of determining appropriate parameters depending on a coefficient matrix. <P>SOLUTION: A coefficient matrix A is incompletely triangular-factorized into a lower triangular matrix L and an upper triangular matrix U, and diagonal elements and off-diagonal elements of L and U are multiplied by respective scaling parameters to form a lower triangular matrix L' and an upper triangular matrix U'. The preconditioned iterative solution is executed by the use of a preconditioned matrix M given by L' and U'. The step of scaling parameter determination estimates such scaling parameters that minimize an evaluation function defined by the coefficient matrix A and preconditioned matrix M. The evaluation function is the Euclidean norm of a vector expressed by a formula whose component is the sum of columns of a difference matrix between the coefficient matrix A and the preconditioned matrix M. <P>COPYRIGHT: (C)2011,JPO&INPIT

Description

  The present invention relates to a method for calculating simultaneous linear equations, and more particularly to a preprocessed iterative solution of simultaneous linear equations whose coefficient matrix is a large-scale sparse matrix.

  Various engineering problems such as the analysis of physical phenomena are ultimately reduced to a multi-dimensional simultaneous linear equation by discretizing a partial differential equation expressed in a simulated manner. In many cases, the approximation accuracy of discretization is at most several orders, so the coefficients of the equations are sparse matrices. On the other hand, in order to solve the problem with high accuracy, it is necessary to increase the number of discretization elements, and thus the coefficient matrix is often a large matrix. Therefore, it is very important to solve simultaneous linear equations having a large-scale sparse matrix as a coefficient matrix.

  The iterative solving method is one of solving methods of simultaneous linear equations, and is known as the most effective method especially when the coefficient matrix is a large-scale sparse matrix. A specific algorithm is to first give an appropriate initial value as an approximate solution, and repeat the update process of the approximate solution to obtain a more accurate solution, thereby gradually approaching the true solution. For example, the CG method (conjugate gradient method) specialized for positive definite symmetric matrices and the BiCGStab method (stabilized biconjugate gradient method) corresponding to asymmetric matrices are often used.

  In the iterative solution method, the convergence point, that is, how close the approximate solution can be to be sufficiently close to the true solution with a small number of iterations is an important point. This convergence is greatly influenced by the condition number of the coefficient matrix, that is, the value obtained by dividing the maximum value of the eigenvalues by the minimum value. The smaller the condition number, the better the convergence. Based on this property, various preprocessing methods have been proposed as means for speeding up the iterative solution.

  The basic idea of preprocessing is to prepare a preprocessing matrix that approximates the coefficient matrix of the equation to be solved, and multiply both sides of the equation by the inverse matrix, thereby reconstructing the equation itself with a coefficient matrix with a smaller condition number In addition, the convergence in the iterative solution is improved. Actually, the reconstruction of the equation is not performed explicitly because of the heavy load, and the reconstructed equation is incorporated by incorporating a vector correction process (hereinafter referred to as preprocessing) by a preprocessing matrix in the iterative routine of the iterative solution method. Construct an algorithm that is mathematically equivalent to the iterative solution to.

  There are two properties required for the preprocessing matrix. One is that the closeness to the coefficient matrix is high. The higher the approximation, the closer the coefficient matrix to the unit matrix by the preprocessing and the smaller the condition number, so the convergence of the iterative solution improves. The other is that the time required for creating and preprocessing the preprocessing matrix is small. This is because even if the convergence is improved by the preprocessing, if these processes take too long, the calculation time may increase as a whole solution. In general, the two are in a trade-off relationship, and various preprocessing methods have been proposed in order to achieve a good balance.

Among them, one of the most popular methods is incomplete trigonometric decomposition preprocessing (see Non-Patent Document 1). This method is a preprocessing method in which triangulation of a coefficient matrix is incomplete and the result is used as a preprocessing matrix. The details will be described below.
Triangular decomposition is a generic term for LU decomposition that decomposes a matrix into a product of a lower triangular matrix L and an upper triangular matrix U, and LDU decomposition that decomposes into a product of a lower triangular matrix L, a diagonal matrix D, and an upper triangular matrix U. It is. Triangular decomposition is a kind of solution method for simultaneous linear equations, which is essentially called direct solution, but is not usually used when the coefficient matrix is a large-scale sparse matrix. This is because there are many f
This is because ill-in (fill-in) occurs and a large amount of calculation time and memory are required when the matrix is large. Here, fill-in means that a new non-zero element is written to a position that was zero in the original matrix during decomposition.
The incomplete triangulation preprocessing is a preprocessing method that uses, as a preprocessing matrix, a matrix that is incompletely decomposed by rejecting the fill-in in the triangulation of the coefficient matrix. The closeness of the preprocessing matrix to the coefficient matrix is higher than, for example, diagonal scaling preprocessing using a matrix obtained by extracting only the diagonal components of the coefficient matrix as a preprocessing matrix, and shows better convergence. On the other hand, it takes some time to create the preprocessing matrix, but it can be done in a relatively short time because the fill-in is rejected. The same applies to pre-processing. Since the pre-processing matrix is decomposed into L * U or L * D * U, simple operations such as forward substitution and backward substitution are sufficient, and L and U are coefficient matrices. Since the sparseness is inherited, the amount of calculation is small, so this can also be executed in a short time. These features enable efficient iterative solution speedup and are used as one of the most popular preprocessing methods.

Several derivation methods have been proposed for incomplete triangulation preprocessing in order to further improve convergence. Both are characterized in that some kind of contrivance is applied when the preprocessing matrix is created.
Both the incomplete triangulation preprocessing with shift parameters (see Non-Patent Document 1) and the incomplete triangulation preprocessing with correction parameters (see Non-Patent Document 2) are both performed at the time of creating the preprocessing matrix. Incomplete triangulation is performed after the change in the diagonal term. The former adds a shift parameter to the diagonal term, and the latter multiplies the diagonal term by a correction parameter. In particular, when the diagonal advantage of the coefficient matrix is low, the objective is to prevent the diagonal term from approaching zero due to the decomposition and to avoid deterioration of convergence.
The modified incomplete triangulation preprocessing (see Non-Patent Document 1) does not simply reject the fill-in in the incomplete triangulation, but subtracts this amount from the diagonal term of the line where the fill-in occurs. Features. This is an operation for forcing the preprocessing matrix to have the same row sum as the original matrix. Actually, the fill-in is often multiplied by a relaxation parameter and then subtracted from the diagonal term, and it is known that the convergence is improved depending on the relaxation parameter.
The pre-processing of incomplete triangulation with a fill-in level (see Non-Patent Document 1) is a technique that does not completely reject the fill-in at the time of incomplete triangulation but allows a fill-in below a certain level. is there. Here, the fill-in level means that a position that is non-zero in the original coefficient matrix is set to level 0, and a fill-in that refers to an element at k level by decomposition is defined as level k + 1. . As the permissible level is increased, the “incompleteness” of the incomplete decomposition is reduced, so that the approximation to the coefficient matrix is increased and the convergence is improved.

J.J. Dongara et al./Satomi Hasegawa Hidehiko Hasegawa Kiyoji Fujino, “Iteration Method Templates”, Asakura Shoten, 1996, pp. 53-57 Isao Kataoka Hideshi Yasuda Naoki Takano, Masahiko Shibahara, “Introduction to Numerical Analysis”, Corona, 2002, pp. 57

  The background art described above has two problems. One is that the effect of speeding up by preprocessing is low, and the other is that each derivative method does not have a step of determining an appropriate value for each unique parameter.

The reason why the effect of speeding up by preprocessing, which is the first problem, is low will be described in detail below. First, in the basic incomplete triangulation, the triangulation of the coefficient matrix is interrupted in an incomplete state and used as it is as a preprocessing matrix. However, there is no guarantee that the preprocessing matrix immediately after decomposition has the best approximation to the coefficient matrix within the range satisfying the desired sparseness, and the improvement in convergence by the preprocessing is limited. For other derivation methods, the basic incomplete triangulation preprocessing is corrected before and during incomplete decomposition, and the matrix immediately after the incomplete decomposition is used as the preprocessing matrix. However, the improvement in convergence is limited.

  In addition, each derivative method has the following factors that impede the speed-up effect. First, for incomplete triangulation preprocessing with shift parameters and correction parameters, changing the diagonal terms of the coefficient matrix before decomposition is a direct means of improving the approximation of the preprocessing matrix to the coefficient matrix. Therefore, the improvement in convergence is limited. In addition, with respect to the modified incomplete triangulation, since the operation of calculating fill-in which is not originally required for calculation and subtracting from the diagonal term is performed, it takes a lot of time to create the preprocessing matrix. Concerning incomplete triangulation preprocessing with a fill-in level, the convergence is further improved as the fill-in permissible amount is increased, but at the cost of the preprocessing matrix creation and the calculation time required for the preprocessing increases. .

  As described above, with respect to all the conventional methods described above, there is a problem that the speed-up effect as a whole solution is low, and a preprocessed iterative solution method with a higher speed-up effect is desired.

  Next, the point which does not have the process of determining appropriately the intrinsic | native parameter in each derivation method which is the 2nd subject is demonstrated in detail below. All derivations shown in the background art have their own parameters. In any case, each parameter greatly affects the degree of speedup, but care must be taken because appropriate values differ for each coefficient matrix. First, for incomplete triangulation preprocessing with shift parameters and correction parameters, if the parameters are too far from 0 and 1, respectively, the matrix to be decomposed is far from the original coefficient matrix, and the approximation of the preprocessing matrix decreases. This leads to deterioration of convergence. In addition, it is known that the modified incomplete triangulation decomposition preconditioning also causes the diagonal term to approach zero due to the decomposition and deteriorates the convergence if the relaxation parameter is too large. Concerning the incomplete triangulation with the fill-in level, the convergence improves as the fill-in tolerance level increases, but conversely, the creation of the preprocessing matrix and the preprocessing take time. As described above, regarding these derivation methods, it is important to set appropriate parameters according to the coefficient matrix. However, it is not clear how to estimate an appropriate value in a time that does not impair the effect of speeding up by preprocessing. As a result, these derivations do not have the process of determining appropriate parameters, and parameter determination is left to the user. Therefore, in order to surely obtain the effect of speeding up by preprocessing, there is no choice but to prepare several appropriate parameters for the coefficient matrix, repeat the trial, and manually search for the appropriate parameters.

  In numerical analysis methods such as the finite difference method and the finite element method, simultaneous linear equations having different coefficient matrices need to be solved one after another for each short time step. When applying each derivation method of the background art to this, it is unrealistic in terms of calculation time to perform manual search of the above parameters at each time step, and it is only possible to perform calculation with the parameters fixed to appropriate values. Absent. However, in that case, the effect of the preprocessing varies depending on the time step, and in the worst case, convergence is not obtained and the calculation may be stopped, which is dangerous. Therefore, a preprocessed iterative solution having a step of determining appropriate parameters according to the coefficient matrix is desired.

  The present invention has been made under such a background, and a first object of the present invention is to provide a numerical calculation method and apparatus for simultaneous linear equations that are highly effective in speeding up by preprocessing. In addition to the first object, the second object is to provide a numerical calculation method and apparatus for simultaneous linear equations capable of determining an appropriate parameter according to a coefficient matrix.

  In order to achieve the first object, the present invention is a numerical solution of simultaneous linear equations, in which a coefficient matrix A is incompletely triangulated to obtain a lower triangular matrix L and an upper triangular matrix U; A scaling step of obtaining a lower triangular matrix L ′ and an upper triangular matrix U ′ by multiplying the diagonal elements of L and U decomposed in the decomposing step by a scaling parameter γ and multiplying the non-diagonal elements by a scaling parameter φ. A preprocessing step of performing preprocessing using the preprocessing matrix M obtained from the L ′ and the U ′, and an approximate solution updating step of updating the approximate solution using the result obtained in the preprocessing step; , Determining the convergence of the solution according to the convergence determination condition, and if it has converged, the process is terminated, and if not, an iterative process that repeats the preprocessing process and the approximate solution update process is included. .

のユークリッドノルムを採用することが好適である。 In order to achieve the second object, the present invention estimates the parameters that minimize the evaluation function defined by the coefficient matrix A and the preprocessing matrix M for the scaling parameters γ and φ. The method further includes a parameter determining step that is determined by: As this evaluation function, a vector whose component is the sum of each difference matrix of the coefficient matrix A and the preprocessing matrix M

It is preferable to adopt the Euclidean norm.

  According to the present invention, by multiplying the diagonal and non-diagonal terms of the lower triangular matrix L and the upper triangular matrix U by the scaling parameter in the scaling step, A direct correction for improving the closeness to the coefficient matrix is possible. Therefore, improvement in convergence can be expected. On the other hand, the amount of calculation in the scaling step is very small, and the sparseness of the preprocessing matrix is maintained even after the scaling step. Therefore, the trade-off factor for improving the convergence, that is, the time required for creating the preprocessing matrix and the preprocessing hardly increases, and the effect of speeding up by the preprocessing is enhanced.

  According to the present invention, in the parameter determination step, an appropriate scaling parameter corresponding to the coefficient matrix can be estimated in a time that does not impair the effect of the preprocessing.

  Furthermore, by applying the present invention to several derivation methods of the background art, the sensitivity to the parameters unique to the derivation method can be reduced and the practicality can be improved.

Flow chart showing the numerical solution of the present invention The figure which shows the structure of the numerical calculation apparatus of this invention Preconditioned conjugate gradient method (CG method) algorithm Relationship between scaling parameter and number of iterations / calculation time when applying the present invention Comparison between the conventional method and the present invention Effects of applying the present invention to the conventional method Effects of applying the present invention to the conventional method

FIG. 2 is a block diagram showing a basic configuration of a computer that functions as a numerical calculation apparatus according to the present embodiment. The computer generally includes a CPU 201, a display device 202, an input device 203, a RAM 204, a hard disk device 205, an NI connected via a bus 207.
C206 and the like.

  A CPU (central processing unit) 201 controls the entire apparatus using programs and data stored in a RAM (main storage device) 204 and executes each process described later. The display device 202 is configured by a liquid crystal screen or the like, and can display the processing result by the CPU 201 with characters as images. The input device 203 is configured by an operation device such as a keyboard and a mouse, and can input various instructions to the CPU 201. The RAM 204 includes an area for temporarily storing programs and data loaded from the hard disk device 205 and a work area required when the CPU 201 executes various processes. The hard disk device 205 stores an OS (operating system) and programs and data for causing the CPU 201 to execute processes to be described later that should be performed by the computer. These programs are loaded into the RAM 204 and executed by the CPU 201, whereby the computer functions as a numerical calculation device for numerically solving the simultaneous linear equations. The NIC (network interface) 206 functions as an interface when performing data communication with an external device, and the computer transmits and receives programs and data to and from the external device via the NIC 206. Note that instead of the programs and data stored in the hard disk device 205, the programs and data received from the external device via the NIC 206 may be processed by the CPU 201.

  Next, a numerical solution method for obtaining a solution of simultaneous linear equations performed by a computer having the above configuration will be described. FIG. 1 is a flowchart of numerical calculation processing for obtaining a solution of simultaneous linear equations performed by the computer according to the present embodiment. Here, the iterative solution method with preprocessing of the present invention is applied to the solution processing of the simultaneous linear equations represented by Ax = b as the coefficient matrix A, the variable vector x, and the right side vector b.

  First, the coefficient matrix A is incompletely triangularly decomposed to obtain a lower triangular matrix L and an upper triangular matrix U (step S101: decomposition step). Here, ILDU decomposition will be described as an example of incomplete triangulation, but this does not limit the present invention, and the incomplete triangulation may be ILU decomposition or IC (incomplete Cholesky) decomposition. Good. In addition, any of the incomplete triangulation preprocessing with shift parameters, the incomplete triangulation preprocessing with correction parameters, the corrected incomplete triangulation preprocessing, and the incomplete triangulation preprocessing with the fill-in level described in Background Art It may be incomplete triangulation processing in the derivation method. That is, the decomposition process here may be any process that incompletely decomposes the coefficient matrix and derives the matrix product of L * U or L * D * U.

Then, the off-diagonal of the diagonal elements L and scaling parameters γ multiplying the diagonal elements U _ii of _ii and upper triangular matrix U, the off-diagonal elements L _{ij (i>} j) of the L and U of the lower triangular matrix L A scaling parameter φ to be multiplied by the element U _ij (i <j) is determined (step S102). The scaling parameters γ and φ are single constants that are uniformly applied to all elements belonging to the diagonal and non-diagonal terms of L and U, respectively. Note that i and j represent subscripts in the row direction and column direction of the matrix. In this parameter determination method, it is desirable that the time required for the determination is such that the speed-up effect by the preprocessing is not impaired. In addition, it is desirable to have a method capable of deriving appropriate values such that the approximation of the preprocessing matrix to the coefficient matrix is improved by applying parameters. Therefore, the scaling parameters φ and γ are determined as values that minimize the evaluation function defined by the coefficient matrix A and the preprocessing matrix M. This evaluation function is preferably a function whose value decreases as the closeness of the coefficient matrix A and the preprocessing matrix M increases. In the following, a method for estimating a parameter that minimizes the evaluation function using the Euclidean norm of the following vector R whose component is the sum of the difference matrix between the coefficient matrix A and the preprocessing matrix M as a component will be described. . That is, the evaluation function is the Euclidean norm of the vector R expressed by the following equation. Note that n represents the order of the matrix.

  The minimization of the evaluation function may be performed by numerically solving γ and φ that make the first derivative of the evaluation function zero by a method such as Newton's method. However, this does not limit the present invention, and the optimization of the multivariable function is not limited. Any of the methods for making it possible may be used.

Next, using the scaling parameters γ _opt and φ _opt obtained in the parameter determination step, L and U are scaled to obtain a lower triangular matrix L ′ and an upper triangular matrix U ′ (step S103). That is, as expressed by the following equation, the scaling parameter γ _opt is multiplied by the diagonal element L _ii of L and the diagonal element U _ii of U, and the scaling parameter φ _opt is _converted to the non-diagonal element L _{ij of} L (i> j) and U off-diagonal elements U _ij (i <j).

  The preprocessing matrix M = L ′ * D * U ′ is determined by the operations of steps S101 to S103 described above. Using this preprocessing matrix M, the iterative solution with preprocessing is performed in subsequent steps S104 to S106. Here, any other iterative solution with pre-processing may be used, such as a steady iterative solution method such as the Richardson iterative method, a non-stationary iterative solution method such as the CG method, BiCG method, CGS method, or BiCGStab method. Hereinafter, the CG method with preprocessing will be described as an example. Since the CG method with preprocessing is well known, detailed description is omitted, but an example of the mounting method is shown in FIG. However, it may be calculated by other methods that are mathematically equivalent.

なお、前進代入・後退代入の詳細については公知の技術であるためここでは記載しない。 First, preprocessing is performed on the residual vector r in the CG method (step S104). That is, a vector z expressed by the following equation is obtained using forward substitution / backward substitution.

The details of forward substitution and backward substitution are well-known techniques, and are not described here.

  Next, the approximate solution in the preprocessed CG method is updated using the vector z obtained in step S104 (step S105). Note that since the approximate solution update of the preprocessed CG method is a known technique, a detailed description thereof will be omitted, but the approximate solution is updated by performing the approximate solution update process shown in FIG.

Next, the convergence determination of the approximate solution is performed. If it has converged, the process is terminated, and if not, the process returns to step S104 (step S106). Here, for example, when the convergence determination condition that the value obtained by dividing the Euclidean norm of the residual vector r by the Euclidean norm of the right-hand vector b is equal to or less than the threshold for convergence determination is satisfied, it is determined that it is not converged. To do. However, this does not limit the present invention, and any method may be used as long as it is a convergence determination of the iterative solution with preprocessing.

  In the above description, the case where the preprocessing matrix M is configured by L * U or L * D * U has been described. However, this is multiplied by an arbitrary scalar value a, and M = a * L * U or a *. You may comprise as L * D * U. For example, the reciprocal of the scaling parameter γ may be used for a.

また、スケーリングパラメタγ，φをこのような評価関数を最小とするものとして決定することが好適であるが、利用者が適当なパラメタを用意して試行を繰り返して適切なパラメタを手動で探索しても良い。 In the above description, the Euclidean norm of the vector R whose component is the sum for each column of the difference matrix is used as the evaluation function, and the scaling parameters γ and φ that minimize this are used. It may be used. For example, the Euclidean norm of a vector whose component is the sum of each row of the difference matrix may be used as the evaluation function. In that case, the evaluation function is the Euclidean norm of the vector R expressed by the following equation. Note that n represents the order of the matrix.

In addition, it is preferable to determine the scaling parameters γ and φ as those that minimize such an evaluation function. However, the user prepares an appropriate parameter and repeats the trial to manually search for the appropriate parameter. May be.

  Also, in order to make the scaling parameter determination process simple and fast, the scaling parameter φ may be fixed at 1.0, and the scaling parameter γ that minimizes the evaluation function under this condition may be determined. . In this case, in the scaling step, the scaling parameters γ and φ can be regarded as being multiplied by the diagonal elements and non-diagonal elements of the triangular matrices L and U after incomplete triangulation, or only the scaling parameter γ is triangular. It can also be understood that the diagonal elements of the matrices L and U are multiplied.

  The effect of the iterative method with preprocessing according to the embodiment of the present invention described above will be described in detail below. Note that the data at a certain time step when the three-dimensional Poisson equation appearing in the fluid analysis was solved by the finite difference method was used for the simultaneous linear equations to be solved. The coefficient matrix is a positive definite symmetric matrix with a scale of 540,000 yuan. In the present invention, the scaling parameter φ of the non-diagonal element is fixed at 1.0. Note that ILDU decomposition is used for incomplete triangulation as a base, and CG method is used for iterative solution.

  First, in FIG. 4, the result when the value of the scaling parameter γ of the diagonal element is manually changed and the result when the scaling parameter γ is determined so as to minimize the evaluation function (Euclidean norm of the vector R) Compare 4A is a graph showing the number of iterations, and FIG. 4B is a graph showing calculation time. In the figure, white circles and white triangles change γ manually, and black circles and black triangles set γ so that the evaluation function is minimized.

As a result of manually changing γ, both the number of iterations and the calculation time showed minimum values at γ = 0.63. Compared to the basic CG method with ILDU decomposition preprocessing, γ = 1.0, the number of iterations is reduced from 171 times to 95 times, and the convergence is greatly improved. On the other hand, the calculation time is also reduced from 12.7 seconds to 7.3 seconds, and it can be seen that the time required for creating the preprocessing matrix and preprocessing hardly increases. As a result, a significant speed increase of about 1.7 times is obtained.
On the other hand, when γ was determined so as to minimize the evaluation function, γ = 0.64, and a value close to the optimum value by manual search was estimated. As a result, the number of iterations was 96 and the calculation time was 7.4 seconds, which was inferior to manual parameter optimization. This result shows that the above parameter determination method estimates an appropriate scaling parameter according to the coefficient matrix in a time that does not impair the speed-up effect of the present method.

  Next, in FIG. 5, various conventional methods described in the background art (incomplete triangulation pre-processing, its derived method with shift parameter, with correction parameter, incomplete triangulation pre-processing with fill-in level, uncorrected The complete decomposition pretreatment) is compared with the method of the present invention. In FIG. 5A, the number of iterations is compared, and in FIG. 5B, the calculation time is compared. In the present invention, the scaling parameter is determined by the parameter determination method described above. On the other hand, with respect to the four derivation methods of the background art, the optimum values that give the best results by trying various values for each unique parameter are searched in advance, and the results using the parameters are shown. From the figure, the results according to the present invention are the best in both the number of iterations and the calculation time. In other words, it can be seen that the solution of the present invention shows superior performance and is more practical than the respective derivation methods of the background art in which the parameters are optimized through repeated trials in a single calculation.

  Note that the method of the present invention in which the lower triangular matrix L and the upper triangular matrix U are multiplied by the scaling parameter can be applied simultaneously to each derivation method of the incomplete triangulation preprocessing. FIG. 6 shows the results when the method of the present invention is applied to each derivation method. That is, not only the basic ILDU decomposition but also the original ILDU decomposition performed in each derivation method is used for the incomplete triangulation in the present invention. A scaling parameter that minimizes the evaluation function by the parameter determination process is adopted.

  FIG. 6A shows the number of iterations, and FIG. 6B shows the calculation time. Similarly to FIG. 5, the four derivation methods of the prior art show the result of optimizing the eigenparameters in advance. From the figure, it can be seen that not only the basic ILDU decomposition preprocessing but also all four derivation methods can be speeded up by applying the present invention. In other words, the present invention does not compete with any prior art, and has the property of enhancing the effect when combined with them.

  FIG. 6 shows only the result of optimizing the eigenparameter for each derivation method, but FIG. 7 shows the result of trying various values for the eigenparameter. FIGS. 7A and 7B show the CG method with ILDU decomposition preprocessing with shift parameters and the results when the present invention is applied thereto. Thereafter, similarly, FIGS. 7C and 7D are ILDU decompositions with correction parameters, FIGS. 7E and 7F are modified ILDU decompositions, and FIGS. 7G and 7H are ILDU decompositions with a fill-in level. . In the figure, white circles and white triangles indicate the number of iterations and calculation time in the solution only by the conventional method, and black circles and black triangles indicate the number of iterations and the calculation time when the present invention and the conventional method are combined.

  First of all, for all derivation methods, it can be seen that the sensitivity to unique parameters (shift parameter, correction parameter relaxation parameter, fill-in level) is very high, and the parameter setting greatly affects the degree of speedup. On the other hand, when the present invention is applied to each derivation method, both the number of iterations and the calculation time are reduced regardless of the value of the eigenparameter, and robustness with respect to the eigenparameter appears.

  In particular, by applying the present invention to ILDU decomposition pre-processing with correction parameters (FIGS. 7C and 7D) and modified ILDU decomposition pre-processing (FIGS. 7E and 7F), the correction parameters and It can be seen that the sensitivity to the relaxation parameter is reduced. As a result, a certain range of parameters are uniformly optimal values, which increases the safety when using fixed parameters. In other words, it can be seen that the solution of the present invention provides two effects of speeding up and improving practicality with respect to the above two derivative methods.

  201 ... CPU, 202 ... display device, 203 ... input device, 204 ... RAM

Claims (6)

A method for calculating simultaneous linear equations in a numerical computer,
The numerical calculator is
A decomposition step of obtaining a lower triangular matrix L and an upper triangular matrix U by incomplete triangular decomposition of the coefficient matrix A;
A scaling step of obtaining a lower triangular matrix L ′ and an upper triangular matrix U ′ by multiplying the diagonal elements of L and U decomposed in the decomposing step by a scaling parameter γ and multiplying the non-diagonal elements by a scaling parameter φ. When,
A preprocessing step of performing preprocessing using a preprocessing matrix M obtained from L ′ and U ′;
An approximate solution update step of updating the approximate solution using the result obtained in the preprocessing step;
It determines the convergence of the solution according to the convergence determination condition, and if it has converged, terminates the process; otherwise, repeats the preprocessing step and the approximate solution update step; and
A method of calculating simultaneous linear equations, characterized in that   2. The parameter determining step of determining the scaling parameters γ and φ by estimating a parameter that minimizes an evaluation function defined by the coefficient matrix A and the preprocessing matrix M. Calculation method of simultaneous linear equations described in 1. The evaluation function is a vector whose component is a sum for each column of a difference matrix between the coefficient matrix A and the preprocessing matrix M.

The simultaneous linear equation calculation method according to claim 2, wherein the Euclidean norm is A numerical calculation device for finding a solution of simultaneous linear equations,
Decomposition means for obtaining a lower triangular matrix L and an upper triangular matrix U by incomplete triangular decomposition of the coefficient matrix A;
Scaling means for obtaining a lower triangular matrix L ′ and an upper triangular matrix U ′ by multiplying the diagonal elements of L and U decomposed by the decomposing means by a scaling parameter γ and multiplying non-diagonal elements by a scaling parameter φ. When,
Preprocessing means for performing preprocessing using a preprocessing matrix M obtained from the L ′ and the U ′;
Approximate solution update means for updating the approximate solution using the result of the preprocessing means;
It determines the convergence of the solution according to the convergence determination condition, ends the process if it has converged, otherwise repeats the processing by the preprocessing means and the approximate solution update means,
A numerical calculation apparatus comprising:   5. The apparatus according to claim 4, further comprising parameter determining means for determining the scaling parameters γ and φ by estimating a parameter that minimizes an evaluation function defined by the coefficient matrix A and the preprocessing matrix M. The numerical calculation device described in 1. The evaluation function is a vector whose component is a sum for each column of a difference matrix between the coefficient matrix A and the preprocessing matrix M.

The numerical calculation apparatus according to claim 5, wherein the numerical calculation apparatus is an Euclidean norm.

Images