JPS63265365A - Analysis of physical phenomenon - Google Patents

Abstract

(57) [Summary] This bulletin contains application data before electronic filing, so abstract data is not recorded.

Description

DETAILED DESCRIPTION OF THE INVENTION Field of the Invention The present invention relates to an improved method for analyzing physical phenomena.

[Conventional technology]

Many physical systems can be described mathematically using a system of linear equations. These linear equation systems are solved by matrix processing methods. Finite element analysis is a method concerned with describing various physical phenomena in terms of systems of equations and developing methodologies for solving those systems of equations. In this specification, the term "physical system" refers to structures, devices (domains, etc.).

nasal paratus), or a material body (solid, liquid,
gas), or simply a region of space in which a particular physical, chemical or other phenomenon occurs. Finite element analysis has its origins as a method of performing structural analysis, but today it is used to analyze electric motors, generators, magnetic resonance imagers, aircraft engine ignition systems, circuit breakers, and transformers, to name a few. It is routinely used in the design of

Its techniques are used to analyze stress, temperature, molecular structure, electromagnetic fields, electrical currents, physical forces, etc. in all kinds of physical systems. It has become a standard part of the design cycle for many products, and the present invention has particular use in the analysis and design of such products.

The systems of linear equations that finite element analysis techniques are required to solve are often so large that they are difficult to solve using computers. For example, a system of equations from a large, but not typical, two-dimensional finite element analysis may have 25,000 unknowns. If such a system of equations has contributions from a large number of nodes (eontribu
tlon), there is no choice but to use tremendous computational power to obtain the solution. However, in some cases, such equations are large and coarse, providing an opportunity to preprocess or transform them so that they do not require extensive computational power to solve. The term "coarse" as used herein refers to the characteristic that only a very small percentage of the elements in the matrix have non-zero values. When extreme roughness exists in very large systems, there are several techniques that can be used to transform the system of equations into one that is more easily handled from a computational standpoint. However, despite such transformations, the size and other characteristics of the resulting matrix equations may make standard computational techniques impractical and the efficiency may be very low.

As can be seen from the foregoing, the field of finite element analysis has advanced greatly due to the availability of larger, more powerful computing machines to solve such large systems of equations. There are now high-performance dedicated computer systems for each scale designed to perform specialized application calculations that are particularly common to perform on general-purpose computers. An example of such a dedicated computer system is a systolic architecture.
It provides a general methodology for mapping high-level computations to hardware structures. sygtolie computer system (sygtolie eonp)
uter sy+stem) Flows periodically from the memory of the Cheetaka computer, in a chain-! ! The memory is pipelined through many processing elements and back to memory, allowing each memory access to perform many calculations without increasing the associated input/output requirements. This results in significantly faster execution of problems that require calculations. Several methodologies for constructing contracted architectures for handling matrix operations are proposed by the Society for Industrial and Applied Matrix
Industrial and Applied Math
emties) 1979, Sparse Matrix Pr
oc, ), 1979, pages 256-282, H.T.

Kung) and C.E. Lazerson (C,E.

Systolic Play (Lsiaerson)
(Systolie A)
rray (for VLSI) J. Another analysis and suggested solution to this problem was published in the journal IEEE Trans.
actions on Computers) J, V
ol, C-32, 3. March 1983, "An Efficient Parallel Algorithm for...
The Solution on Large Sparse Linear Matrix Equation (AnEfficie
nt Parallel Algorithm for
the solution of Large 5pa
rse Llndar Matrix Equatio
ns) described in J.

U.S. Patent Application Nos. 870,467 and 870,56
The specifications of No. 6 disclose techniques and devices that use parallel processors to execute various matrix operations common to finite element analysis among various analysis methods and to increase calculation speed. ing.

This is primarily done by arranging the solution to be performed as a large number of repeated, nearly identical steps executed in parallel by a series of parallel processors.

For example, U.S. Pat. Disclosed. One such operation used to solve a system of linear equations is backward substitution (back 5u
bstition), but other operations will be obvious to those skilled in the art.

Solutions to linear equations generated as a result of the techniques of the present invention can be implemented on a variety of parallel multiprocessor architectures, as discussed below. However, for illustrative purposes, the solutions of the equations generated as a result of the method disclosed herein were presented at the IE International Conference held March 26-29, 1985.・On Acoustics, Speech, and Signal Processing (Part I)
1ons I Conference@nceon Acou
sticks, 5peech, and Signal P
``ASystolie Array Computer'' by E. Arnold et al.
Computer) J, pages 232-235. The method of the present invention is described in U.S. Pat. No. 4,493.048 [Sytolle Array Apparatus
uses for Matrix Computer
ions) J, which patent is frequently cited herein.

[Problem that the invention seeks to solve]

Backward substitution or forward elimination (forward @11
The main disadvantage of the method and apparatus disclosed in the aforementioned U.S. Pat. This arises from data dependence that forces

[Purpose of the invention]

, It is therefore an object of the present invention to obtain an improved apparatus and machine-implemented method for analyzing physical systems represented by systems of linear equations.

Another object of the invention is to operate a multiprocessor in such a way that manufactured products can be analyzed with higher efficiency and physical systems represented by systems of linear equations can be analyzed with higher efficiency.

Another object of the invention is to process a system of linear equations by a parallel processor.
In particular, it is an object to provide a method and apparatus for increasing the speed and efficiency of solving on parallel, pipelined processors.

Another object of the present invention is to be simple so as to be able to perform a large number of calculations simultaneously and to enable highly efficient implementation.
It is an object of the present invention to provide a method for processing large coarse matrices in a contraction play computer system using only periodic 2 communications and control.

Another object of the invention is to provide a "bandwidth-maximized" triangulated matrix, i.e., a preselected triangulated matrix that separates non-zero elements on the main diagonal from non-zero elements in the upper right or lower left corner of the matrix. A methodology for increasing the speed of computing solutions to systems of linear equations by using backward substitution or forward elimination techniques performed on transformed triangulated matrices to obtain bands of zero elements of minimum width. It's about getting.

Another object of the present invention is to obtain a technique for renumbering the nodes of a finite element network used in performing the finite element analysis method so that a decomposed mental matrix characterized by the above-mentioned bands is generated. That's true.

Yet another object of the invention is to obtain a parallel processor;
and to obtain a methodology for operating the parallel processor to solve systems of linear equations characterized by system matrices of the type described above in a new and more efficient manner.

These and other objects of the present invention are to provide a triangulated matrix for a physical system with a band of zero-valued elements adjacent to the main diagonal of the matrix, i.e., a predetermined minimum. The physical This is achieved by methods and techniques for performing finite element analysis methods to analyze systems. The existence of bands of zero-value data in the decomposed matrix as described above allows a novel way to perform many calculations of the unknown vector components of the system simultaneously and with very high efficiency. Computation of the values of the components is performed by backward substitution, independent of previous solutions of other values of the components of the unknown vector. The overall speed of solving the system is significantly higher since it is no longer necessary to delay the backward substitution method for a given vector component until the calculations previously performed for the other components of the vector are completed.

〔Example〕

Hereinafter, the present invention will be described in detail with reference to the drawings.

As mentioned above, the method of the present invention relates to improvements in finite element analysis methods used to analyze physical systems. In those physical systems, field variables (fl・ld
Such a field is The variables have infinitely many component values. In the first step of the method of the present invention, the problem contains a finite (but large) number of unknowns by dividing the solution domain into elements and representing the unknown field variables by approximate functions in each element. can be separated (dlsc)
r@tizsd).

In this way, complex problems are reduced to solving a series of greatly simplified problems. The approximation function is defined by the unknown values of field variables at specified next points called nodes. Matrix equations are then written that describe the properties of the individual elements.

In general, the outcome variable or forcing
g) Let the variable R be the field variable Y at a given node
Each node is characterized by a linear equation relating K to a "stiffness" variable K. More specifically, for a solution domain under investigation with a large number of such nodes,
The component ``R□, R3,...RNJ (known or having a zero value at the boundary node) result variable or forcing variable rRJ is replaced by the component [Y□, Y2...
・The vector field variable rYJ with ``YN'' is multiplied by the stiffness matrix rKJ, which is made up of constants or coefficients of the field variables at their nodes. Ru. For example, in the analysis of a linear spring system, the result variable components represented by (R□, R1,...) can be the values of the forces applied to selected nodes of the system. The variable components of (Yl, Y2... can be the displacement values at the nodes, and the coefficients can be the spring stiffness values of the spring element being investigated, which relate the forces to the displacements at the nodes. The constants form a coefficient matrix or stiffness matrix (K).

The next step in the finite element analysis method is to solve the above Marx equation. The subject of this application is a mechanical technology that handles matrix equations in order to quickly and efficiently automatically obtain solutions to systems or equations representing complex physical systems analyzed by finite element analysis. A more detailed explanation of the principles of finite element analysis is provided by K.H.
) Iuebun@r) [The Finite Element Method for Engineers (The Fi
nite ElementMethod for En
gine@rs) J, and Klaus; Jurg@n-Klaus Bathe
tProcedur@in Eng1ne@ring
Analysit) J.

Many techniques or methods are available for solving the above system of equations. The nature of these techniques in general is often very complex and involves the use of decompositions or transformations of systems of equations that are derived from and are equal to or approximate the original equations.

Therefore, for example, the original system equation can be converted to the form %expression%(). Here, [A] has dimensions NXN
is a known decomposed matrix, Q is a known transformed N vector, and an unknown transformed N vector.

The purpose of such transformation and decomposition of the original system of equations is to be able to apply various known matrix techniques to the transformed system in order to arrive at the final solution.

Backward substitution and forward elimination are two such numerical techniques used to solve systems of equations obtained using finite element analysis methods.

The backward substitution technique uses the decomposed or transformed matrix (
A) is used when the matrix is triangular (upper or lower) in which all diagonal elements are non-zero. Backward substitution techniques can be used once to arrive at a definitive final result, or repeatedly as part of a more complex approximation technique to arrive at a final solution. Techniques for converting linear systems of general form to triangular systems are well known.

A system of such triangles has the following shape:

A11Xt+A1gXs” “”l, N-lXN-1+
AINXN“QIAllXl”””1. N-lXN-
1+AINXN ” QIAN-I N-1xN-1+
AN-I NXN = QN-IANNxN= QN To solve this system of equations by backward substitution, ANN (
Note first that the last equation can be immediately solved for XN since Q9 (the diagonal element of the stiffness matrix) and Q9 (the component of the forcing vector variable) are known. If we know XNN, there is only one unknown in the penultimate equation, namely 4-□, so if we know XN, we can solve the penultimate equation.

In particular, . Since XN and xN-1 are known, the third-to-last equation has only one true unknown, i.e., XN-1
Since it contains, we can solve the equation. Therefore, it should generally be considered that i=N, N-1, . . . 1. This is interpreted as a sum with no terms and gives the value O by convention.

Several features of the backward substitution method allow for the large number of simultaneous processing possibilities realized by the multi-processing apparatus described herein. First, once calculated, XN must be multiplied by each element of the nth column of the matrix during the process of solving for the component values of X. Therefore, in the calculation of XN-□, XNKAN-1,N(N
- the coefficient from row 1 and column N). Similarly, in calculating XN-I N...Xl, ~
must be multiplied by AN-1, Nl...XIN (remaining coefficients of N columns).

Similarly, once calculated, the next term XN-1 is multiplied by each coefficient in the N-1 column.

An apparatus and method for efficiently performing the above process on parallel processors is disclosed in the aforementioned US patent application Ser. No. 870,566. However, the apparatus and method disclosed in that U.S. patent application are inherently limited in that the data dependencies in the operations being sought limit the available parallelism in solving finite element equations. Therefore, the processing speed of the multi-gross processor of this application is reduced. In particular, from the above equation, the previous X1+1
+Xl+1...XN cannot be computed until all of XN have been obtained, their associated column-wise multiplications have been performed, and the partial results have been accumulated, so XN-11XN-2...
The process of obtaining etc. is essentially serial.

It has become quite common to use decompositions and transformations that are performed when the physical system to be analyzed is characterized by large and coarse system equations (with a high proportion of zero coefficients).

Band*dness means confining the nonzero class within a particular portion or band of the decomposition matrix, and can be better understood with reference to FIGS. 5A and 5B. 5A and 5B, the coefficient structure of the two matrices, i.e., the fixed pattern (fo
otprlnt) is shown. Solid diagonal lines represent non-zero elements, cross-hatched areas represent areas containing zero data elements and non-zero data elements, and white or black areas represent only zero elements.

Thus, in FIG. 5A, the non-zero elements of the matrix are confined (along with a certain number of non-zero elements) to the width bands on either side of the diagonal. Only zero elements are placed in the upper right and lower left corners of the matrix. This type of summation is called ``band width minimization'' because nonzero elements are confined to small diagonally adjacent bands. As those skilled in the art will appreciate, such a matrix with minimized band width reduces the number of calculations required to perform the backward substitution method. The reason is the known calculation (blank (b
This is because it is inevitable that the number of times (including all elements in the rank) area will yield a zero result. However, since the unknown component values of the field variables must be calculated in series, the data dependence problem cannot be solved by minimizing the band width.

The formation of a system of equations with a minimized band width and its corresponding minimized decomposition matrix is such that the difference in the numbers assigned to adjacent nodes in the mesh of the finite element system is minimized. is accomplished during a finite element analysis method as a result of known techniques for numbering or renumbering nodes within a finite element system. Current techniques for numbering nodes (including appropriate programs) to achieve band width minimization are discussed in detail in the following publications: Earl Roces (R,
Ro■・n) [Matrix Bandwidth Minimization (Matrlx BandwidthMl
nimlzation) J., Proced.
edings.

National Conf@rence A, C, M
, ) # 1968. pages 585-595, and E.H. Cassill (iH.

Cutllll) and G.M.
J.M.

M6 @e), rReducing the Bandwidth on Sparse Symmetric Matrix (R
Educating the Bandwidth of
5paras SyrrnatrlcMatrlce
s) J. Proaaedings, National Conference N.C.

National Conference A, C, M
), 1969. pp. 151-172.

To solve the finite element equation, remove the above data dependence,
In order to solve these equations more efficiently on parallel processors, we decided to "maximize the bandwidth".
A proposed technique is used herein, ie, the use of a decomposed matrix with a coefficient structure as shown in FIG. 5B. Such a matrix (and its corresponding system of equations) has non-zero data on its main diagonal, a band of only zero data elements adjacent to the diagonal, and other bands grouped in the upper right and lower left corners of the matrix. all data (including all other non-zero data). It is easy to create a matrix with a Q1-maximized pattern by numbering or renumbering the nodes of a finite element network so that adjacent nodes differ by at least a certain minimum value. It will be done. Therefore, there are only zero data elements in the KICK bands on either side of the matrix diagonal.

When a system of equations characterized by a decomposed matrix of the form shown in FIG. High efficiency can be obtained compared to other devices. These efficiencies derive from the modification of data dependencies when performing backward assignment or forward elimination in parallel architectures according to the method of operation of the present invention.

The backward substitution to solve for the first component of the unknown field variable now looks like this:

and similarly, etc.

The reason for this is that, as shown in Figure 5B, the generalized bandwidth of the solution matrix is maximized in the ninth AN-
1,N' 1AN-1, da and AN-j, N-0 are all zero. Generally 1 If a matrix with maximized bandwidth is adopted, backward assignment is 1=N, N-1
,..., for l, X, = (Qi-Σ Xj
AI, j)/Al, ij-gl +1 will be performed. Therefore, any given X, such that j-1≧, can be computed independently of the previous KXJ. For example, in performing the backward substitution method using a bandwidth-maximized matrix described below, all the data on which a given K is chosen to correspond to the longest time required to obtain

A brief description of matrix algebra as a computational device for solving systems of linear equations may be useful as background for understanding the present invention.

The most useful form of matrix is a rectangular array of m rows and n columns of scalar fields arranged orthogonally to each other. Since the order of a matrix is determined by its number of rows times its number of columns, the matrix shown below is called an mxnJ matrix. The usual concise way to represent a matrix is shown below.

The rows of a matrix are horizontal lines or one-dimensional arrays of quantities, and the columns are vertical lines or one-dimensional arrays of quantities. The quantities A□□, A□3, etc. are matrix (A
) elements. A matrix with n-1 is a column matrix or column vector J), and a matrix with m=1 is a row matrix or row vector. In general, a vector is an ordered value consisting of n elements (ord@red n-tup).
l@of value@a). Lines in a matrix are rows or columns.

A square matrix is a matrix in which the number of rows and the number of columns m are equal. The diagonal of the square matrix is composed of elements A1□1AllA88...ANN.

A triangular matrix is one that contains all its nonzero elements in, above, or below its main diagonal. The upper triangular matrix has all of its nonzero elements in or above its main diagonal, and the lower triangular matrix has all of its nonzero elements in or above its main diagonal. Has it on the bottom.
・If [A] is mxn order and (X) is HXp order, (A
) In the operation of X(X) - (Q), (Q) is m
It is a matrix of order Xp. Therefore, if [X] is a column vector, then [Q] is a column vector with the same number of rows m as the number of rows in [A]. It should be noted that matrix multiplication is defined in the above example only if matrix (A) has the same number of rows as the number of columns of matrix (A). In the above multiplication, in general, [Q] is obtained by multiplying the column of [X] by the row of [A], such that Q,j=Σ AlkxXkj k■l .

According to the above, vector (X) to matrix (A)
The multiplication of is normally performed as follows.

Due to the processing time required, it is most efficient to perform all multiplications involving each specific vector element simultaneously, as this reduces intermediate memory accesses. It is true. Therefore, all operations Al lXl ! A11Xl' and A8
It would be more efficient from an input/output standpoint to perform the above multiplications by doing □X0 simultaneously, ie, in parallel. However, this approach allows us to identify various contributing elements of rQJ.
This results in scattering.

Those elements must be collected later, i.e.
have to be recombined.

The second feature of the above multiplication method is that each element Ql of the product vector
+QB +Qs is the result of an accumulation or addition based on the original row (or row index (1nd@x)) of the matrix element. Specifically referring to the example above, Q
Note that □ is the sum of the partial results obtained by multiplying the various elements from row 1 of the matrix by the relevant or corresponding elements of the vector (X).

Thus, the main observation used in the methodology described in this specification is to increase the concurrency of operations and minimize the number of input/output requests in the multiplication. This means that we should temporarily proceed as follows.

First operation → X1* (A1□ Afil A81
) Second operation → X, * (A, A,, A,
, ) 3rd operation → X, *(A,, A,, A
38) Following the above operations, add the partial products as follows based on the row indices of the matrix elements included in these partial products.

A more general opinion, concluded from the above, is to perform a first common operation across a series of different elements emanating from one line of the matrix, followed by the ultrathickness ( By performing the second operation based on ``orlgin'', many matrix operations can be performed in a highly parallel manner. Matrix multiplication is one of those operations.
It is one.

More specifically, from the above explanation, the first operation in the multiplication process is the first element K of the multiplication vector.

Element A□□IAm l which is an element of a matrix column
The second vector containing IAs□ is multiplied.

In this way, multiplication is performed by assigning each element of the multiplication vector to the transposition vector (tranrp
This is done by multiplying ose v@ctor). Vector elements X and l have elements Al I I Am
11'41 is similarly multiplied by the transpose of the second column of the matrix containing 11'41. The result of each such multiplication includes partial results. Assuming that those partial results are accumulated based on the original rows of the stored matrix elements, the result spec.

This results in the generation of torque. This means that the original row
orlgin), i.e. NA11XllA11
It follows from the previous explanation that Q is generated by accumulating partial results based on XllAllX8.

The same applies for the remaining elements of the result vector.

It turns out that the above matrix multiplication can proceed simultaneously among other multiplications as above, but
The following observations apply to large and coarse matrices encountered in finite element analysis.
on) remains to be added.

For such large, coarse matrices, zero elements can be discarded or ignored in the procedure for mapping the matrix to storage, since they do not contribute to the final result. Although zero elements are always ignored or discarded when mapping the matrix to be processed onto parallel processors, it is equally possible to ignore or discard other elements that are not exactly zero but are not important in the computation. .

For example, once it is determined that including elements that fall near zero does not change the final result, those elements can be ignored. As a precaution, the term zero should be interpreted in this sense.

Once it is recognized that matrix operations have the potential to be performed in many parallel ways, it remains to develop techniques for mapping matrices onto parallel processors in a way that allows optimization of parallel processing. ing. Such mapping methods are described in the aforementioned U.S. Patent Application No. 870,46.
This is the subject of issue 7.

Reference will now be made to a parallel processing apparatus having a network of interconnected processing cells each having an associated memory and in communication with each other, as shown generally in FIG.

Reference is made to FIG. 1, where a basic deflation device is shown. In the contraction device, a host device 10 receives data;
The results are output via inter-78 unit 12 to a contracting array of pipelined processing cells 15.

Of course, host device 10 can include a computer, memory, real-time devices, and the like. As shown in FIG. 1, the 'input and output terminals are coupled to a host device 10. Input terminals can be from one physical host and output terminals can be directed to another physical host. An important goal of the contraction architecture is to perform a large number of computations on the host for each force/output access of the shear array. A large number of those calculations are
When cells are configured to continue processing received input in a "pipeline" fashion, processing is done in parallel. In this way, while a cell is receiving and processing new data, other cells continue processing data or partial products of data received in a previous input access. The challenge of using a shrinkage device is to divide the problem to be solved into subprocesses that can be processed periodically in parallel by different cells of the device.

A procedure for loading and storing a large coarse matrix into a series of parallel processing cells is described in U.S. Pat.
No. 457. In order to fully explain the invention, the procedure will now be repeated with reference to FIG. FIG. 2 is a flow diagram of the contraction play architecture shown in FIG. Matrix loading and storage takes place under the control of the interface unit 12 shown in FIG. The interface unit 12 communicates with the host 10 during its processing operations. The host already has the matrix stored in an array or other suitable data structure.

As mentioned above, the purpose of the mapping procedure shown in Figure 2 is to store a matrix in the memory of a multiprocessor so that a very large number of matrix operations can be performed simultaneously in general, and in particular a large number of backward assignment operations can be performed simultaneously. and store it. This is done by storing the elements of one line of the matrix, allowing simultaneous operations with common operators. More specifically, as will be explained later, in the 1-backward assignment operation, the elements of the matrix sequence are continuously assigned by a common operator that is the associated and previously 1-increased component of the variable vector of the field to be solved. Stored for multiplication or simultaneous multiplication.

A second feature of the storage method is to generate a ninth index that is associated with each stored matrix element that identifies the origin (relative to another line of the matrix) of each element. In particular, each element from a column of the matrix has
It is combined with an index that identifies the row of the matrix from which it came. Additionally, the stored elements from each column are formed into a permutation vector. Common operators can be combined with each permutation vector to perform operations regularly and simultaneously. Several permutation vectors are thus formed. Because each permutation vector is placed in a given memory location within multiple cells, a common operator can operate on the permutation vectors simultaneously. The origin and character of the common operators are highly dependent on the particular matrix operation being performed. Also, invalid elements, zero elements, and computationally unimportant elements are discarded during the store operation so that mah IJ lux calculations involving the stored elements can be performed more quickly. Finally, given any permutation vectorization, there is a diagonal element contained in any given permutation vector that is easily available for periodic and iterative processing during matrix operations (i.e., are collected for storage in one processor cell.

The above object is achieved by converting the elements of the matrix columns into permutation vectors constructed across the processing cell. Since zero elements of the matrix are ignored during the storage operation, the number of processing cells required to implement the above storage method is quite reasonable. A large number of physical systems can be represented by a matrix of the order of 25000x25000. In those matrices, the number of nonzero elements in any given column is on the order of 20 or less. Therefore, as a result of the above storage operation, each non-zero element of a column of the matrix is a permutation vector or array stored in the same memory location of each processing cell. Also, all diagonal numbers are stored in one processor cell for ease of access. Due to the non-zero elements of the column that can be accessed simultaneously by this procedure), a common operator can be applied to the permutation vector in a number of parallel operations, as will be explained in more detail below.

Next, a method of mapping a matrix to a multiprocessor array in order to achieve the above objective will be described with reference to FIG. Starting at step 200, a constant rCJ is set equal to the total number of processing elements in the parallel processor array to be utilized. The number of processing cells can vary, but efficiency is highest when the number is approximately equal to the maximum number of nonzero elements in any one column of the large coarse matrix to be mapped. Procedural variable “ROWJ (
The row index for the matrix element in the processing field) and [C0
Since LJ (column index for the matrix element being processed) is set equal to 1, processing starts at the matrix element in row 1, column 1 and proceeds from there. The variable 'PROcJ (the next processing cell to load the non-zero, non-diagonal element) is set equal to 2. Of course, l'-PROCJ can vary between 1 and C, but later As will be explained in detail, cell 1 is marked to contain only the diagonal elements of the matrix.The memory variable (the memory location in PROQ being loaded) is set to [MEMJ. M.E.
MJ is the first memory location to be placed in each processor cell.

The host sets variables M (number of rows in the matrix) and variables N (
(number of columns in the matrix) is also set, or has already been set by a previous operation. Finally, step 200
A table marked table (COL) (this table is provided in the host) is used to store one column of the matrix (the particular permutation vector to be created) into a particular column to be used later in the matrix operation. All are initialized to 1 in the preliminary step of correlating operators. This initialization will be explained in detail later.

Next, in step 201, the matrix element person (
A determination is made as to whether ROW) (COL) (whose element is matrix element A(1)(1) at the beginning of the procedure) is zero. If the element is zero, the element is discarded without being loaded into any processor cell. Therefore, the process is step 2
Proceed to step 02 and increase the variable ROW in that step. Processing then proceeds to step 203 where it is determined whether the variable is greater than M. If the variable is not greater than M, processing returns to step 201 to determine a new matrix element.

If the variable ROW is greater than M, it is determined that the matrix is a complete column, and the process proceeds to step 205.

In step 205, zero is stored in cells in memory locations MEM of any remaining processors from PROC to C. The purpose of this step, as described above, is to enable each particular memory location within each processor to contain stored elements associated with one and only one column of the matrix. be. In this way, the elements associated with a given column contain a permutation vector (each stored in a known memory location across multiple cells) and multiple such permutation vectors (one for each column of the matrix). ) is created.

After all processing cells in a particular MEM are filled,
Processing continues scanning the matrix starting at the top of the next column. Therefore, in step 207, the variable ROW is reset equal to 1 in preparation for testing the matrix columns in the new column of the matrix, and the variable C
Increase OL. However, if it is determined in step 209 that the variable COL is greater than N (the number of columns in the matrix), the process of loading the matrix ends (step 210). Also, if the variable COL becomes greater than N, processing returns to step 201 where a new matrix element is tested.

If it is determined in step 201 that the matrix element is not zero, processing proceeds to step 230 to determine whether the element is a diagonal element, i.e., ROW-C
Determine whether or not you are an office worker.

If the element is a diagonal element, the element is PROC-1
(processor cell 1), an exponent is generated,
It is accompanied by a stored diagonal element indicating the row of the matrix that generated it. After storing the diagonal element in cell 1, processing proceeds to step 202 and to subsequent steps similar to those described above.

If it is determined in step 230 that the matrix element is not a diagonal element, the matrix element is stored in a memory location of processor cell PROC (step 220). An index equal to ROW is then generated and stored in association with the newly stored element in the same memory location MEM (step 221). , the process increments PROC" after storing a non-zero element (step 225
), the next non-zero element is stored in the next processor cell, and then a determination is made in cell 227 whether its PROC is equal to C (the total number of processor cells). A positive determination indicates that all memory locations within the processor cell are filled. Processing then proceeds to step 228 where PROC is reset to 2 and memory is increased to the new value. Along with that, the entries in the previously created table are increased to allow the host to correlate each column of the matrix to the given length of the permutation vector during later processing. . This factory allows the host or interface unit to associate the non-zero element corresponding to this current COL with a permutation vector covering more than two memory locations (the exact number is included in the table). be able to do so. If a large number of memory locations are required, the host takes appropriate actions later in the matrix arithmetic procedure so that the corresponding operator is provided to the processor the correct number of times. Processing then proceeds to step 202, as noted above.

If it is determined in step 227 that PROC is not greater than C, processing proceeds directly to cell 202 since additional memory location MEM is present in the processing cell.

From step 202, processing continues as described above.

Referring now to FIG. 3, a particular upper triangulated matrix [A] (shown at the top of FIG. 3) whose bandwidth is maximized in accordance with the present invention is An example of the operation of the procedure shown in FIG. 2 will now be described to map to exemplary memory locations of processor cells 1-4 of a multiprocessor similar to the multiprocessor shown in the figure. As can be seen from the figure, the decomposed matrix (A) is of order 12x12 and has the shape of an upper triangle. In particular, this matrix has only zero values below the diagonal, and all diagonal elements are non-zero elements. Also, all non-zero data in matrix (3) are grouped together in the upper right corner, and a band the width of the zero element is provided between the diagonal and the data group in the upper right corner. In this way, non-zero data in the matrix are placed diagonally or in the upper right corner, with only zero elements placed between them.

Techniques for creating such a bandwidth-maximized matrix and selecting the value K will be described in detail below. When mapping this specimen, it is meant to be an example of an upper triangular matrix that maximizes the bandwidth of the large, coarse matrices commonly encountered when characterizing physical systems using finite element methods. You should remember one thing. The typical order of such a matrix is 25000x25000, with rows/
Only 0.1% of the elements in the column are non-zero.

Starting at step 200 of the flowchart shown in FIG. 2, constants and variables are initialized to match the properties of the matrix (A) to be mapped as follows. C-4, ROW-1, C0L-1, N-12
, M-12, memory mL-MEM, PROC=2 is set to 2 to set aside the first processor cell to hold the diagonal elements, as will be explained in more detail below. In step 201, the first element A(ROW)(
COL) is zero. Since its value is not equal to zero, it is then mapped in step 230.
IJ Lux diagonal element or not, especially RO
A determination is made as to whether it is W-COL. Since element 1 is a diagonal element, the element is stored in location of processor cell 1 with index 1 indicating its superthickness from row 1 (step 231). Processing then proceeds to step 202 where ROW is increased. Since ROW is small by M, step 203
Then return to processing step 201 and perform the next element AC2]
Determine [1]. Since there are no other non-zero elements in the first column of the matrix, the process repeats steps 201, 202, and 203 until there are more than 12 rows. If the row becomes larger than 12, processing proceeds from step 203 to step 205. In step 205, a zero is stored in memory location MEM of processor cells 2-4 to indicate that column 1 of the matrix has no non-zero elements other than the diagonal element stored in cell 1.

The loading of these zeros into memory locations MEM of cells 2-4 is shown in FIG.

After scanning and loading column 1 of the matrix,
COL is increased to 2, the row is reset to 1, and CO
If L is smaller than N, scanning of column 2 is started (step 209).

Columns 2 to 7 are scanned in the same way as column 1, and each column also has only one nonzero element B-G, and since their nonzero elements are diagonal a elements, their nonzero element is cell 1
All are stored in the memory locations MEM+ 1 to MEM+6,
Zero is stored in MEM+l of cells 2-4.

Processing then returns to step 201 to scan the eighth column of the matrix. Therefore step 20
1, A[8](1) is determined.

Since its value (T() is non-zero and not a diagonal element (step 230), it is stored in memory location MEM+7 of cell 2 in step 220.

An index is also stored with rI(J in the same memory location in cell 2 to indicate its superthickness in the first column of the matrix. PROC is then incremented and determined in steps 225, 227.

Since PROC is less than 12, the row is incremented, determined, and the new element [8][2] is processed. Since this element is zero, the process returns to step 201 again and element A[8][3] is determined. Zero is 8 in rows 2-7
Because they are arranged in columns, processing is returned sequentially, but data is not stored in between. Since element A(8) (8) is nonzero and a diagonal element, it is stored in the same way as the previous data element was stored in cell 1.
As shown in FIG. 3, it is stored in cell 10MEM+7. Since all other elements of the matrix in column 8 are zeros, a zero is stored in memory location MEM+7 of cell 3.4 in step 205 in preparation for scanning column 9 of the matrix.

Processing then continues iteratively as described above, storing each nonzero diagonal element of the matrix in a different memory location in cell 1, storing the nonzero elements in each column in cells 2-4, and Discard zero elements and finally store a zero in any given memory location of the remaining memory cells before scanning a new column. As a result, all non-zero elements of the first column of matrix A are transformed into permutation vectors across the processor cell at memory locations identified as MEM. Furthermore,
All diagonal elements contained in such a permutation vector are stored in a particular processor cell 1. In fact, the memory location MKM acts as an identifier for the newly created, i.e. permuted, vector containing all non-zero elements of column 1 of the matrix, the diagonal elements being stored in cell 1. . In the case of column 1 of matrix (A) in FIG. 3, only one element is contained therein. Since it is a diagonal, it is stored in cell IK and all other elements are set to zero values and the remaining storage locations contain zeros as shown in FIG. The importance of the newly created permutation vector MEM with special treatment of diagonal elements in allowing simultaneous processing of matrix elements is described in the aforementioned US patent application Ser. No. 870,566.

As stated earlier, the number of processing cells is chosen equal to the maximum number of nonzero elements in any given column of matrix A, but may be 2 or 2 or more without departing from the spirit of the invention. Can be equal to any large number. However, regardless of the number of processing cells used, the number of non-zero elements in a particular column of the matrix may be greater than the number of processor cells.

To deal with this situation, step 227°22
8 is included in the mapping process shown in the flowchart of FIG. For example, if the first column of matrix (A) has more nonzero elements than the processor cell has available, i.e., if PROC is greater than C before all elements of a column are stored, Then (step 227
), its PROC is reset to 2 and MEM
+11c (step 228) and “MEM
The permutation vector identified as J continues to be stored in the MEM+1 memory location, ie, expanded.

As a result of this, column 1 of the matrix is moved to memory location ME
It is converted into a permutation vector located at M and MEM+1.

In that case, the table (COL) entry in column 1 is "2" to indicate this.

After determining all elements from a given column of the matrix, any remaining processor cell corresponding to a given memory location in which no non-zero elements are stored is determined in step 20.
Store zero according to 5. For example, matrix (A
) has only three non-zero elements M, N, and P in column 10, and zero is still stored in MEM+9 of processor cell No. 4, completing the transformation vector in MEM+9. A similar situation is shown in FIG. In this figure, column 4 of the matrix has three nonzero elements (S
.

R, T) exists, so as a result of step 205, zero is stored in memory location MEM+10 of processor cell number 4.

The function of that portion of the memory location, designated q in FIG. 3, will be described below in connection with performing a backward substitution operation to solve the system of linear equations containing the matrix A of FIG. .

A decomposed matrix whose bandwidth is maximized according to the invention
The execution of a backward substitution operation on a decomposed system of equations characterized by IJTxA will now be described with reference to FIG. To briefly recall the explanation given earlier, it is important to note that before the start of the backward assignment operation, the product or system of equations describing the physical system being investigated is generated in the form (K)(Y) = [R). should be understood. From what we have explained above, (Y) is a vector representing an unknown field variable that describes the properties of the physical system being investigated, such as the displacement of a spring system, with components Y,...YN Consists of.

The system of equations is solved for its components. On the other hand, the vector (R), like the force in the spring system being analyzed,
Known components R0, representing the resulting vector variables
- Has an RN. (K) is a stiffness matrix, as described above, containing values or constants relating known result vectors to unknown field vector variables at particular nodes of the finite element network.

According to the invention, the above system of equations is replaced by another system of equations (A) (
It is converted or decomposed into X)=(Q). Here, [A] is the matrix with the maximum bandwidth as described above! 0,
[Q] is a known vector, and [X] is an unknown vector.

Once the system is transformed into this new bandwidth-maximized form, our solution can be applied to a parallel processor with the structure shown in Figure 4 (at least for some parts of the solution).
A backward substitution technique or a forward substitution technique can be performed on the contraction play. Therefore, the purpose of this solution is
Assuming that [A] and (Q) are known, the value of Xl...XN is calculated by backward substitution technique. The decomposed matrix A, in its basic form, and as mentioned above, is usually large and coarse. especially,
The matrix may have elements of the order of 25000x25000, of which the number of non-zero elements only accounts for a very small percentage of the total number of elements. As shown in FIG. (containing zero elements) are grouped in the upper right corner of the matrix, and a wide band is placed on the diagonal containing only zero value data. The matrix in Figure 3 (
A) is just such an upper triangle bandwidth-maximized matrix.

Once the upper triangle bandwidth has been maximized, the first step of the solution is explained with reference to Figures 2 and 3. , to load and store matrix (A). The results of loading and storing matrix (A) shown in FIG. 3 are shown in FIG. In FIG. 4, memory 56 of cell 1 stores diagonal elements A...Y in locations MEM-MEM+11 of memory 56, respectively.

Each stored diagonal value is accompanied by an index stored in memory portion 59 indicating the superthickness of the row of diagonal elements associated with it. Calculated value of vector X X1...XHl
A partial result store 57 is also provided in cell 1 for storing . The operation of the partial result storage device will be explained in detail later.

Similarly, memory 56 in cells 2 through is shown storing the remaining non-zero values of the matrix. Each stored non-zero element contained in memory portion 58 is accompanied by an index stored in portion 59 to identify the row of the matrix in which the associated non-zero element is generated.

Each cell 2-4 also includes a storage device 57 for storing partial results accumulated according to stored indices with stored matrix elements, as will be explained in more detail below.

The elements of each column of matrix A are permuted during the mapping process to form what can be called a "permutation vector" in each memory location MEM-MEM+11.
You should also be careful. The diagonal elements are collected in cell 1. For example, column 9 of matrix A is memory location ME
Converted to a permutation vector at M+8, the diagonal elements from row 9 of the matrix are placed in cell 1, with the remaining non-zero elements (J, K) in column 4 indicating their respective rows of the matrix from which they occurred. They are included in cells 2 and 3 together with index (12), respectively. As mentioned above, zero data values are used to store unfilled locations in each permutation vector that remain after all nonzero elements in the associated column have been stored.

Given the triangular matrix (A) in the upper part of Figure 3, the system of equations to be solved for Cx) is as follows.

7・' (1) AX1+... ...to Shizuku+
...=Q engineering (211%+...
... Pigeon → [1, + [1□=Q.

(3) α, +......
Hara, -■1□ =Q.

+41 [X4+...
...■, -to, =Q.

(5) IX, +・-S −VJXl
, =Q.

+6) F′x6+−−−−−−
= Q.

(7) α, +...
...=Q.

+8) IX, +...
...=Q.

+91 LX9+Yu...
= Q.

(10) P)'1
0+...=QIQ(11)
Txll””” = Qt, (
12) y
x,=Q.

As previously mentioned, in order to solve the above equation system by the backward substitution method, X by solving equation (12).

First calculate. Similarly, in a similar way, the diagonal elements T, P correspond to the related known components of [Q].

Simply divide by L, I, G, and F to create X□1. X engineering. ,
X, +x8+X, , X6 can be solved at once. It should be noted that the ability of the present invention lies in the ability to calculate several components of the unknown field variable [X) simultaneously without regard to any other values of the field variable. This)
, as will be explained in detail later, can greatly increase the computational speed of parallel processors. This means that the field variable 1
This is clearly different from the conventional method of calculation disclosed in the prior art and pending U.S. patent application, in which the value of the remainder is always calculated and the subsequent values are known, except for the calculation of one component. Contrasting. This data dependency problem significantly slows down the speed at which these equations can be solved by backward substitution techniques on parallel processors, as will be readily apparent to those skilled in the art.

As mentioned above, we calculate the value of It can be used to calculate the values, ie, X, ~X in equations (5) to (1). It should again be noted that the above system of equations is an example of a very large and coarse system encountered in typical physical problems studied using finite element analysis methods. . The above system of equations is here for illustrative purposes only, and the use of a maximized bandwidth stratified matrix of the form shown in FIG. It should again be noted that this is the ninth detailed description of the resulting invention.

The actual method of solving the above system of equations by backward substitution techniques in a contracting multiprocessor system will now be described in detail with reference to FIG. 4 and FIGS. 4A-4L. This multiprocessor device is composed of four processor cells 1-4. , each processor cell is assigned to perform a particular operation that contributes to the solution. Each operation performed by a processor cell is represented in FIG. 4 by a functional unit. The functional units illustrated in FIG. 4 do not refer to the various hardware components of the machine, but rather to the various data processing operations that can be performed by each cell of the device, as will be readily understood by those skilled in the art. It should be remembered that it refers to The operations performed by each cell are arranged to be periodic and repetitive during each cycle of the machine. The flow of information through the processor will now be described with reference to Figures 4A-4L. The figures illustrate the flow of operators and results through the machine during the first 12 cycles of a typical backward substitution method for the above example system of equations.

The multiprocessor device shown in FIG.
It can be seen that the memory 56 associated with .about.4 is included. Memory 56 can be thought of as having a first array 58 that stores the values of the elements of solution matrix A, as described above. Each element stored in memory array 5B has associated with it an index indicating the row of the matrix that is stored in array 59 and yields the stored value associated with it. Thus, looking at cell 1, matrix element 1 arises from row 1 of the matrix, matrix element 1 arises from row 9 of the matrix, and so on. Storing in cell 1 means that the elements being stored are diagonal elements based on the storage technique shown in FIG.

Each cell also includes another storage device 57. This storage device 57 stores the results of the indexed calculations performed by the device, as will be explained in more detail in the description of the calculation process that follows. The storage device 57 functions as a storage device for storing the last component values X1 to X1□ of the unknown field variables. The storage devices [57 in cells 2, 3. gives the stored value. For this reason, the memory devices 57 of cells 2 to 4 are
produces two different outputs. The first output provided on line 57' is provided to a series of accumulators 53. In these accumulators, different values associated with the same index are added together before being recovered in memory in order based on that index.

The output on line 5r represents the value stored in memory associated with the diagonal index and processed simultaneously in accumulator 53. For this reason, the output at 11157' is controlled by the exponent input provided from memory 59 to line 5g. A second output from each memory device 5T of cells 2 to 4 to line 5r is provided to adder 65. The indexed output on line 5r represents the accumulation QN' of the previous calculations. Their accumulations are added together when performing the backward assignment method. The output on line 5r is selected by the simultaneous output to each memory 5T of cells 2-4 from the counter 61, which starts from index X1 during cycle 10 and counts down to index 1 during cycle 11 (4A-4L). figure).

The values stored in memory 58 of cells 1-4 are provided to the functional units in the remainder of the device associated with each cell. The values stored in memory 58 are applied step by step from the top or highest numbered address (FIG. 4) via line 58' to functional unit 52.51. This step-by-step application means that the value stored in memory location MEM+11 of memory 5B of cell l is applied to divider 52 in a given machine cycle. The given machine cycle is
The corresponding value dfi11 at the same memory location in the memory 58 of cell 2.3.4! 58' before being applied to the associated multiplier 51. This will be explained in detail later. Cell 1 includes a subtracter 60 and a divider 52. The subtracter 6° receives the first signal from the adder 65.
is supplied via line 65'. The second input to the subtractor 60 is the value of one of the known vector components in the physical system Q1-Ql,l. Divider 52 divides the result of the subtraction performed in subtractor 60 by the appropriate value from cell 1's storage. It is stored (ma)
IJ socks are the selected diagonal elements.

Therefore, each result generated from the divider 52 of cell 1 represents the solution of one component of the unknown field variables of the system of equations being solved. According to the backward substitution method, those results are
It is stored in the memory device 57 of cell 1. Also, the divider 52
The result of each division in M) must be multiplied by a permutation vector representing the corresponding column of IJ socks. For example, X12 must be multiplied by each nonzero element from column 12 of the matrix. Collecting the results of those operations,
Must be added later according to the row index.

For this purpose, each result from the divider 52 is fed in a pipeline manner to the multiplier 51 of cells 2-4. Indexed results are calculated in those multipliers. The output of each multiplier 51 is accumulated in an accumulator 53 according to the row index, and the result of the accumulation is stored again in the storage device 5T according to the row index, and is sent to the adder via #5r under the control of the counter 61. 65.

The apparatus shown in FIG. 4 will be explained in detail by explaining the operation when solving the above-mentioned simultaneous equations by the backward substitution method. In this description, for this purpose we will use the fourth A
~4L Refer to the diagram.

Generally, the processing apparatus shown in FIG. 4 operates in a pipelined manner. At the beginning of the calculation process, the storage device 56
The values stored in and the known values Q~Q12 can be used. These values Q0-Ql11 are stored in an associated processor, such as host 1 (FIG. 1), for proper input to subtractor 60.

It will be done. In a subsequent machine cycle, a selected one of those values is made available to the multiprocessor cell for generating a result operator. Processing proceeds by moving data through the processors generally from left to right in FIG. 4, gradually filling the processors until each cell performs a significant calculation during each cycle. The first few machine cycles are performed while simultaneously computing the first few values for the components of the unknown field variable (X).
, fills the pipeline as described below.

Referring first to Figure 4A, during the first machine cycle the value Q
, is provided to one input terminal of subtractor 600 from a suitable storage device or host computer (not shown).

Another input Q given to the subtracter 60 from the adder 65,
゛ is zero at this dark spot in the process. I mean,
This is because the operators to be acted upon do not reach the other functional units of the device at this point in the backward assignment method. Inputs to all other functional units are all zero.

During machine cycle 2 (Figure 4B), the equation 12
is solved for X engineering. This is done by dividing the output Q1m from the subtractor 60 by the first stored diagonal element YK of the matrix A, and then returning the first solution X1m1 to the divider 52.
This is done by causing a signal to appear at the output terminal of. At the same time, the output 901' from the adder 65 during one cycle 20
is subtracted from the value Q1□ in the subtracter 60 to provide a new input to the divider 52 for the next cycle from the value Q1□ in the subtracter 60 to the output terminal Q□□ from the adder 65'. Deducted. Also, the value X1fi is line 6
4 and is applied to multiplier 51 of cell 2 via line 63.

During cycle 30 (Figure 4C), divide the value Q□1 by T,
By storing them in memory 57 of cell 1 via these @64, equation 11 is solved for X□1. Also, during cycle 30, the value W (one of the matrix elements from the 12th column of the matrix) is added to the value
Multiplied in multiplier 51, accumulator 530 of cell 2
One input terminal has a partial product W-X1. to occur. Also, during cycle 30, the output of adder 65 (which is still zero) is the value Q in subtracter 60 in cell 10 subtracter 60. will be deducted from

In a somewhat similar way, during cycle 40, divider 52
Equation 10 is solved for X1o by dividing the υ value Q1° by the next stored diagonal element P inside.

At the same time, in cell 2, the previously calculated
The previously generated partial product W-X, is multiplied by 1 K to generate the partial product S-X, , and stored as Q,1 in the next cycle in the accumulator 53 of cell 2. is accumulated with another Q, partial product (at line 5T). In addition, the X-manufacturer, which has passed the one-cycle delay 55, is cell 3.
is reached and multiplied by the matrix element V to generate a partial product. The partial products are accumulated in subsequent cycles,
It is stored in the memory device 57 of cell 3.

During cycle 5, the output Q, + (still equal to zero) from adder 65 is subtracted from the value Q8 in subtractor 60, and Q, (the output from subtracter 60 in the previous cycle) is subtracted from divider 52. is divided by L to generate X. The previously calculated X. is stored in the storage device 57 of cell 1, and multiplied by N in cell 2.

The previous output SX1 from the multiplier 51 of cell 2 is similarly indexed and stored in the storage 57 (Q,' (if any) in #i57' and the accumulator 53 and stored in memory f157 in the next cycle.
is stored in The previously calculated Xl, is multiplied by R in the multiplier 51 of cell 3, and v-Xls is similarly indexed to the previously stored value Q,′ (if any) and the accumulator of cell 3. It is accumulated in the calculator 53 and stored in the storage device 57 in the next cycle. Finally, X
o, moves to cell 4 and is multiplied by U in multiplier 51, and the product M-X1. is sent to accumulator 53 in cell 4 to be accumulated and stored during the next cycle.

The subsequent machine cycles 6-13 for carrying out the backward substitution process proceed in a manner similar to that described above.

Figures 4F-4M detail the flow of data through the machine in solving subsequent equations. To briefly summarize, during each subsequent cycle another of the above equations becomes Xl1
1...Xo and the value of this component is stored in memory 57 of cell 1 along with the previously calculated value. Each newly calculated position is also sent sequentially to the right to each subsequent stage of the processor, as seen in Figure 4, and is accumulated as a partial product by an index according to the row of origin of the matrix element involved in the operation. Calculate QN'. For example, during cycle 6 (Figure 4F), three partial products are accumulated simultaneously for later addition in adder 65 as follows.

(1) The first partial product N-X1o is accumulated in the memory 57 of cell 2 (along with any other previous partial product sharing the same row index as element N). In this way, element N
Since the origin is from row 3 of the matrix, N−X
1oFiQ,'.

(2) The second partial product R-X1 is also Q under the exponent 3, 1
It is stored in the storage device 51 of the cell 3 as a.

(3) Simultaneous trap, since matrix element U is from row 2 of the matrix, the third partial product “Xl l is Q
A third partial product, V-X□, is accumulated along with other partial products representing the row superthickness shared with , t.

The first two partial products remain stored in storage 51 until the machine cycle during which counter 61 accesses all partial products associated with row 3 of the matrix. This occurs in cycle 9 (fourth engineering drawing). The row 30 partial product Q3' previously accumulated during that cycle
are applied to adder 65. During cycle 100, during cycle 11 the difference Q, -Q,' is expressed as a diagonal element O
Addition (N-X engineering.+R
-X1□) can be used. As a reminder, it can be seen that equation 3 above can be solved during cycle 11. The reason is only that in previous cycles the accumulation was done on the basis of the row index of the product of all the matrix elements of row 3 with the previously calculated value of the unknown vector.

Several important features of the methodology of the present invention will be apparent to those skilled in the art. The main feature is that the above-mentioned 12 simultaneous equations can be solved in 13 steps. This is because data dependence is removed as a result of the bandwidth-maximized form of the matrix used in the solution;
As described herein, this is only possible in combination with a unique pattern for storing and accessing matrix values on a multiprocessor. More specifically, if during a machine cycle the first few values of the unknown vector Operations are performed simultaneously in other processor cells to generate partial products (QN') using previously calculated components.

Proceeding in this manner, when partial products that depend on previously calculated unknown vector components are needed, those components can be used for calculations. In this way, the one-backward substitution method proceeds without delay until K are available so that these partial products can be utilized.

As mentioned above, the minimum bandwidth K of the band of zeros located between the main diagonal of the matrix and the upper right corner collection of nonzero elements υ is selected based on the specifications of the multiprocessor system. In general, K should equal the maximum number of cycles or steps required to gather and make available the data necessary to compute any given unknown vector component Xl. For the apparatus shown in FIG. 4, for example, during cycle 2
. 6 more to give
It requires steps or 6 cycles (3-8).

Therefore, a band@K of 6 is sufficient to give one space and time to fully utilize the particular contraction play shown in FIG. The particular bands required for other hardware devices will of course depend on the characteristics of that hardware device. Given a narrower bandwidth, some efficiency may be sacrificed depending on the composition of the data in the system of equations being manipulated.

As commonly understood in this technical field,
breadthdness (bkaend@dness) roughness is the connectivity of the finite element network (conn*ctivity)
reflect. Therefore, the system matrix (K)
It is well known in the art that it is well within the skill of those skilled in the art to convert A into a matrix [A] that maximizes the bandwidth of the upper triangle having a coarse structure as shown in FIG. 5B. It is about.

Minimizing the bandwidth as shown in Figure 5A involves renumbering the nodes of the finite element network such that the minimum difference between any two adjacent nodes is significant. However, maximizing the matrix in FIG. 5B is accomplished by renumbering the nodes of the system to maximize the difference between adjacent nodes. This can be done in accordance with the present invention by trial and error techniques, such as using random numbers to assign the first pass to the nodes, without requiring a lengthy period of time.

Because typical meshes are large, this technique typically removes 90% of the nonzero data from bands near the diagonal. Further non-zero data removal is performed by switching the numbers between the nodes of non-zero data inside and outside the band.

By continuing this trial and error technique a little longer, bands will be obtained that are wide enough to allow the decomposed matrix to operate efficiently in accordance with the techniques of the present invention. Furthermore, since a relatively small number of processors are used, the band size is generally smaller than the total number of finite elements in the system network.

As those skilled in the art will understand, by using the matrix A with the maximum bandwidth as in the present invention, the plural unknown components XN of the unknown field variable vector can be
It can be calculated without relying on other values of field variables. if,
Assuming that the calculations are performed sequentially in one processor (cell 1 in the above example), the other processors of the device (cell 1 in the above example)
Cells 2-4) calculate and calculate the partial products QN' that are later required to achieve the calculation of all components of XN (especially those that depend on previously calculated values) within the same time. Accumulation can be performed. Therefore, by increasing the number of XN calculations that can be performed without data dependence, time is given to the multiprocessor device to perform intermediate calculations without increasing idle time. All processors in the device are forced to wait for some part of the backward assignment method.
K so that meaningful calculations can be made during each machine cycle.

Although the matrix calculation method specifically mentioned here is a backward substitution method, a forward elimination method can also be performed in a similar manner. The forward elimination method can be implemented by one skilled in the art using the principles of the present invention. However, the forward elimination method converts the original system matrix into the upper triangle used in the backward substitution method.
Need to convert to lower triangle. Forward elimination methods are also described in the above documents.

Although the invention is described in terms of a small matrix and a small system of equations in order to briefly explain the invention, the advantages of the invention are We emphasize here again that the large and coarse associated system of equations follows more clearly from the operation on the IJ lux.

Although the invention has been described in terms of preferred embodiments including use in backward assignment operations performed on contracted computer architectures, the invention is not limited to the operation of any particular parallel multiprocessor matrix.

[Brief explanation of drawings]

FIG. 1 shows the basic contraction architecture used according to one method of the invention to process large, coarse matrices, and FIG. 2 shows how to process large, coarse matrices into parallel multiprocessor storage. FIG. 3 is a flowchart illustrating a method for mapping a particular matrix, e.g., a large, coarse, bandwidth-maximized matrix, to a parallel multi-gross processor array in accordance with one embodiment of the present invention. 4 shows backward substitution in the solution of a triangulated system of linear equations characterized in part by a bandwidth-maximized particular decomposition matrix stored in a parallel multiprocessor device according to FIG. Showing how to carry out the law, 1st 4A~
The 4M diagram is shown in FIG. 4 during a subsequent machine cycle performing a backward substitution method on the decomposition matrix shown in FIG. 4 with exemplary stored bandwidth maximized. FIG. 5A shows a prior art bandwidth-minimized decomposition matrix, and FIG. 5B shows a bandwidth-maximized decomposition structure of the present invention. The rough structure of the decomposed matrix is shown. 10... Host, 12... Inter 7 Eighth Units), 15... Cell, 51... Multiplier, 5
2...Divider, 53...Accumulator, 5B, 57.
58.59...Storage device, 6°...Subtractor, 6
1...Counter, 65...Adder. Sub-agent Masaki Yamakawa (#ν12 people) Engraving of drawings (;']'l't? No changes) Procedural amendment ζ Date 1. Indication of the case 1932 patent No. 22'360 2, name of Yun's eye r'+4I77 7 beggar elephants) λ method 3, person making amendment

Claims (2)

[Claims] (1) In order to analyze physical phenomena, unknown N vector Y
a first linear equation characterized by a large and coarse known system matrix of order N×N that relates R to a known N vector R by a vector representation KY=R, the unknown vector, and the known vector. generating a system having only nonzero elements on a main diagonal, a band of zero elements extending from said main diagonal by a preselected bandwidth, and no nonzero elements in said band; transforming said first system of linear equations into a second system of equations characterized by a decomposed triangular matrix such that all remaining nonzero elements are located in the outer part of said band; and said decomposition. loading the elements of the decomposed matrix into a memory associated with a plurality of processors such that the diagonal elements of the decomposed matrix are contained in a memory associated with one of the processors; solving the second system of equations by calculating values of a plurality of components of the unknown N vector independently of previously calculated values of any components of the N vector of A method for analyzing characteristic physical phenomena. (2) A process of dividing the solution region into a large but finite number of elements, and an approximate function determined by the value of the unknown field variable at the specified node existing on the boundary of the elements. and assembling the element functions into a first system of equations that describes the behavior of the solution domain, the first system of equations having the general matrix form KY=R; In the finite element analysis method of analyzing physical systems, where Y represents the unknown value of the field variable at a node in the system of equations, R represents the known boundary value, and K is a large and coarse matrix. an improved process for solving the system of equations by converting the system into a large, coarsely decomposed triangular matrix with only nonzero elements on the main diagonal; A finite element analysis method for analyzing physical phenomena, characterized in that a band of zero data is placed at a right-angled corner of a matrix, and a band of zero data is placed between the diagonal and the cluster.