JP2010122850A - Device and method for calculating matrix equation - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a device and method for calculating a matrix equation, which enable high-speed calculation of a large sparse matrix equation appearing as an equation to be solved finally in numerical simulation, such as electromagnetic field analysis and which can be applied to both cases that the large sparse matrix of the large sparse matrix equation is a symmetric matrix and an asymmetric matrix. <P>SOLUTION: The matrix equation calculation device 1, which solves the large matrix equation by an iterative method, includes: a memory which includes a predetermined number of memories for storing only non-zero components for each row among the components of the large sparse matrix of the large sparse matrix equation; and one or more operation sections which perform at least part of operations of the iterative method with a data flow format. The device is configured such that data required for the operations of the iterative method for each line may be loaded to the operation section from the memory at once. <P>COPYRIGHT: (C)2010,JPO&amp;INPIT

Description

  The present invention relates to a matrix equation calculation apparatus and a matrix equation calculation method, for example, for calculating a large (large-scale) sparse matrix equation that appears as an equation to be finally solved in a numerical simulation such as electromagnetic field analysis, fluid analysis, and structural analysis. It is suitable for application.

  In many cases, scientific and technical calculations for solving a discretized partial differential equation in electromagnetic field analysis, fluid analysis, structural analysis, etc. ultimately result in a problem of solving a matrix equation with respect to unknowns to be obtained. In these scientific and technical calculations, in many cases, the matrix of the matrix equation to be solved is a sparse matrix whose order is several thousand to several tens of millions, and whose components (elements) are almost zero, That is, a large sparse matrix. Actually, in the electromagnetic field analysis, for example, as shown in FIG. 10, the matrix equation to be solved is different for each type of electromagnetic field such as electrostatic field / static magnetic field, quasi-static field, high frequency field, etc. There are numerical methods such as the difference method (FDM), the finite element method (FEM), and the boundary element method (BEM), and there are a lot of calculation methods by combining these methods. Most of the calculation methods result in the problem of solving large sparse matrix equations.

  In these science and technology calculations, the larger the calculation scale, the larger the ratio of the calculation time of the matrix equation to the total calculation time, often 95% or more. The demand for scientific and technological computations has become stronger and the problems to be handled have been increased in scale due to the recent improvement in performance of supercomputers and the enhancement of performance at the personal computer (PC) level. In particular, high performance computing (HPC) technology has become a very big issue in industrial applications where product development in a short period of time is required. Against this background, it is no exaggeration to say that speeding up the calculation of matrix equations, which account for the majority of scientific and technical calculations, is one of the biggest challenges of scientific and technical calculations overall.

  In general, there are two methods for solving matrix equations. In principle, direct solutions such as the LU decomposition method, in which unknowns are deleted one by one and the solution is obtained, and initial values are set and a true solution is obtained. This is roughly divided into an iterative solution method in which the iterative approximate solution is updated and the convergence value is obtained. Among them, the iterative solution method is a method of replacing linear simultaneous equations with a problem of obtaining a minimum value of an equivalent N-dimensional quadratic form (N is the order of a matrix), such as a conjugate gradient method (CG method), and such a problem. The matrix is classified into a simple iterative method in which the matrix is separated into a diagonal component and a non-diagonal component, and the approximate solution is updated only from the diagonal component. Furthermore, the method of replacing the linear simultaneous equations with the problem of obtaining an equivalent N-dimensional quadratic local minimum value is based on the CG method that applies a symmetric matrix to the biconjugate gradient method that extends it to an asymmetric matrix ( BiCG method), pre-processing the matrix before performing the iterative calculation, or configuring the optimal search direction vector in each iteration process to stabilize or speed up the iterative calculation, ICCG method, CGS method, Many derivative methods such as the BiCG-Stab method and the BiCG-Stab2 method have been proposed.

  Conventionally, it is considered that the above iterative method is effective for solving large sparse matrix equations among these methods. Large sparse matrices appearing in various numerical simulations such as electromagnetic field analysis, fluid analysis, and structural analysis. Often used to calculate equations.

  In an iterative solution of a matrix equation such as the CG method, an approximate solution of a linear (primary) simultaneous equation in a matrix form to be solved is obtained by the following procedure. First, with N being the number of unknowns or the order of a matrix, the linear simultaneous equations to be solved are replaced with a problem of obtaining an equivalent N-dimensional quadratic minimum value. Then, an appropriate N-dimensional solution vector is set as an initial value. The initial value (k = 0) is used as a starting point, and the following approximate solution (k = 1) is determined according to each iterative method such as the CG method. Further, the next approximate solution (k = 2) is obtained from this approximate solution (k = 1), and the procedure for obtaining the (k + 1) th approximate solution from this kth approximate solution is repeated. Then, this iterative calculation is performed until the (k + 1) -th approximate solution is within the preset accuracy, and when it is within this accuracy, that is, when the solution has converged, the above iterative calculation is terminated and this converged Let the value be the final approximate solution as a numerical calculation.

  As a technique for speeding up the iterative solution of the matrix equation, running on a parallel processing and high-speed CPU such as a PC cluster or a larger supercomputer is most commonly used. It is not a general environment that anyone can use easily because it can be used only in relatively limited research institutes, or it can be used jointly.

  On the other hand, against the backdrop of the recent emergence of rewritable LSIs such as FPGA and CPLD and their easy design / use environment, as well as the enhancement of the printed circuit board design / development environment, HPC technology options for scientific calculation As one of them, a dedicated computer that can be realized at low cost and in a familiar PC environment by narrowing down the target and configuring hardware with an architecture specialized for the calculation method is newly attracting attention. .

  At present, a dedicated computer for performing a linear matrix calculation (hereinafter also referred to as a matrix calculation dedicated computer) is targeted at a dense matrix, and a complicated LU decomposition process is performed in order to perform the LU decomposition method at high speed. A type that performs partial simple processing such as forward substitution processing with dedicated hardware, targeting a sparse matrix, and iterative processing and convergence determination processing of a simple iteration method with a relatively simple hardware configuration And a type in which the inner product operation of a vector or the multiplication of a matrix and a vector is implemented in hardware among the LU decomposition method or CG method for a sparse matrix (non- (See Patent Documents 1 to 13.)

Hideki Asai, Mitsuo Asai, Mamoru Tanaka, Dedicated Parallel Computer for LU Decomposition of Large Sparse Matrix, Science, Vol.J69-D, No.7, pp.1044-1053 (1986) Noboru Tanabe, Yasutaka Dohi, Construction of sparse matrix computer with associative switch, theory of theory, Vol.J70-D, No.12, pp.2393-2401 (1987) Noboru Tanabe, Kenichi Ishizaka, Hideki Murakoshi, Shoji Izumiya, Yasutaka Dohi, Construction of Parallel Pipeline Sparse Matrix Dedicated Computer RAMP, Science Theory, Vol.J71-D, No.10, pp.1939-1948 (1988) Yasushi Kiyoki, Toshiyuki Fukushige, Masato Taiji, Shinichiro Makino, Masaki Ogawa, Shunichi Kaizaki, Development of a dedicated matrix computer GENERAL-1, Computer Architecture, 111-9, pp.65-72 (1995) Masahiro Taiji, Minoru Fujita, Research on the development of quasi-parallel CPUs suitable for scientific and engineering computation, NEDO 2000 Proposal Public Proposal Report Meeting, Proposal Open: 99S29-019-1 (2000) Yosuke Ohno, Shunichi Kaizaki, Masahiro Taiji, Akihiko Konaseya, Effective Performance of Matrix Computation Engine (MACE) University of Hasegawa, Development of a dedicated computer for matrix operation for circuit simulation, Technical Report of IEICE, CA2002-116, pp.51-56 (2003) Yoshitaka Higuchi, Toshiaki Matsuoka, Hideki Kawaguchi, Research on the development of Laplace / Poisson equation finite difference solver dedicated computer, 2005 Japan Association of Electrical and Information Engineering Hokkaido Branch (2005), No.121. GRMorris, VKPrasanna, RDAnderson, A Hybrid Approach for Mapping Conjugate Gradient onto an FPGA-Augmented Reconfigurable Supercomputer, Proc. Of the 14th Annual IEEE Symposium on Field-Programmable Custom Computing Machines (FCCM'06), Napa, California (2006 ) R. Strzodka, D. Goeddeke, Pipelined Mixed Precision Algorithms on FPGAs for Fast and Accurate PDE Solvers from Low Precision Components, Extended technical report of IEEE Proc. On Field-Programmable Custom Computing Machines (FCCM'06), Napa, California (2006 ) O. Maslennikow, V. Lepekha, A. Sergyienko, FPGA Implementation of the Conjugate Gradient Method, LNCS 3911, pp.526-533 (2006) Tahara, Ippei, Kawaguchi, Hideki, Matsuoka, Shunsuke, Research on the development of computers dedicated to linear computation and conjugate gradient method, 2007 Hokkaido Association of Electrical and Information Society (2007), p.147. D. DuBois, A. DuBois, Sparse Matrix-Vector Multiplication and Conjugate Gradient Algorithm on Hybrid Computing Platforms, LALP-07-041 (2007)

  In a dedicated computer for matrix equation calculation, it is difficult to implement exactly the same processing as that performed by a software program-based computer. For this reason, conventionally, it is suitable for hardware processing in the overall processing and appropriately selects processing that can be effectively speeded up by hardware processing and executes them on a dedicated computer, while interlocking with the host computer A method for speeding up the processing as a whole has been adopted.

  Non-Patent Document 4 proposes a dedicated computer for Gaussian elimination with a dense matrix as a target. In particular, a method for hardware implementation of the forward erasure inner product operation, which accounts for most of the computation amount in the Gaussian elimination method, is proposed, and a computer system configuration is proposed. However, due to memory capacity limitations, the target is a matrix of 1000 dimensions or less, and further studies are necessary to perform large-scale calculations.

In Non-Patent Document 5, for the LU decomposition method dedicated computer that also targets a dense matrix, the calculation amount is N ² , N ^1.5 and the memory access for the data amount (matrix dimension) N in the LU decomposition method calculation. Focusing on the property that is difficult to become a bottleneck, we are investigating the method of a dedicated computer capable of high-speed processing in parallel, specializing in the inner product operation appearing in the LU decomposition method, and proposing the system configuration of the computer.

  Further, in Non-Patent Document 6, with regard to a dedicated computer for LU decomposition method, which also targets a dense matrix, only the four arithmetic operations of matrix components for LU decomposition of a matrix among the operations of LU decomposition method are hardwareized. We consider a method of processing by selection and forward substitution / backward substitution by a host computer, and propose a system configuration of the computer.

  At this time, in the attempt to speed up the matrix calculation by the dedicated computer for the above-described dense matrix, a large amount of memory is inevitably required, and a large amount of hardware is required to realize higher speed calculation. Due to the scale, it is limited to handle a matrix size of about several thousand dimensions at most, and there is a problem that it is difficult to apply to larger-scale calculations that often appear in scientific and technical calculations.

  On the other hand, Non-Patent Documents 1 to 3 propose a dedicated computer for a sparse matrix LU decomposition method targeting a sparse matrix. Utilizing sparse matrix properties, an arithmetic unit having a local memory corresponding to a much smaller number of non-zero matrix components than the dimension N of the matrix is prepared, and addition / subtraction / multiplication / division for one row can be performed in parallel. A system configuration is studied, and a method for executing all the processes of the LU decomposition method is studied, and a system configuration of the computer is proposed.

  In Non-Patent Document 7, two general-purpose arithmetic operators are prepared for a dedicated LU decomposition method computer that also targets a sparse matrix, and pipeline processing is possible with two channels of memory access. A small hardware is mounted on an FPGA evaluation kit having an actual PCI interface function, but its operation is confirmed.

  However, in an attempt to make the above LU decomposition method hardware, LU decomposition of a complex matrix including pivot selection and subsequent multiple processes of forward and backward substitution are required. The CG method and the BiCG method There was a problem that it was necessary to perform calculations with a large amount of calculation compared to the iterative solution method.

  Non-Patent Document 8 proposes a dedicated computer specialized for the Laplace-Poisson equation. In this dedicated computer, the number of non-zero components per row of the large sparse matrix appearing in the Laplace-Poisson differential equation is at most 5, and the matrix components appearing in the same row are the addresses of adjacent grids. A simple iterative method specializing in this structure has been proposed. However, although it has been found that the entire system can be realized with a simple architecture, the simple iterative method still has poor convergence of the solution, and is not necessarily appropriate as an HPC implementation of the matrix solution.

  Non-Patent Document 9 proposes a dedicated computer for the CG method targeting a large sparse matrix. In particular, only the multiplication operation of the matrix and the search direction vector in the CG method is hardwareized, and other processing is executed by the host computer, and the calculation of the CG method is performed while the host computer and the dedicated computer are linked. We are investigating a method for speeding up and proposing a computer system configuration.

  In Non-Patent Document 10, we examined the trade-off between accuracy and hardware size / calculation speed when assuming implementation on FPGA for a dedicated computer of the CG method targeting a large sparse matrix, and used pipeline processing together Appropriate calculation methods are proposed.

  Non-Patent Document 11 examines an efficient calculation method for division that appears in the calculation of the CG method for a dedicated computer of the CG method targeting a large sparse matrix, and the calculation part is actually implemented in an FPGA and effective. Showing sex.

  Non-Patent Document 12 discusses hardware implementation of one iteration of the CG method for a dedicated computer for the CG method targeting a large sparse matrix, and its feasibility while actually being implemented on a small scale FPGA. Is shown.

  Non-Patent Document 13 examines the HPC of matrix equation calculation using a commercially available reconfigurable computer for a dedicated computer for the CG method targeting a large sparse matrix, and proposes a system specialized for the CG method. is doing.

  As described above, the conventional computer dedicated to matrix calculation is not necessarily suitable for a large sparse matrix that often appears in scientific and technical calculations, or it is a hardware of a CG method suitable for a large sparse matrix. Even in the case of hardware, since the calculation of the CG method includes complicated processing such as division, the calculation method is partially hardwareized and operated in conjunction with the host computer. None of these dedicated computers for matrix calculation are not necessarily a complete dedicated computer, and are not necessarily sufficient in terms of operation speed and large-scale calculation such as data exchange with the host computer. The CG method is basically intended for a symmetric matrix, but an asymmetric matrix often appears in scientific and technical calculations, so that there is a problem that the scope of application is limited in a dedicated computer specialized for a symmetric matrix.

  Therefore, the problem to be solved by the present invention is that high-speed calculation of large sparse matrix equations that appear as equations to be finally solved in various numerical simulations such as electromagnetic field analysis, fluid analysis, and structural analysis, In addition, the present invention provides a matrix equation calculation apparatus and a matrix equation calculation method that can be applied to both cases where the large sparse matrix of the large sparse matrix equation is a symmetric matrix and an asymmetric matrix.

In order to solve the above problem, the first invention is:
A matrix equation calculator that solves a large sparse matrix equation of interest by an iterative method,
A memory unit having a predetermined number of memories for storing only non-zero components for each column among the components of the large sparse matrix of the large sparse matrix equation;
One or a plurality of calculation units that execute at least a part of the calculation of the iterative solution method in a data flow format,
The memory unit is configured to load data necessary for the calculation of the iterative solution for each row to the calculation unit at a time.
This matrix equation calculation apparatus is typically connected to a host computer such as a personal computer and obtains necessary data from the host computer, but is not limited to this. You may give the function of a computer.

The second invention is
A matrix equation calculation method for solving a target large sparse matrix equation by an iterative method,
Storing only non-zero components among the components of the large sparse matrix of the large sparse matrix equation in a predetermined number of memories in the memory unit for each column;
Loading data necessary for iterative solution operation for each row from the memory unit to one or more arithmetic units at a time;
And a step of executing at least a part of the calculation of the iterative solution in a data flow format by the calculation unit.

  In the first and second inventions, the large sparse matrix means a dimension of several thousand to several million or more appearing in a matrix equation to be finally solved in numerical simulation such as electromagnetic field analysis, fluid analysis, and structural analysis. However, most of its components are zero, and a non-zero component means a matrix having several to a dozen or more components per row. The large sparse matrix may be a symmetric matrix or a non-target matrix. For the iterative solution method, for example, CG method, BiCG method, ICCG method, BiCG-Stab method, CGS method, BiCG-Stab2 method and the like can be used.

から未知数ベクトルｘを求める問題において、これをＣＧ法の計算方法

を用いて解く場合を考える。（２−２）式の反復部分の計算には、

なる行列とベクトルとの乗算がある。このベクトルｑのｉ番目の成分ｑ_i の計算には行列Ａのｉ行目の１行分の成分の値が必要であるが、行列の成分の値を１つのメモリに格納しては、必要な値を１度にこれらのメモリから取得することはできない。このため、第１および第２の発明では、メモリ部が、行列の非ゼロ成分のみを１列ごとに格納する複数のメモリを有する構成とし、１行ごとの計算に必要なデータをメモリ部から演算部に並列に一度にロードすることとしている。また、後述の複数の演算部による並列処理の場合は、これらの列ごとにそれぞれ設けられたメモリをさらに並列数分設け、それぞれの並列処理が分担する行の成分の値を格納する。 In the first and second inventions, a predetermined number (plural) of storing only non-zero components among the components of the large sparse matrix for each column (selected as necessary, for example, 8 or 16) This memory is used for the following reason. FIG. 1 shows an example of a large sparse matrix equation. As shown in FIG. 1, a large sparse matrix equation consisting of a matrix A and an inhomogeneous vector b

In the problem of obtaining the unknown vector x from the above, this is calculated by the CG method

Consider the case of solving using For calculating the repetitive part of the equation (2-2),

There is a multiplication of a matrix and a vector. The calculation of the i-th component q _i of this vector q requires the value of the component for the first row of the i-th row of the matrix A, but it is necessary to store the value of the matrix component in one memory. A single value cannot be obtained from these memories at a time. For this reason, in the first and second inventions, the memory unit has a plurality of memories for storing only the non-zero components of the matrix for each column, and data necessary for calculation for each row is stored from the memory unit. The calculation unit is loaded in parallel at a time. Further, in the case of parallel processing by a plurality of arithmetic units described later, memories provided for each of these columns are further provided for the number of parallels, and the values of the row components shared by the respective parallel processings are stored.

Preferably, the matrix equation calculation apparatus includes a plurality of arithmetic units of the same hardware circuit. Then, the number of rows N of the matrix is divided by the number of these arithmetic units, and parallel processing is performed while the respective arithmetic units operate in conjunction with each other. By doing so, the calculation of the matrix equation can be speeded up by the parallel processing in a scalable manner by the number of operation units.
Preferably, iterative solution operations are performed in data flow form as much as possible.

  Further, preferably, the CG method dedicated computer for a symmetric matrix has an even number and a plurality of arithmetic units of the same hardware circuit. Then, by linking two of these arithmetic units in pairs, it is made to function as a BiCG method dedicated computer for one asymmetric matrix. Alternatively, the number of rows of the matrix may be divided by the number of pairs of the calculation units, and the parallel processing may be performed while each calculation unit operates in conjunction. In this way, one dedicated computer can be applied to both the symmetric matrix and the asymmetric matrix, and in the case of the symmetric matrix mode, it is possible to realize a computing device that operates at twice the speed of the asymmetric matrix mode. .

According to the present invention, a large-capacity memory can be mounted in a matrix equation calculation apparatus as long as it is allowed by hardware without software restrictions such as an OS, so that a large-scale calculation can be executed. In addition, the throughput is one time in the CG method or BiCG method, for example, by the data flow architecture circuit configuration of the iterative solution such as the CG method in the arithmetic unit and the memory configuration for each non-zero component of one column of the matrix in the memory unit. Can be processed with 3N + 2 clocks (N is the order of the matrix) and with the BiCG-Stab method or CGS method with 5N clocks, and the high-speed operation can be realized with the operation time minimized to the limit of the calculation method. Can do.
In addition, by operating in parallel with a configuration having a plurality of arithmetic units, it is possible to achieve high speed with scalable parallel processing.
Furthermore, since the BiCG method computing unit can be realized by linking two CG method computing units, efficient calculation suitable for the application can be performed.
In addition, a plurality of iterative solution calculation units are prepared, the same matrix calculation is performed simultaneously by these calculation units, and the host computer uploads the solution of the calculation unit that completed the earliest calculation, which automatically becomes a problem. The calculation of the matrix equation can be performed by selecting the optimal iterative solution method.

Embodiments of the present invention will be described below with reference to the drawings.
First, a matrix equation calculation apparatus according to a first embodiment of the present invention will be described.
This matrix equation calculation device is constituted by a computer dedicated to the CG method. This matrix equation calculation device calculates the matrix equation of equation (1) by the CG method shown in equations (2-1), (2-2), and (3).

FIG. 2 shows the overall configuration of the matrix equation calculation apparatus 1, and FIG. 3 shows the details of the configuration of each part of the matrix equation calculation apparatus 1. As shown in FIGS. 2 and 3, the matrix equation calculation apparatus 1 includes a memory module 11 that stores the values of all components of the large sparse matrix A, and an iterative calculation for one time of the CG method expressed by equation (2-2). the arithmetic module 12 to run with 3N clock data flow format, between the memory module 11 and the calculation module 12, (3) a large sparse matrix a and vector p _k not combining circuits clock delay multiplication of The matrix-vector multiplication circuit 13 that provides the result value Ap _k to the arithmetic module 12 and the (k + 1) th iterative solution obtained by the arithmetic module 12 have converged and become the final solution. And a convergence determination unit 14 for determining whether or not. The master controller 15 controls operations of a series of data flows in the memory module 11, the matrix-vector multiplication circuit 13, the arithmetic module 12, and the convergence determination unit 14.

  The memory module 11 includes m memories M1-1 to M1-m that store, for each column, the components of a matrix A ′ having a size of N × m obtained by compressing the large sparse matrix A into only non-zero components. Specifically, the large sparse matrix A is compressed as follows. For example, a matrix of N × m size is prepared, where m is the maximum of all N rows out of the number of non-zero components for each row, and the non-zero components are stored left-justified for each row. At the same time, another memory having exactly the same size as these memories M1-1 to M1-m is prepared (memory M2-1 to M2-m). In these memories M2-1 to M2-m, indexes indicating in which column of the original large sparse matrix A the non-zero components of the corresponding memories M1-1 to M1-m are stored. As a result, normally, in a large sparse matrix, N is several tens to several tens of millions, and m is at most a dozen, so that the memory capacity can be saved significantly without losing information of the large sparse matrix A. By configuring the memory module 11 in this way, the memory capacity can be greatly saved, and all the non-zero components of the matrix for one row can be accessed in parallel at a time. The number of accesses to the memories M1-1 to M1-m in which the matrix components that greatly affect the memory are stored can be reduced to a minimum of one.

The memory module 11 further stores memories M3-0 to M3-m for storing the values of the same vector p _k individually for each m in order to temporarily store values during execution of the CG method iterative calculation. memory M3-0 for storing the updated value p _{k + 1} of p _k ', the memory M4 for storing the vector x _k, the memory M4 for storing the updated value x _{k + 1} vector x _k', memory M5 for storing the vector r _k, memory M5 for storing the updated value r _{k + 1} vectors r _k ', the memory M6 for storing the vector q _k, 2 squared norm of the vector r _k ‖ r _k ‖ ² = (r _k, r _k) the memory M7 for storing the square ‖R _k ‖ ² update value of the norm ‖r _{k + 1} ‖ _{^{2 = (r k + 1,}} r k + 1 ), The order N of the matrix, the maximum number m of non-zero components per row, the non-homogeneous term vector b Memory for storing iteration convergence conditions such as the square of the norm of ‖b‖ ² = (b, b), the convergence criterion value ε of the relative residual, the maximum number of iterations Km, and the residual Rs of each iteration step It has M8.

  Memory M1-1 to M1-m, M2-1 to M2-m, M3-0 to M3-m, M3-0 ′, M4, M4 ′, M5, M5 ′, M6, M7, M7 ′, M8 Each of which may be an independent memory, or one or two or more of them may be obtained by dividing a memory area of one independent memory into a predetermined number. The memory areas are memories M1-1 to M1-m, M2-1 to M2-m, M3-0 to M3-m, M3-0 ′, M4, M4 ′, M5, M5 ′, M6, M7, M7 ′, It may be used for any of M8. These memories M1-1 to M1-m, M2-1 to M2-m, M3-0 to M3-m, M3-0 ′, M4, M4 ′, M5, M5 ′, M6, M7, M7 ′, M8 For example, SRAM or DRAM can be used.

The matrix-vector multiplication circuit 13 refers to the components of the large sparse matrix A divided and stored in the memories M1-1 to M1-m for each column as described above, while referring to the memories M2-1 to M2-m. simultaneously accessed in parallel to all of the components of one row every, also simultaneous access to the address of the required values for the calculation of each equation (3) from the memory M3-1~M3-m for storing the vector p _k, The calculation of the expression (3) is executed in the data flow format, the vector q _k is calculated for each component in one clock, and the result is delivered to the arithmetic module 12.

Calculation module 12, a matrix - vector q _k from vector multiply circuit 13, a vector p _k of memory M3-0, vector x _k of the memory M4, the vector r _k of the memory M5, the vector q _k and the memory M7 of the memory M6 (r _k, r _k) using, for performing the calculation of (2-2) formula (i) vectors p _k and vector q _k in (= Ap _k) the inner product of the (p _k, q _k) Arithmetic circuit 31, (2-2) arithmetic circuit 32 for executing the calculation of division (r _k , r _k ) / (p _k , q _k ) in formula (i), formula (2-2) ( (ii) an arithmetic circuit 33 for executing the calculation of the vector x _{k + 1} , (2-2) an arithmetic circuit 34 for executing the calculation of the vector r _{k + 1} of the equation (iii), (2-2) operation circuit 35 for performing the calculation of the square of the vector r _{k + 1} in formula _{(v) (r k + 1} , r k + 1), (2 Vector of 2) (iv) the division in _{(r k + 1, r k} + 1) / (r k, the arithmetic circuit 36 for performing the calculation of r _k), (2-2) Formula (v) An arithmetic circuit 37 for executing the calculation of p _{k + 1} is included.

  As can be seen from the equation (2-2), the calculation of the arithmetic circuits 31, 33, 34, 35, and 37 can be executed in N clocks, and the calculation of the arithmetic circuits 32 and 36 can be executed in one clock. 12 is configured so that these calculations can be executed in parallel as much as possible. Specifically, after the calculation of the N + 1 clocks of the arithmetic circuits 31 and 32 is completed, the calculation of the arithmetic circuits 33, 34, and 35 are collectively executed with N clocks using the data flow property of the arithmetic circuits 34 and 35. Therefore, the calculation of the N + 1 clocks of the remaining arithmetic circuits 36 and 37 can be combined and executed in a total of 3N + 2 clocks.

The convergence determination unit 14 calculates a relative error from the k + 1-th residual vector r _{k + 1} obtained from the matrix-vector multiplication circuit 13 that executes the calculation of the CG method and the iterative calculation for one time in the calculation module 12. Then, it is determined whether or not the iterative calculation is converged within the preset convergence criterion value ε of the relative residual in the memory M8.

If the calculation has converged, the master controller 15 ends the calculation with x _{k + 1} as the final approximate solution, and notifies the host computer (not shown) that the convergence value has been obtained normally. If the calculation does not converge, if the preset maximum number of iterations Km in the memory M8 is exceeded, the host computer is notified that the iterations have not converged, and the maximum number of iterations Km is not exceeded. Sends the vectors p _{k + 1} , x _{k + 1} , r _{k + 1} and (r _{k + 1} , r _{k + 1} ) obtained in the iteration calculation to the next iteration calculation, and again, the matrix-vector multiplication circuit 13 to instruct the execution of the iterative calculation in the arithmetic module 12.

FIG. 4 shows how the matrix equation calculation apparatus 1 is used. The matrix equation calculation apparatus 1 is connected to a host computer 42 by a PCI interface 41. The values necessary for executing the calculation of the CG method are a large sparse matrix A that is compressed and stored, vectors p ₀ (= r ₀ ) and x ₀ that are initial values for iterative calculations, and (b, b ), Convergence criterion ε, and maximum number of iterations Km. These values are stored in the host computer 42, downloaded to the matrix equation calculation apparatus 1 before the calculation start instruction, and stored in the memories M3, M4, M5, and M8 of the memory module 11. The host computer 42 downloads these values to the matrix equation calculation apparatus 1 and then instructs the matrix equation calculation apparatus 1 to start calculation. The matrix equation calculation apparatus 1 repeats the iterative calculation and convergence determination of equation (2-2), and approximates the vector x _{k + 1} in the memory M4 ′ if the solution falls below the convergence criterion ε within the maximum number of iterations Km. If it is delivered to the host computer 42 as a solution and does not converge even if the maximum number of iterations Km is calculated, the calculation is not converged, and this is notified to the host computer 42.

5A to 5E show implementation examples of the matrix equation calculation apparatus 1. As shown in FIGS. 5A and 5B, in this example, on one small printed circuit board 51, a memory 52 for storing compression matrix components for two columns of the memory module 11, and a search corresponding to the matrix of the two columns A matrix-vector multiplication circuit 13 that performs multiplication with the component of the vector p _k is mounted. As shown in FIGS. 5C and 5D, a total of 16 small printed circuit boards 51 are connected to the sub printed circuit board 53. As shown in FIG. 5E, the arithmetic module 12, the convergence determination unit 14, and the master controller 15 are mounted on the main printed circuit board 54, and the sub printed circuit board 53 is connected thereto. Thus, the matrix equation calculation apparatus 1 as a computer dedicated to the CG method is configured. Most of these mounted components can be manufactured by, for example, FPGA or ASIC. With this configuration, even when the matrix equation calculation device 1 is configured for a predetermined number m of columns, if a sufficient socket is prepared in advance, the small printed circuit board 51 is additionally connected later. It is possible to have flexibility for changing the matrix size.

According to the first embodiment, the following advantages can be obtained. That is, the matrix equation calculation apparatus 1 is configured by a computer dedicated to the CG method, and the memory module 11 and the arithmetic module 12 are separated from each other, so that the memory module 11 is appropriately configured according to the large sparse matrix used for the numerical simulation. The arithmetic module 13 having a relatively complicated configuration can be replaced without changing, and can flexibly cope with the change of a large sparse matrix, and has high flexibility. Further, the memory module 11 is provided with m memories M1-1 to M1-m for each column of the matrix, and is necessary for calculating one row of the equation (3) from these memories M1-1 to M1-m. Data is loaded into the matrix-vector multiplication circuit 13 in one clock, and the calculation for one row of the equation (3) is executed in parallel, thereby minimizing the number of accesses to the memory module 11. Accordingly, high-speed operation can be achieved, and the matrix equation can be calculated at high speed in a short time. In addition, the matrix equation calculation apparatus 1 can be easily implemented and scaled up. As described above, it is possible to realize a high-performance matrix equation calculation apparatus 1 with extremely high practicality.
The matrix equation calculation apparatus 1 is suitable for application to solving a large sparse matrix equation that appears as an equation to be finally solved in various numerical simulations such as electromagnetic field analysis, fluid analysis, and structural analysis.

Next explained is a matrix equation calculation apparatus according to the second embodiment of the invention.
In the second embodiment, in addition to minimizing the number of accesses to the memory module 11, a plurality of arithmetic modules 12 are prepared to further increase the calculation speed of 3N + 2 clocks in the arithmetic module 12, and 3N + 2 A method is used in which the calculation for the clock is executed in parallel by the plurality of arithmetic modules 12. The calculation of 3N + 2 clocks in the arithmetic module 12 may be all processed in parallel originally, and by performing this part of the calculation with a plurality of arithmetic modules 12, the number of arithmetic modules 12 is larger than when one arithmetic module 12 is used. The speed can be increased by that amount. Other than the above are the same as in the first embodiment.
As shown in FIG. 6, the matrix equation calculation apparatus 1 can be easily mounted by mounting a plurality of sub print boards 53 on a main print board 54 similar to that shown in FIG. 5E.

  According to the second embodiment, in addition to the same advantages as the first embodiment, the operation speed can be further improved, and the matrix equation can be calculated at higher speed. The advantage that can be obtained.

Next explained is a matrix equation calculation apparatus according to the third embodiment of the invention.
In the third embodiment, in addition to speeding up by parallel processing of a plurality of arithmetic modules 3, the CG method is limited to application only when the large sparse matrix is a symmetric matrix. In order to provide a function that can be applied even in the case of an asymmetric matrix, a method of operating two arithmetic modules 12 in an interlocked manner and operating as a BiCG method machine that can also be applied to an asymmetric matrix is used. The BiCG method is originally obtained by adding processing of the CG method to the transposed matrix to the CG method. By causing the other computing module 12 to process the transposed matrix, the two computing modules 12 are linked and operated in pairs. The behavior of the law can be realized. Other than the above are the same as in the first embodiment.
As shown in FIG. 6, the matrix equation calculation apparatus 1 can be easily mounted by mounting a plurality of sub print boards 53 on a main print board 54 similar to that shown in FIG. 5E.

となる。この方程式を半径方向、軸方向ともにΔｌのサイズの一様なグリッドで分割した上で軸対称性を考慮して差分法に基づき差分式で表すと、半径方向にi 番目、軸方向にj
番目のグリッドでは、

となる。ただし、μ_i,j 、Ａ_i,j 、Ｊ_i,j はそれぞれ、半径方向にi 番目、軸方向にj 番目のグリッドでの透磁率μ、ベクトルポテンシャルＡのθ成分、電流密度ベクトルＪのθ成分である。このとき、未知数Ａ_i,j の数は、全てのグリッド数の５０×１００＝５０００であるので、差分法により数値解析する際に現れる大型疎行列は５０００×５０００のサイズとなる。この大型疎行列計算を、非対称行列用のＢｉＣＧ法に基づいて構成した行列方程式計算装置１により実行する場合のＶＨＤＬ論理回路シミュレーションにより行った。その解から得られた静磁場分布の結果を図９に示す。実際の計算では、（５）式からわかるように、１行に現れる未知数Ａ_i,j の係数は高々５個であるので、この疎行列性を利用し、行列は５０００×５のサイズに圧縮したものを用いている。図９に示す結果は、同じ数値モデルに基づいてＣ言語ソフトウェアシミュレーションを行った結果と一致しており、この行列方程式計算装置１の妥当性が保証される。このＶＨＤＬ論理回路シミュレーションでの動作を参考に、この行列方程式計算装置１がアクセス速度２６６ＭＨｚのメモリと同程度の速度で動作した場合の性能を見積ると、２．６ＧＨｚのＣＰＵを有するパソコン上でＣ言語を用いて計算した時の約５０〜６０倍の性能に匹敵する。 Next, the result of logical simulation of large sparse matrix calculation by the matrix equation calculation apparatus 1 using VHDL (VHSIC Hardware Description Language) which is a hardware description language will be described.
FIG. 7 shows a model of the magnetic field portion of the linear motor used as an example of simulation, and FIG. 8 shows its numerical model. This model relates to a phenomenon in which a static magnetic field excited by three exciting coils 61 to 63 is shielded from being exposed to the outside by a magnetic shielding material 64 made of iron (μ _r = 5000). 8 has a two-dimensional 50 × 100 grid size in consideration of axial symmetry. The governing equation describing this magnetic field is that for a given current density vector J and permeability μ, the vector potential A is an unknown,

It becomes. When this equation is divided by a uniform grid of size Δl in both radial and axial directions, and expressed by a differential equation based on the difference method in consideration of axial symmetry, it is i-th in the radial direction and j in the axial direction.
In the second grid,

It becomes. However, μ _{i, j} , A _{i, j} and J _{i, j} are the magnetic permeability μ, the θ component of the vector potential A, the current density vector J, and the i-th grid in the radial direction and the j-th grid in the axial direction, respectively. The θ component. At this time, since the number of unknowns A _{i, j} is 50 × 100 = 5000, which is the number of all grids, the large sparse matrix that appears when performing numerical analysis by the difference method has a size of 5000 × 5000. This large sparse matrix calculation was performed by VHDL logic circuit simulation in the case of executing by the matrix equation calculation apparatus 1 configured based on the BiCG method for an asymmetric matrix. The result of the static magnetic field distribution obtained from the solution is shown in FIG. In the actual calculation, as can be seen from the equation (5), since the coefficient of the unknown A _{i, j} appearing in one row is at most 5, this sparse matrix property is used, and the matrix is compressed to a size of 5000 × 5. We use what we did. The result shown in FIG. 9 matches the result of the C language software simulation based on the same numerical model, and the validity of the matrix equation calculation apparatus 1 is guaranteed. Referring to the operation in the VHDL logic circuit simulation, the performance when the matrix equation calculation apparatus 1 operates at the same speed as a memory with an access speed of 266 MHz is estimated on a personal computer having a 2.6 GHz CPU. It is comparable to about 50-60 times the performance when calculated using language.

Here, for the sake of simplicity, an axisymmetric two-dimensional model is used, but basically the same applies to a three-dimensional model. That is, for example, in the same three-dimensional model of FIG. 7, the grid size is 100 × 100 × 100, so the large sparse matrix that finally appears is 1000000 × 1000000, and the three-dimensional difference equation corresponding to equation (4) Then, since there are 16 coefficients, the result is finally to solve a matrix compressed to a size of 1000000 × 16. For this reason, in this matrix equation calculation apparatus 1, in the case of the two-dimensional axisymmetric case, N = 5000 and m = 5, but in the three-dimensional case, only N = 1000000 and m = 16. The operation itself is exactly the same.
According to the third embodiment, in addition to the advantages similar to those of the first embodiment, the matrix equation calculation apparatus 1 can be applied even when the large sparse matrix A is an asymmetric matrix. Can be obtained.

Although the embodiment of the present invention has been specifically described above, the present invention is not limited to the above-described embodiment, and various modifications based on the technical idea of the present invention are possible.
For example, the numerical values, circuit configurations, arrangements, and the like given in the above-described embodiments are merely examples, and different numerical values, circuit configurations, arrangements, and the like may be used as necessary. The matrix equation is not solved by the CG method or BiCG method, but other iterative methods such as the ICCG method, BiCG-Stab method, CGS method, BiCG-Stab2 method are adopted, and the calculation module 12 also uses these matrix equations. You may comprise for the solution of these.

It is a basic diagram which shows one example of a large sparse matrix equation. 1 is a schematic diagram showing an overall configuration of a matrix equation calculation apparatus according to a first embodiment of the present invention. It is a basic diagram which shows the detail of a structure of each part of the matrix equation calculation apparatus by 1st Embodiment of this invention. It is a basic diagram which shows the state which connected the matrix equation calculation apparatus by 1st Embodiment of this invention with the host computer. It is a basic diagram which shows the example of mounting of the matrix equation calculation apparatus by 1st Embodiment of this invention. It is a basic diagram which shows the example of mounting of the matrix equation calculation apparatus by 2nd Embodiment of this invention. It is a basic diagram which shows the numerical model used for the numerical simulation performed using the matrix equation calculation apparatus by the 3rd Embodiment of this invention. It is a basic diagram which shows the numerical model used for the numerical simulation performed using the matrix equation calculation apparatus by the 3rd Embodiment of this invention. It is a basic diagram which shows the result of the numerical simulation performed using the matrix equation calculation apparatus by 3rd Embodiment of this invention. Taking electromagnetic field analysis as an example, this is a schematic diagram for explaining that many scientific and technical calculations ultimately result in calculations of large sparse matrix equations.

Explanation of symbols

  DESCRIPTION OF SYMBOLS 1 ... Matrix equation calculation apparatus, 11 ... Memory module, 12 ... Arithmetic module, 13 ... Matrix-vector multiplication circuit, 14 ... Convergence determination part, 15 ... Master controller, 42 ... Host computer, 51 ... Printed circuit board, 53 ... Subprint Substrate, 54 ... main printed circuit board, M1-1 to M1-m, M2-1 to M2-m, M3-0 to M3-m, M3-0 ', M4, M4', M5, M5 ', M6, M7 , M7 ', M8 ... memory

Claims (8)

A matrix equation calculator that solves a large sparse matrix equation of interest by an iterative method,
A memory unit having a predetermined number of memories for storing only non-zero components for each column among the components of the large sparse matrix of the large sparse matrix equation;
One or a plurality of calculation units that execute at least a part of the calculation of the iterative solution method in a data flow format,
A matrix equation calculation apparatus configured to load data necessary for iterative solution operation for each row from the memory unit to the operation unit at a time.   The matrix equation calculation apparatus according to claim 1, comprising a plurality of the arithmetic units.   3. The matrix equation calculation apparatus according to claim 2, wherein the plurality of arithmetic units are operated in parallel.   4. The matrix equation calculation apparatus according to claim 3, wherein the iterative solution method is a conjugate gradient method or a biconjugate gradient method. A matrix equation calculation method for solving a target large sparse matrix equation by an iterative method,
Storing only non-zero components among the components of the large sparse matrix of the large sparse matrix equation in a predetermined number of memories in the memory unit for each column;
Loading data necessary for iterative solution operation for each row from the memory unit to one or more arithmetic units at a time;
And a step of executing at least a part of the calculation of the iterative solving method in a data flow format by the calculation unit.   6. The matrix equation calculation method according to claim 5, wherein a plurality of the arithmetic units are used.   7. The matrix equation calculation method according to claim 6, wherein the plurality of arithmetic units are operated in parallel.   8. The matrix equation calculation method according to claim 7, wherein the iterative solution is a conjugate gradient method or a biconjugate gradient method.

Images