JPS61175774A - Computer for analysis of simultaneous equations - Google Patents

Abstract

PURPOSE:To attain the LU analysis of a spurse matrix with high efficiency by storing only the non-zero elements into a local memory of each local unit after ignoring the greater part of the zero elements of an original matrix Y to convert the matrix Y into a band matrix and flowing the elements successively in parallel and every line to (q) pieces of processors corresponding to the band width. CONSTITUTION:This processor contains (q) pieces of local units. Then, the local units 1-3 are all connected with each other via a common data bus 10, a control line 11, a priority circuit 12, etc. All non-zero elements only including the fill-in are stored in each local memory unit {Aj} in the row direction of a matrix after ignoring all zero elements. The operations of the processor are controlled in a division mode (a), a multiplication mode (b) and a gauss erasion mode (c) respectively. In the division mode the unit {Aj} receives an access by a counter 13. Only the local bus corresponding to a pivot aii is enabled by an enabling circuit 15. Then the pivot data aii are transferred to all local processors via a common bus.

Description

Detailed Description of the Invention [Technical Field of the Invention] The present invention relates to a method of configuring a parallel computer for LU decomposition of a large-scale sparse matrix.

[Prior art and its problems] In parallel with the development of circuit network analysis algorithms based on the concept of parallel processing, many inventions have been made of dedicated computers for executing these algorithms on large-scale systems. , is making rapid progress.

The fundamental reason for the increasing need for dedicated computers is that the sequential processing of general-purpose computers requires too much analysis time for analysis of large-scale systems.

In recent years, the concept of systolic array processors has become popular as a dedicated computer for increasing parallelization and pipelining of processing associated with direct methods such as Gaussian elimination and LU decomposition.
g [1], and its applications have been frequently described [2-5].

Theoretically, array processors can be considered as special-purpose machines that can make the most effective use of their parallelism in data processing, but they cannot always be considered without the concept of infinite hardware. In order to efficiently analyze a large number of sparse matrices, there is a serious problem that data transfer becomes complicated. Even if you target 0, you need too many processor cells, for example.

The LU decomposition processor utilizes the parallelism inherent in Gaussian elimination processing, but the sparsity of the matrix cannot be ignored, and in order to handle a square matrix with nXn dimensions, Requires n2 processor cells. This can be said to be an unrealistic number for analysis of large-scale circuit networks where n is very large.

Therefore, in order to simplify such problems, many 7
The scope of application of the ray processor is limited to solving nodal equations in which the coefficient matrix has a band structure. This greatly simplifies data communication and allows parallel processing and pipeline processing to be performed very efficiently.

For a band matrix with a bandwidth of 2p+1.

If approximately 22 processor cells are prepared for H>)p, it becomes possible to use the systolic data method,
Therefore, n-dimensional simultaneous equations can be analyzed in O(n) time.

However, the structure of the coefficient matrix Y of the nodal equation depends not only on the topology when the circuit network is represented graphically, but also on the method of labeling the nodes.
It generally does not have a band structure.

In order to shorten analysis time and save memory capacity, there are two labeling methods: band method and backing method [6-6 packing method reduces fill-in during matrix calculation process. Labeling is a technique, and by reducing the number of fill-ins, the calculation of zero elements can be efficiently omitted, resulting in a reduction in analysis time. One of them is the optimization of the pivot order between 1g [7
] There is. In this method, the number of fill-ins is extremely small, but the positions of non-zero elements in the matrix are completely irregular.

On the other hand, in the band method, the labeling method is such that all non-zero elements of the coefficient matrix Y are contained within a certain bandwidth, so the number of fill-ins increases compared to the backing method, but the number of fill-ins increases. Since the locations of non-zero elements including ins are limited, data communication is simplified and analysis time is eventually shortened.

In a sparse matrix that follows the backing method, if non-zero elements appear in positions far from the diagonal elements of the matrix, if we dare to define a bandwidth of 2p+1, the result will be p/%-n, which will eventually form an array. The number of processor cells required is 22,112.

Therefore, if an array processor were applied to such a matrix, an impractical number of processor cells would be required, and too many cells would be idle for efficient analysis. Therefore, the elements of the coefficient matrix in large-scale simultaneous equations exist irregularly,
Moreover, a dedicated computer that can handle extremely strong sparsity is required.

[Object of the Invention] The present invention presents a method of configuring a parallel computer for LU decomposition of a large-scale sparse matrix. In this processor, nine local units corresponding to the number q of non-zero elements including fill-ins existing in each row of the coefficient matrix Y are operated in parallel. In the local memory in each local unit, ignoring most of the zero elements in the original matrix Y,
Memory capacity is saved by storing only non-zero elements. Furthermore, the data mismatch that occurs during the Gaussian erasure calculation process is eliminated by a data shift operation or by sequential writing and parallel reading of the memory. The results show that LU decomposition of sparse matrices according to the backing method, etc., which is practically difficult in systolic array processors, can be performed as efficiently as sparse matrices using the band method.

[Summary of the Invention] In order to achieve the above object, the present invention provides the element at2 (i=
1.2. ,,,,N;j=1,2. , , , , N), in each i-th row of the sparse matrix Y, the non-zero elements including Lyle in are left-justified in the row direction, and the value and column number of each element aLk of each column vector are as data.

The first local memory to be stored at address i is set up to have at most q the maximum number of non-zero elements including fill-ins that exist in each row of matrix Y, and In each j-th column of the triangular matrix part, at least the row number statement of each element a of each row vector created by top-filling non-zero elements including fill-in in the column direction is stored at address i as data. 2 local memories are installed in a number smaller than the number of the first local memories, and in the division mode, each element read from each of the first local memories by the address i and the column number are A dividing means detects the pivot element aLL which satisfies i=, and simultaneously divides each of the read elements aLk to obtain the pivot element ar<1' to obtain aYO/aLL, and in the multiplication mode, the address i Each row number sentence of each element aeJ read from each second local memory by
Read the contents aHa of the local memory, detect the element a41ε where β=i, and multiply each value aLk/att obtained by the dividing means by the element ali.
a multiplication means for calculating i'aik/a, ·;
Using a shift means to match column number j of element a2j read from the local memory of , each value "Jl(-aLk' al, Provided is a computer for simultaneous equation analysis, characterized in that it has a Gaussian elimination calculation means for subtracting . , wherein the first local memory is accessed by the contents of the second local memory, and the dividing means.

When the number p of processors that execute the multiplication means and the Gaussian elimination operation means is smaller than the maximum number q, the two first local memories are divided into r so that r×p exceeds 9.
Provided is a computer for analyzing simultaneous equations that is characterized by being hierarchically connected to levels.

[Action Co.] In this processor, in order to operate q local units in parallel corresponding to the number q of non-zero elements including fill-in existing in each row of the coefficient matrix Y, the local By storing only non-zero elements in memory, ignoring most of the zero elements in the one-element matrix Y.

It is based on the idea of converting the matrix Y into a band matrix and passing the elements line by line in parallel and sequentially to the q processors corresponding to the band width.

[Embodiments of the Invention] Examples of the present invention will be described below with reference to the drawings.

First, the LU decomposition method will be described.

Generally, the I'I n solution for the one-node equation Yv = j forv is
It is performed according to the following steps using LU decomposition, forward substitution, and backward substitution.

1) Y=LU for L, U2) L
ξ=j for ξ3) U v r,
=ξ forvnTo complete these processes, O(n3) multiplications and additions are required, but
Most of the time is spent in the LU decomposition process, so in order to solve a one-node equation in a short time, it is necessary to perform the LU decomposition efficiently.

In the i-th row pivot operation, the i-th row and column in the matrix Y are transformed according to the first-order equation (1).

In these processes, the elements ., j are updated to a kj) according to equation (2).

a −' = a −a , −a−7m” (2)
k, 1 kJ kL'
JLL That is, the fill-in occurs at the i-th pivot point j). In a banded matrix, all fill-ins will occur within the band.

Below, we will propose a method for configuring an LU decomposition processor for a matrix with a sparse structure, which is difficult to practically apply to a systolic array processor.

This processor is composed of q local units, and as shown in FIG. There is.

The interior of each local unit is shown in FIG.

The details of one processor will be explained below.

1) Local memory and data structure A) Local memory unit (A, 2) Each local memory unit (A21 j=1.2°3y ---
q) stores only all non-zero elements, including fill-ins, in the row direction of the matrix, with all zero elements ignored.

Therefore, the address of each local memory (A2>) is equal to the row number in the matrix Y of the data stored at that address, and the data format is (k, number), where is the original is the column number in the matrix,
The "number" is 1. For example, in a resistor network, it corresponds to the branch conductance.

B) Local memory unit (Bj) Each local memory unit (Bjl j=1.2゜3+, −-* q
) has 1 for the lower triangular matrix excluding the diagonal elements of the matrix Y.
Only all non-zero elements, including fill-ins, are stored in the column direction of the matrix, with all zero elements ignored.

Therefore, the address of each local memory (B, j) is
The only data format that is equal to the column number in matrix Y of the data stored at that address is i. Here, 1
is the row number in the original matrix of the data.

The data structure of the local memory is shown using the graph shown in FIG. 3 as an example.

The structure of the nodal conductance matrix for the graph structure shown in FIG. 3 is shown in FIG.

Here, the symbols "1" and "f" mean a non-zero element and fill-in, respectively.

For the matrix structure shown in Figure 4, the local memory (
The data structures in A,i) and (Bj) are.

The result will be as shown in Figures 5 and 6.

In other words, for non-zero elements due to fill-in, only the label number of that data is stored.

The numerical value follows matrix Y and is set to zero.

2) Operations by local processor Operations in this processor are controlled by three operation modes. The calculation modes are a) division mode (÷),
b) Multiplication mode (X), Q') Subtraction (Gaussian elimination) mode (G).

an accumulator A inside each local processor;
The content of C [ACc] is 1 expression (2
), updated as follows:

a) Acc″”Lj '” LLb) A
4-[A(c)''a,;C(c) A 4-a --[ACC]CCJ Below, the operation of the processor in each calculation mode will be explained.

a) Division Mode In FIG. 1, the symbol "i" means the pivot row (column) number for the pivot a cj. First, in the local memory (Aj), the address i for the pivot atE is the counter 1 that stores the pivot number i.
Accessed by 3. Data (k, numeric value) is set to register (R°) and (R, j co-). The column number k of a R:LJ is compared with the pivot number i, and the local to be assigned to the pivot aLi is The unit is detected by exclusive OR 14. As a result, only the local path corresponding to the pivot aLL is enabled by the enable circuit 15,
Pivot data ati is transferred to all local processors via a common bus.

The processor P of each local unit executes division of the above "LL" by the contents of R3i. The result is written to the local memory (Aj), and (R3j)
is set to

b) Multiply mode In the divide mode, address i of local memory (Aj) is accessed and at the same time address i of local memory (B
In zero multiplication mode, address i in j) is accessed by pivot number i and its data p is set in register (R,j).

Local memory (Aj) is accessed with the contents of register (R,j). Its content is the row number of the nonzero element of the lower triangular matrix in column i. Register (R, Δ1J=2
The contents of #3t ---) are sequentially shifted to the next higher register (R, .), and the data accessed with the contents of register RXl is set to registers R4 and R4A. . Then, as in the division mode, the column number data of R24 is compared with the pivot number i, and one column element is detected. The detected value is transferred to all other processors on a common bus and multiplied with the contents of R.

The calculation result is a in R33. Set.

C) Gaussian elimination mode Next, Gaussian elimination operation 1I is executed using the contents of R and the contents of R9j, but the pair of registers (R4J + R
y, 4) data [Rb2], [: Rg;
] may not match. In other words, in this processor, which is designed for analyzing large-scale circuit networks, data is stored in a form that ignores all zero elements in order to save memory capacity, so in the Gaussian elimination process, registers (R, , At the initial stage when data is set in i), the label number of the addition data (column number in the 7 matrices) and the label number of the augend data (column number in the 7 matrices) are often different. Therefore, in order to perform addition (subtraction) between elements at the same position in the original matrix, the contents of register R5a are shifted,
The contents of register Rqj are compared with the label number.

When all pairs of label numbers in ([R4・]) and ([R, j) match, that is, when the start signal becomes "1", addition (subtraction) is executed in Gaussian elimination mode. be done. The result of the operation is then written to the local memory (A,i) at the original address using the address from the contents of R+a.

3) Shift operations Here, the shift operations mentioned in 2) are explained and the influence of these shift operations on the calculation time is considered.

Before the start of the shift operation, the add data and the augend data in the Gaussian elimination operation are set in register (R, j) and register (R, 2), respectively. Here, the label numbers of the contents of register R source R and 5j are compared, and if they do not match, the contents of register R width are right shifted to register R6 layer 1. In principle, this operation is repeated sequentially.

During the right shift process, what is the maximum number of shifts required until the label numbers of all elements match?

It is clear that it depends on the number of non-zero elements, including fill-ins, contained within each matrix of the matrix Y. In the right shift operation, the contents of registers (Raj) and (R,j) have the following properties.

[Property 1] If the matrix Y is structurally symmetric, ([R smell] ) will include all the fill-ins in the i-th row that may occur during the multiplication process of the i-th row and the i-column, so the register The contents of ([R,,il) are the contents of the register ([R4,i
l). That is, ([R・] l ,)
= l, 2e 3+ ---t q) d Ho ([R・] IJ"1+ 213901. r q) It is a smell.

[Property 2] The number of shifts of the initial value [R5-, lr] of each register is T・F ''It 2+ 3+ , -, v
r≦q), then from C property 1, T : ii; T ≦01. ≦T1, where r is the number of non-zero elements included in the i-th row of the upper triangular matrix of the matrix Y.

[Property 3 [Property 2] From property 2, the contents of each register [R
The number of times that 5jl should be shifted is at most T1 times.

[Property 4] Let the numbers of non-zero elements in the i-th row and the i-th row in the matrix Y be N ZL and N 71, respectively.

T =lNZ・-NZflll+1 r　　　　　　　　　L It is.

The above properties can also be applied to band matrices.

That is, for the LU decomposition of a band matrix, since NZj = NZ1 in most cases, the maximum number of shifts is T7 = 1, and one shift can shift the label number in all register pairs. Alignment will take place. For dense matrices T7=0.

If the label numbers in register R (j = m) and register R4j (j = J m) match, register (R5
・l j=rn + 1 v m + 2 t −1,
Only the contents of t r ) are shifted. And then, Jister R4ff...

The numerical value "0" and a label such as R9 are forcibly set to .

By the way, there may be a case where the label numbers of two or more pairs of data in registers (R4j) and (Rrj) coincide at the same time, but in this case, the following property exists.

[Gt-51 When two or more pairs of label numbers match at the same time, those registers are always adjacent.

In this case, only the data in those registers contained in (R,j) needs to be right shifted.

In this processor, the result of detecting a match between the label numbers of the registers (R4j) and (Rrj) is used as a flag, latched into a flip-flop, and passed through the right priority circuit [8] to determine the match. This is accomplished by detecting the rightmost register among the registers and using its output to enable the registers to the right.

If the number of L branches connected to the node nL corresponding to pivot a, - is Inc (n.), then the number of non-zero elements NZL in the i-th row is equal to Inc (n2)+1.

That is, when the graph structure of one circuit network is a regular graph, the number of non-zero elements in each row is equal, and operations are performed most efficiently.

In this processor, the calculation time including shift operations is not determined by the structure of the matrix (i.e., the labeling of graph nodes), but by the number of nonzero elements in a row. Even in the case of a backing method that reduces the number of fill-ins,
The number of shifts is not increased, and therefore efficient LU decomposition is performed.

In principle, a right shift aligns the label numbers, but for general sparse matrices, a left shift may be required. In this case, as shown in FIG. 7, label numbers can be matched by a left shift operation using a left priority circuit.

The evaluation of the efficiency of a single processor will be described below.

The calculation efficiency E of the computer is E = O (time) X N (co-proce
ssor). Here, O(time) and N
(co-processor) is the time required for solution solving, the number of processor cells, etc., respectively.

For the LU decomposition of a sparse matrix with dimensions nXn,
The efficiencies Ea and Es of the array processor and this processor are respectively.

Eo+! It can be expressed as n×p”=np” R3=nqXq/2+α=nq”/2+a. Here, α is the time required for data communication, and is considered to be the time required for the shift operation. If a barrel shifter or the like is used, it may be possible to match the label numbers instantaneously without a shift operation, but if INZL-NZQI is not that large, α will not take a very large value.

In the band matrix, q:2p. Array processors cannot take advantage of the sparsity of matrices, so if labeling such as the backing method is used,
It may also be p /% Q, in which case q<<p, hence.

]1:, << EQ Therefore, one processor is considered to be effective regardless of the band method or backing method, as long as the number of fill-ins is suppressed to some extent and the labeling is used to suppress the spur input matrix. Furthermore, if the data in the register R55 is shifted to the pivot position and then the Gaussian operation is performed, the number of local processors can be reduced to about half, although the number of shifts increases. The efficiency in this case is E
5'.

In the following, using a concrete example, we compare the efficiency of a single processor and a 7-ray processor through simulation.The structure of the coefficient matrix Y of the circuit network shown in Reference [7] is shown in Figure 8. shown.

According to the optimization of the pivot order shown in Reference [7], the structure of the matrix Y becomes as shown in FIG.

Consider using this processor for the matrices shown in FIGS. 8 and 9. 808 for each local processor
6, and the number of clocks required for each calculation is set as follows.

Computation efficiency E, E for the solution of nodal equations having coefficient matrices with matrix structures shown in FIGS. 8 and 9
and E, J are shown in Table 2.

α 5 Also, find the number of nonzero elements in each row of the original matrix.

It has been confirmed that when simple optimal ordering, in which ordering is performed from the smallest number of rows, is applied to the matrix shown in FIG. 8, the computational efficiency is almost equal to that achieved by optimal ordering.

Next, another embodiment in which the shift operation is performed using memory will be explained, and further, an explanation will be given of an embodiment regarding the system configuration when the number of non-zero elements included in a row exceeds q. Do this using diagrams.

Local memory (AIMAl□t A13) e (
AzlpA, . . . A23) is filled with zero elements other than fill-in in the row direction of the matrix as shown in FIG. 5, as will be described later. The data format is (k, number), and k
is the column number, and the numerical value corresponds to the branch conductance.

The local memory Bj a (Bj S j = 1.2) is filled in the column direction for the lower triangular matrix excluding diagonal elements.
Taking the graph shown in Figure 3 as an example, which is stored with zero elements other than IN, but does not require numerical values, the data structure of each local memory will be as shown in Figures 5 and 6. . The operation of the processor will be explained below. Local memory Ajk, where i indicates a pivot, is accessed with address i, set to R21, a JE element is detected from the column number, and transferred to all processors to perform division. The calculation result is stored in the local memory Cjε(Cjjj=1*21
, , I") with the column number as the address, and is also copied to other local memory Cj. Local memory BJ is also accessed at address i, and the row number of the non-zero element in column i is set to RIJ. When local memory AJI is accessed with the contents of RlJ in multiplication mode, one column element is detected in the same way as in division mode.When R4j is accessed with the column number of local memory Cj, its contents are multiplied by the detected value. It is set to RlJ.The contents of R%5 and Rtej are subjected to Gaussian elimination calculation in Gaussian elimination mode.

After this operation is executed for the number of nonzero elements in column i of the lower triangular matrix, i is counted.

If the number of non-zero elements included in a row exceeds the number of local processors, each local memory is hierarchically divided into 1, for example, q==5. p=2. If r=3, the first
As shown in Figure 0, the local memory consisting of (Au*A1□#A+3) and the local memory consisting of (A21yA2□+A23) are moved to the third floor (r
=3) in a hierarchical manner, A, , Am, A3 in Fig. 5.
．． A4. A! − content respectively AI Hs Ax+
, A1□l A221 Just store it in Alff,
That is, with two (p = 2) local units (
This means that the case where the number of non-zero elements is 1=5 can be handled.

[Effects of the present invention] In the present invention, a system for configuring a parallel computer dedicated to LU decomposition of large-scale sparse matrices is described.
In systolic array processors that have been proposed for such purposes, the structure of the coefficient matrix Y handled is limited to a band structure. However, since the coefficient matrix of a general circuit network does not have a band structure, it is difficult to apply an array processor in practice.

Therefore, an LU decomposition processor has been proposed which processes operations on columns corresponding to each non-zero element in parallel by preparing q local processor units equal to the number of non-zero elements in one row of a matrix.

Each local unit is connected by a common data bus, priority circuit, etc. Only the non-zero elements of the matrix Y are stored in the local memory, ignoring all zero elements, and during the Gaussian elimination process, the contents of the registers are shifted while the data label numbers match. By doing this, it is possible to execute the calculation. Furthermore, some properties of the shift process were demonstrated, and it was shown that even if a method of matching label numbers by shift operations is used, it does not require much time for data communication. Finally, by simulation.

It has been shown that the efficiency of this processor is superior to that of a 7-ray processor.

In this manner, the present invention has the advantage that an nXn large-scale sparse matrix can be decomposed into I, U in approximately O(n) using processor units whose number is less than or equal to the number of nonzero elements in each row in one matrix.

[Brief explanation of drawings]

Figure 1 shows a local unit of a simultaneous equation analysis computer according to the present invention, Figure 2 is a system configuration diagram of the present invention, Figure 3 is an example of a graph, and Figure 4 shows the matrix structure for the graph. Figure 5 shows the data structure in the local memory (A2), Figure 6 shows the data structure in the local memory (Bj), Figure 7 shows the circuit configuration for left shift, and Figure 8 shows the data structure for one matrix Y. Figure 9 shows an example of the structure of the matrix Y according to optimization of the pivot order.
Figure 0 shows another embodiment of the local unit of the computer for simultaneous equation analysis according to the present invention. 1.2.3...Local unit to, 1G...Common bus 11/Fist/Control line 12...Priority circuit 13...Counter 14...Exclusive OR 15...Enable circuit A ,,A,1...first local memory B, B2.
...Second local memory P,,p,...Local processor Aoi20 澾3邑'! 414 Figure 1, Section C, Section r? 80 yakusen 90 arithmetic

Claims (4)

[Claims] (1) Element a_i_j (i=1, 2, ..., N; j=
The value and column of each element a_i_k of each column vector created by left-justifying non-zero elements including fill-in in each i-th row of a sparse matrix Y consisting of 1, 2, ..., N) in the row direction. The first one stores the number k as data at address i.
A local memory of at most q, the maximum number of non-zero elements including fill-ins, existing in each row of the matrix Y, and each j column of the lower triangular matrix part excluding the diagonal elements of the matrix Y. In the second local memory, at least the row number l of each element a_l_j of each row vector created by top-filling non-zero elements including fill-in in the column direction is stored at address i as data. Install a number smaller than the number of local memories, and in division mode,
detecting a pivot element a_i_i such that i=k among the respective elements a_i_k and the column number k read from each of the first local memories according to the address i;
Each of the read elements a_i_i is simultaneously divided by the pivot element a_i_i to obtain a_i_k/a_i_i.
and a division means for calculating, in the multiplication mode, reads the contents a_l_d of each of the first local memories using each row number l of each of the elements a_l_d read from each of the second local memories according to the address i as an address, l=i. Detect the element a_l_i, and multiply each value a_i_k/a_i_i obtained by the dividing means by the element a_l_i to obtain a_l_i・a_i_k/
multiplication means for determining a_i_i; and an element a read from each first local memory by said address i;
Each first column number k and row number l of _i_k
Using a shift means to match the column number j of the element a_l_d read from the local memory of the element a_l_d, each value a_l_i・a_i_k/a_i_
A computer for analyzing simultaneous equations, comprising a Gaussian elimination calculation means for subtracting i. (2) A computer for simultaneous equation analysis according to claim 1, characterized in that the shift means uses a random access memory whose contents are exactly the same and which executes sequential writing and simultaneous reading. (3) The simultaneous equation analysis computer according to claim 1, wherein the first local memory is accessed by the contents of the second local memory. (4) If the number p of processors that execute the division means, the multiplication means, and the Gaussian erasure calculation means is smaller than the maximum number q, the p first local memories are divided into r×
A computer for analyzing simultaneous equations according to claim 1, characterized in that the computer is hierarchically connected to the rth order so that p exceeds q.