JPH0752748B2 - 3D device simulation device - Google Patents

Description

Description: TECHNICAL FIELD The present invention relates to a three-dimensional device simulation apparatus for performing characteristic analysis of a semiconductor device on a computer without manufacturing an actual device, for example.

[Prior Art] In device simulation in which characteristics of semiconductor elements are analyzed on a computer without using actual elements, it is necessary to execute large-scale sparse matrix calculation at high speed.

Numerical simulation technology, which seeks to quantitatively elucidate phenomena, is becoming more important due to advances in computers. Techniques such as process simulations, device simulations, and circuit simulations have already been put into practical use when designing semiconductor integrated circuits, and there is a strong demand for shorter simulations in order to shorten the design period.

Device simulation ultimately results in a solution problem for large-scale matrix equations, and most of the simulation time is spent in this solution part. This tends to become stronger as the device structure becomes more complicated and three-dimensional simulation becomes indispensable, and there is a demand for speeding up the solution-solving portion.

[Problems to be Solved by the Invention] It is an object of the present invention to provide a simulation apparatus for a three-dimensional device that can achieve high speed using a systolic array. The systolic array is one of the VLSI-oriented architectures proposed by HTKung and has been studied from a logical and practical side. However, when processing a dense matrix proposed so far on a systolic array and performing a large-scale sparse matrix calculation such as device simulation, millions of arithmetic cells are required.

Therefore, the present invention is intended to provide a systolic array capable of effectively handling a large-scale matrix with a small number of arithmetic cells.

[Means for Solving the Problem] The three-dimensional device simulation apparatus according to the present invention has five types including a basic cell, a lower triangular element output cell, a diagonal element output cell, an upper triangular element output cell, and a time adjustment cell. Of hexagons, each element of the matrix is input from above, and each element of the decomposed matrix is output from below. A systolic array for incomplete LU decomposition and a plurality of multiply-accumulate operation cells are arranged in a one-dimensional array, and a vector element is specified from one end of the one-dimensional array of the plurality of operation cells at time intervals. Each element of each matrix is input to each operation cell at a specified time interval, and the product of each matrix element and each vector element is sequentially accumulated to obtain an element of an output vector. A sparse banded matrix systolic array for vector product which is adapted to output a plurality of product-sum operation cell is arranged one-dimensionally,
The elements of the vector are input from one end of the one-dimensional array of the plurality of operation cells at specified time intervals, and the elements of each matrix are input at the specified time intervals to each operation cell, And a systolic array for forward / backward substitution in which the product of each matrix element and each vector element is sequentially subtracted and output as the element of the output vector. The rick array, the sparse vector product systolic array, and the forward / backward substitution systolic array are connected by a data bus and a control bus, respectively, so that arithmetic operations can be executed based on commands from the main calculation unit. .

[Operation] According to this simulation device, the problem can be divided based on the characteristics of the device simulation problem, and the problems can be solved in parallel. In particular, sparse matrix vector products, forward / backward substitutions, vector sums, and inner products are repeatedly calculated in a loop, and therefore a large difference occurs in the calculation time depending on whether or not the speedup of this part is possible.
This simulation apparatus can provide a systolic array that can efficiently handle a large-scale problem with a small number of cells in order to speed up device simulation.

[Embodiment of the Invention] An embodiment of the present invention will be described below with reference to the drawings. FIG. 1 shows the configuration of the device simulation apparatus, which is an incomplete matrix decomposition (Imd = I) of a sparse band matrix.
ncomplete Matrix Decomposition) module 11, matrix-vector product (Mvm = Matrix-Vector Multiplication) module 12, and forward / backward (Fbs = Foward-Backwar)
d Substitution) module 13, and these modules 11 to 13 are connected to the host computer 15 via the global data bus 14. Where Imd module 11, Mvm module 12, and Fbs module 13
There are Imd controller 16 and Mvm controller 1 respectively
7 and an Fbs controller 18 are provided, and these controllers 16 to 18 are controlled by the host computer 15 via the global control bus 19.

Here, considering a semiconductor device simulation, this simulation seeks the electrical characteristics of the element based on information such as the impurity distribution of the semiconductor element and the lightness of the junction surface. , The electron density n, and the hole density p can be appropriately expressed by the Poisson equation, the current continuity equation for electrons and holes, and the like.

In device simulation, these three expressions are linearized and discretized, and then numerically solved by giving appropriate boundary conditions. Here, as the discretization method, there are a finite difference method, a finite element method, and the like, but it is effective to adopt the finite difference method in consideration of consistency with the systolic algorithm.

In the three-dimensional simulation, the attention point D shown in FIG. 2 and grid points (A, G, C, E, B,
Discretize the above three equations for a total of 7 points in (F) (7
Point difference approximation).

Serial numbers are assigned to grid points with grid coordinates (k, m, n) according to the following formula.

i = k + x (m−1) + xy (n−1) (1) (1 ≦ k ≦ x, 1 ≦ m ≦ y, 1 ≦ n ≦ z) where x, y, and z are discretization grids X-axis at,
It is the number of grid points on the Y axis and the Z axis.

At this time, the physical quantity Ψ _i such as the electric potential at the i-th point is expressed by the following equation.

a _i Ψ _i-xy + b _i Ψ _ix + c _i Ψ _i-1 + d _i Ψ _i + e _i Ψ _{i + 1} + f _i Ψ _{i + x} + g _i Ψ _{i + xy} = h _i (2) where a _{i to} g _i are coefficients, and h _i is a constant term.

When this is expressed by a matrix equation (AΨ = h), it becomes as shown in FIG. 3, and the coefficient matrix A (the number of dimensions N = xyz) is a regular, large-scale, sparse band matrix. For example, when each axis is divided into 40 lattice points, it is a simultaneous linear equation of 64,000 dimensions.

As the miniaturization of devices progresses and the device structure becomes more complex,
In three dimensions, the number of dimensions of the coefficient matrix in the simulation is from hundreds of thousands to millions, and it is indispensable to speed up the solution portion.

The solution method is mainly the iterative method of the conjugate gradient method (CG method) system, and the basic operations are incomplete decomposition of sparse band matrix, sparse band matrix vector product, forward / backward substitution, vector Is the dot product of and the sum of the vectors. FIG. 4 is a flow chart showing an overview of this solution.

Regarding the systolic array for sparse band matrix, the incomplete decomposition of the sparse band matrix will be described below.

The coefficient matrix A obtained by the device simulation may be a symmetric matrix or an asymmetric matrix. In either case, a systolic array that can perform incomplete decomposition is required. Now, if the coefficient matrix A is incompletely decomposed as follows, it becomes as follows.

A = LDU + R where L is a lower triangular matrix, D is a diagonal matrix, U is an upper triangular matrix, and R is an error matrix generated by incomplete decomposition.

The set (i, j) of the row number i and the column number j of the non-zero elements of A, L, D, and U in this equation are included in the sets S, ^SL , ^SD , and S ^U , respectively (1≤i , J ≦ N).

S = {(i, j) | j = i-xy, ix, i-1, i, i + 1, i + x, i + xy} S ^L = {(i, j) | j = i-xy, i-xy + 1, i -Xy + x, i-x, i-x + 1, i-1, i} S ^D = {(i, j) | j = i} S ^U = {(i, j) | j = i, i + 1, i + x-1 , i + 1, i + xy−x, i + xy−1, i + xy} In the following description, diagonal columns including non-zero elements designated by S of matrix A are set as (a, b, c) as seven N-dimensional vectors, respectively. , d, efg), and the row numbers of the matrix A are used as subscripts. Also, with respect to L, D, and U, diagonal columns including non-zero elements are used as vectors, and (a ^L , a ^L ′,
a ^L ″, b ^L , b ^L ′, c ^L , d ^L ), d ^D , (d ^U , e ^U , f ^U , f ^U , g ^U ″, g ^U ′,
g ^U ) and also the line number as a subscript.

The calculation formula for the incomplete LU decomposition (d ^D _i = d ^U _i = 1) is shown below (1 ≦ i ≦ N).

_{^{a L i = a i a L}} 'i = -a L i e U i-xy a L "i = -a L i f U i-xy -a L' f U 'i-xy + 1 b L i = b _i −a ^L _i g ^U ″ _i-xy b ^L ′ _i ＝ −a ^L ′ _i g ^U ″ _{i-xy + 1} −b ^Li e ^U _ix c ^L _i ＝ c _i −a ^L _i g ^U ′ _{i -xy} −b ^L _i f ^U ′ _ix d ^L _i ＝ d _i −a ^L ′ _i g ^u _i-xy −a ^L ′ _i g ^U ′ _{i-xy + 1} −a ^L ″ g ^U ″ _{i-xy + x} −b ^L _i f ^U _ix −b ^L ' _i f ^U ′ _{i-x + 1} −C ^L _i e ^U _i-1 e ^U _i = (e _i −a ^L ′ _i g ^U _{i-xy + 1} − b ^L ′ _i f ^U _{i-xy + 1} ) / d ^L _i f ^U ′ _i ＝ (− a ^L ″ _i g ^U ′ _{i-xy + x} −c ^L _i f ^U _i-1 ) / d ^L _i f ^U _i = (f _i −a ^L ″ _i g ^U _{i-xy + x} ) / d ^L _i g ^U ″ _i = (− b ^L _i g ^U _ix −b ^L ′ _i g ^U ′ _{i-x + 1} ) / d ^L _i g ^U ′ _i = −c ^L _i g ^U _i-1 / d ^L _i g ^U _i = g _i / d ^L _i (3) Summarizing the data dependence of this equation (3), each data Focusing on the regularity of the reference relation between the two, a systolic array that performs incomplete decomposition of the sparse matrix is obtained, as shown in FIG. In this figure, the circled arrows indicate the portions that require delay between adjacent cells, and the delay times are δ1, δ2, and δ3, respectively.

This systolic array includes a basic cell shown in FIG. 6A, a lower triangular element output cell shown in FIG. 6B, a diagonal element output cell shown in FIG. The upper triangular element output cell shown in d) and the time adjustment cell shown in (e) of the figure are arranged in a hexagon. Then, as shown in FIG. 5, each element of A is input from above the array, and each element of L, D, and U is output from below.

Here, the timing of each input data is as follows.
The first element of each of the seven vectors of matrix A, a1, b1, c1, d
1, e1, f1, g1 are time T = 1,5,7,8,10,14,22, respectively.
Input to each array, and the time interval for subsequent input is 3
Is.

At this time, the output d ^L _N is output at time T = 8 + 3 (N-1) + 7 + 1.

Of the five types of cells that make up the array, the time adjustment cell shown in (e) of FIG. 6 is a cell for delaying data without performing actual calculation. However, the data flow in the array is regular. Therefore, it is necessary to correctly calculate the equation (3). Further, in this systolic array, the following three types of decomposition can be performed by slightly changing the cell function.

I. Incomplete LU decomposition (d ^D _i = d ^U _i = 1) II. Incomplete LDL ^T decomposition (d ^D _i d ^U _i = 1, d ^L _i = d ^U _i ) III. Incomplete LDL ^T decomposition ( d ^L _i = d ^U _i = 1) The calculation formula for the incomplete LDL ^T decomposition (II, III) may be obtained by changing part of the formula (3) as follows.

II. Let a ^L _{i to} c ^L _i on the right side of the equation (3) be the following ° a ^L _i to ° c ^L _i
In addition, the expression without division by d ^L _i .

III. Let a ^L _{i to} c ^L _i on both sides of equation (3) be the following ° a ^L _i to ° c
To ^L _i, further replacing the d ^L _i to d ^D _i Equation.

However, ゜ a ^L _i ＝ a ^L _i d ^D _i-xy゜ a ^L ′ _i ＝ a ^L ′ _i d ^D _{i-xy + 1}゜ a ^L ″ _i ＝ a ^L ″ _i d ^D _{i-xy + x}゜b ^L _i = b ^L _i d ^D _ix ° b ^L ′ _i = b ^L ′ _i d ^D _{i-x + 1} ° c ^L _i = c ^L _i d ^D _i-1 By the way, this incomplete decomposition will be described later. -The original purpose is for use in back substitution, and the calculation is performed in a loop in the iterative method. Therefore, the sheet strip lower triangular element output cell in the array, and pairs diagonal elements devised to forward and backward substitution functions of the output cell, the data that can be used directly (e.g., d ^L-1 _i Ya゜a ^L _i ~ ° c ^L _i ) and output.

FIG. 7 shows a specific example of the configuration of the Imd module 11 that performs incomplete matrix decomposition in this way. That is, an input is given from the decomposed matrix (A) buffer memory 111 to the incomplete matrix decomposition array 112, and an output from this decomposition array 112 is supplied to the decomposition matrix (L, D, U) buffer memory 113. become.

Next, consider the sparse matrix-vector product. Each component q _i of the sparse band vector product (q = Ap) is expressed by the following equation.

q _i = Σ ^N _{j = 1} A _ij P _ij , (i, j) εS (4) The array is a one-dimensional array consisting of seven product-sum operation cells shown in FIG. 8, and this is shown in FIG. Configure as shown in. Then, each element q _i of the vector q is input from the cell at the left end with an initial value of 0, each element p _i of the vector p is input from the cell at the right end at a time interval of 2, and the calculation result q _i is expressed by the formula (4 ), Go to the right in the array. Where δ
_{1 to} δ ₁₂ is a delay time to be set between the cells, and according to the equation (4) and the time relationship of data movement in the array,
It is expressed by the following conditional expression.

Σ ⁶ _{i = 1} δ _i = Σ _{i = 7} δ _i = 2xy-6, Σ ⁵ _{i = 1} δ _i = 2x-4 + δ ₇ , Σ ¹⁰ _{i = 7} δ _i = δ ₁ + δ ₂ , Σ ⁶ _{i = 1} δ _i = δ ₇ δ ₈ , Σ ¹¹ _{i = 7} δ _i = 2x-4 + δ ₁ , Σ ³ _{i = 1} δ _i = Σ ⁹ _{1 = 7} δ _i , (5) This equation (5) and δ _i ≥0 And lead to the next result.

δ ₁ + δ ₁₂ = δ ₆ + δ ₇ = 2xy−2x−2, δ ₂ + δ ₁₁ = δ ₅ + δ ₈ = 2x−4, δ ₃ = δ ₄ = δ ₉ = δ ₁₀ = 0, δ ₁ + δ ₂ + δ ₅ + Δ ₆ = 2xy−6 (6) Further, when the time when a _{1 to} g ₁ is used for q ₁ is Ta _{1 to} Tg ₁ , Ta ₁ = 1 and Tb ₁ = 2 + δ ₁ and Tc ₁ = 3 + δ ₁ + δ ₂ , Te ₁ = 5 + Σ ⁴ _{i = 1} δ _i Tf ₁ = 6 + Σ ⁵ _{i = 1} δ _i Tg ₁ = 1 + 2xy (7), and the timing of each input data is expressed as a function of the delay time. It However, when p ₁ and q ₁ are input to the rightmost cell and the leftmost cell at time T = 1 respectively, the time at which p _N and q _N are output is T = 1 + 2xy + 2 (N-1) +1, and the delay time is It can be seen that the selection method of does not affect the total calculation period. Therefore, the delay time and the timing of the input data are not uniquely determined, and can be arbitrarily selected within the range satisfying the expressions (6) and (7). Although this is not a direct effect in the case of this sparse band-vector product, it comes to play an important role in the next forward / backward substitution.

Forward / backward substitution will be described. This part requires the most calculation time in the loop and is (LU) ⁻¹ r = t
Is calculated in two steps.

In forward substitution L ⁻¹ r = s, r is input to obtain s. Further, in backward substitution U ⁻¹ s = t, s is input to obtain t.

First, in forward substitution, each component s _i of s is calculated by the following equation.

s _i = (r _i −a ^L _i s _i-xy −a ^L ′ _i s _{i-xy + 1} −a ^L ″ _i s _{i-xy + x} −b ^L _i s _ix −b ^L ′ _i s _{i- x + 1} −c ^L _i s _i−1 ) / d ^L i (8) This equation (8) includes division by d ^L _i ,
Since the reciprocal has already been calculated at the stage of incomplete decomposition, multiplication is actually performed. Also, rather than a ^L _i ~c ^L _i is incomplete LDL ^T decomposition, ° a ^L _i ~ ° c ^L _i and d ^L-1 _i, or it is necessary to d ^D-1 _i, still incomplete decomposition At the stage of, the data is processed and only multiplication is performed.

The forward / backward substitution array is the same one-dimensional array as the sparse matrix-vector product array and can be constructed with the same number of cells. The function of cells other than the right end of this array is that the addition of α _out = α _in + β _in γ _in is replaced by subtraction _in the product-sum operation cell shown in FIG. The cell at the right end is a cell that outputs the result of forward / backward substitution, and its function is α _out = α _in γ _in β _out = α _out . Here, the time interval of the input data is 2, and the following result is obtained as the delay time between cells from the conditional expression of δ ₁ to δ ₁₂ as in the case of the sparse band matrix / vector product.

δ ₁ = δ ₁₂ = δ ₄ = δ ₉ = δ ₆ = δ ₇ = 0 δ ₂ + δ ₁₁ = 2x-4, δ ₅ + δ ₈ = 2x-6 δ ₃ + δ ₁₀ = 2xy-4x-2 ...... (8 ) Input data timings are Ta ^L ₁ = 1, Ta ^L ′ ₁ = 2 + δ ₁ , Ta ^L ″ ₁ = 3 + δ ₁ + δ ₂ , Tb ^L ₁ = 4 + Σ ³ _{i = 1} δ _i , Tb ^L ′ ₁ = 5 + Σ ⁴ _{i = 1} δ _i , Tc ^L ₁ = 6 + Σ ⁵ _{i = 1} δ _i , TdL ^L-1 = 7 + Σ ⁶ _{i = 1} δ _i (10), and r ₁ was input to the leftmost cell at time T = 1. Then, S _N is output to T = 7 + Σ ⁶ _{i = 1 1} δ _i +2 (N-1) +1. That is, unlike the case of the sparse band matrix vector product, the method of setting the delay time between cells is calculated. Time will be affected, so to minimize T, δ ₂ =
Let δ ₃ = δ ₅ = 0. The calculation time at this time is the shortest, and ends at time T = 2M + 6. By the way, the longest calculation time is when δ ₈ = δ ₁₀ = δ ₁₁ = 0, and T = 2N
The process ends at + 2xy-6.

In device simulation, xy is large and can reach tens of thousands, and since it is an iterative calculation, the choice of delay time is important.

On the other hand, the calculation formula for backward substitution is as follows (i = N ~
1).

t _i = s _i −g ^U _i t _{i + xy} −g ^U ′ t _{i + xy-1} −g ^U ″ _i t _{i + xy-x} −f ^U _i t _{i + x} −f ^U ′ t _{i + x -1} −e ^U _i t _{i + 1} (11) Therefore, the input data set (r _i , ...,
By replacing c ^L _i , d ^L-1 _i ) with (s _i , ..., e ^U _i , 1), the same array can be used and run under the same delay time settings. This is also the case when based on incomplete LDL ^T decomposition.

Here, FIG. 10 shows a specific configuration of the Mvm module 12 that executes the matrix-vector product, and FIG. 11 shows a configuration sequence of the Fbs module 13 that performs forward / backward substitution.

Further, FIG. 12 shows an example of cell configuration of the incomplete matrix factorization array shown in FIG. 5, where (A) is a basic cell, (B) is a lower triangular element output cell, (C) is a diagonal element output cell, (D) is an upper triangular element output cell,
(E) figure shows each time adjustment cell.

Further, FIG. 13A shows a concrete configuration example of the inner product operation cell constituting the product sum operation cell shown in FIG. 8, and FIG. 13B shows the sum operation cell. .

Considering means for speeding up this in device simulation, the problem can be divided based on the characteristics of the problem of the device simulation, and these can be solved in parallel. In particular, the sparse band-vector product,
Since forward / backward substitution, vector sum, and inner product are in a loop and are repeatedly calculated, there is a large difference in calculation time depending on whether or not this part can be accelerated. In the following description, a means for speeding up sparse matrix vector product and forward / backward substitution will be shown. However, the means for accelerating the forward / backward substitution is the same as the means for accelerating the sparse band matrix vector product, so only this will be described.

When the discretization is performed by the method described above, some of the neighboring grid points do not exist at the points located on the boundary of the simulation target region. With this in mind, the coefficient matrix A is divided into a square matrix in the xy dimension, and each is made a block matrix. Due to the above-mentioned property, this block matrix is either a diagonal matrix, a pentadiagonal matrix, or a zero matrix.

For each point (i, j), (1≤i, j≤z), let ° A ^ij be the xy-dimensional block matrix and ° p ^j , ° q ⁱ be the corresponding xy-dimensional vector. At this time, the sparse matrix-vector product (q = A
p) is divided by each block matrix and expressed by the following equation.

゜ q _i = Σ ^Z _{i = 1}゜ Ai _j゜ p _j (12) However, if s ^A = {(i, j) | j = i-1, i + 1}, then (i,
j) For s ^A゜ A ^ij = 0. Therefore, each degree q ⁱ of the equation (12) is specifically expressed as follows.

゜ q ¹ = A ¹¹゜ p ¹ + ゜ A ¹²゜ p ²゜ q ² = ゜ A ²¹゜ p ¹ + ゜ A ²²゜ P ² + ゜ A ²³゜ p ³゜ q ³ = ゜ A ³²゜ p ² + ° A ³³ ° p ³ + ° A ³⁴ ° p ⁴ 1 ° q ^Z = ° A ^ZZ-1 ° p ^z-1 + ° A ^ZZ ° p ^Z …… (13) where ° A ^ij-1 , ° A ^{ij + 1} is a diagonal matrix and ° A ^ij is a pentadiagonal matrix.

A high speed can be achieved by causing each array q q _{i in} equation (13) to be calculated by one array and operating z arrays in parallel. Furthermore, synchronizing with the output time of each data,
By calculating the inner product (q, p) of the vector in a pipeline, the overall speed can be further increased. This is called block parallel calculation, and the array that performs the processing for ° q _i is called [Array_i].

As shown in FIG. 14, this [Array_i] is realized by adding the inner product operation cell shown in FIG. 15 to the sparse band matrix vector product array. Further, for the pipeline processing of the sum, a cell in which the function of the product-sum operation cell shown in FIG. 7 is slightly modified is required.

This focuses on the fact that input data is input at a time interval of 2, and utilizes the time when there is no data in the data stream. The input data flow for γ _in is that each element of ° A ^ij flows at a time interval of 2, and the state of the time is (Red) and the state of the time between is (White). In FIG. 13, cells requiring a function change are distinguished by a + mark, and the function is as follows.

if (White) Register ← β _in + γ _in else if (Red) α _out = α _in + register × γ _in First, set the initial values of the accumulator and counter of the inner product operation cell to 0 and 1, respectively, and set the counter from xy When it is small, partial sums of inner products are sequentially accumulated in the accumulator in the block, and when the value of the counter becomes xy, that is, after the processing for the last element of xy-dimensional vector ° q ⁱ ° q ⁱ _xy is performed, The partial sum of inner products is [Array_i + 1]
Sent to γ _in .

The operation of this [Array_i] is performed as follows, for example. In equation (13), the first used in the calculation of ゜ q _i
p ^(i-1) is the sum of the vectors ゜ p ^{(i-1) a} and ゜ p ^{(i-1) b} ,
° p ^{(i-1) a,} ° p ^(i-1) the first element of ^b ° _{^{p (i-1) a 1}} , ° p ^{(i-1) b} ₁ is a product-sum in T = T _i, respectively Β _in the operation cell,
Input to γ _in as White column. The result of the addition is stored in a register in the cell and used in the calculation of ° q ⁱ ₁ at T = T _i +1. From T = T _i +2 to T = T _i + 2 + 2x, ゜ q ⁱ ₁ moves to the right while receiving the result of ゜ A ^ij゜ p ⁱ , and T = T ⁱ
At + 5 + 2x, ° q ⁱ ₁ ° p ⁱ ₁ is added to the accumulator, and at the same time, ° q ⁱ ₁ is output.

The whole operation starts with the processing of ° q ¹ at T = 1 in [Array_1], and then delays by one unit time [Array_2] [A
rray_3] and the sequential processing is started. As mentioned above, [Ar
ray_i] processing progresses, ゜ q ^Z _xy is T = 2xy + 2x + z
It is output from [Array_z] to +3, and at the same time, the result of the inner product (q, p) is also output.

In this way, in the block parallel calculation (T _bp = 2xy + 2x + z + 3) using multiple [Arrays], 8z cells are used, but T ₀ / T _bp compared to the case where it is not (To = 2xyz + 2xy)
= Z, which is about z times faster.

Considering the delay processing between cells, as a method of performing the delay processing between cells, a means using a shift register and a means using a memory can be considered. It is better to use a shift register for speeding up, but it becomes large as hardware. Using existing memory device
Means for realizing FIFO are more realistic.

The maximum delay time between cells is about 2xy (≧ 2 ¹¹ ), and the numerical data of the device simulation is 64 bits, so at least 16 kilobytes of memory is required. On the other hand, the block parallel calculation has a delay of about 2x, and if the delay memory is about 2 kilobytes, it can be used for simulations of several hundred thousand dimensions.

[Effects of the Invention] As described above, according to the simulation of the three-dimensional device according to the present invention, a fractional number of cells (for example, 64 for matrix decomposition, 7 for matrix-vector product, forward / backward) for speeding up the simulation. (7 for substitution) makes it possible to obtain a systolic array that can efficiently handle large-scale problems. Therefore, when N is the total number of lattice points in the target area of the three-dimensional simulation, the calculation time is proportional to N, and the speed can be greatly increased. In addition, a speed-up means that takes advantage of the characteristics of the coefficient matrix obtained by discretization was obtained. In this case, the number of required cells is N ^1/3 times, but the calculation time is N
It becomes proportional to ^2/3 .

Here, although the cells required for each array are of several types, the required arithmetic processing is similar. Therefore,
In order to realize this, it is possible to introduce all the functions and introduce a mechanism for switching the functions at the time of use to realize them as one cell. Figures 12 and 13 show the outline.

[Brief description of drawings]

FIG. 1 is a configuration diagram for explaining a device simulation apparatus according to an embodiment of the present invention, FIG. 2 is a diagram for explaining a 7-point difference approximation in a simulation, and FIG. 3 is a diagram for explaining a mode of a matrix equation. , FIG. 4 is a flow chart for explaining the solution of the basic operation, FIG. 5 is a diagram showing the configuration of the incomplete matrix factorization array, and (a) to (e) of FIG. 6 are the cells configuring the array, respectively. FIG. 7, FIG. 7 is a diagram for explaining the incomplete matrix factorization module, FIG. 8 is a diagram for explaining the function of the product sum operation cell, FIG. 9 is a diagram for showing a sparse band matrix / vector product array, and FIG. Is a diagram illustrating the configuration of a matrix-vector product module, FIG. 11 is a diagram illustrating the configuration of a forward / backward substitution module, and (A) to (E) of FIG. 12 are each the incomplete matrix factorization array. Example of specific configuration of calculation cell To figure FIG. 13 (A) and (B) is a diagram showing an example of the configuration of each product-sum operation cell 14
FIG. 15 is a diagram for explaining a block parallel calculation array, and FIG. 15 is a diagram showing inner product calculation cells added to this parallel calculation array. 11 …… Incomplete matrix decomposition module, 12 …… Matrix-vector product module, 13 …… Forward / backward substitution module, 14…
… Global data bus, 15… Host computer,
16-18 ... Controller, 19 ... Global control bus.

─────────────────────────────────────────────────── ─── Continuation of the front page (51) Int.Cl. ⁶ Identification number Office reference number FI technical indication location H01L 29/00 (56) Reference JP-A-1-129434 (JP, A)

Claims (1)

[Claims] 1. A hexagonal arrangement of five types of arithmetic cells consisting of a basic cell, a lower triangular element output cell, a diagonal element output cell, a lower triangular element output cell, and a time adjustment cell,
Each element of the matrix is input from above, and each element of the decomposed matrix is output from the bottom.A systolic array for incomplete LU decomposition of a sparse band matrix and multiple products The sum operation cells are arranged in a one-dimensional array, and the elements of the vector are input from one end of the one-dimensional array of the plurality of operation cells at a specified time interval. A systolic array for sparse matrix vector multiplication, in which elements are input at specified time intervals and the product of each matrix element and each vector element is sequentially accumulated and output as the element of the output vector, and Multiply-accumulate operation cells are arranged in a one-dimensional array, and the elements of the vector are input from one end of the one-dimensional array of the plurality of operation cells at a specified time interval. A systolic array for forward / backward substitution in which each element of each matrix is input at a specified time interval and the product of each matrix element and each vector element is sequentially subtracted and output as the element of the output vector. Wherein the systolic array for incomplete LU decomposition of the sparse band matrix, the systolic array for sparse band vector product, and the systolic array for forward / backward substitution are coupled by a data bus and a control bus, respectively, and a main calculation is performed. A simulation device for a three-dimensional device, characterized in that a calculation operation can be executed based on a command from the section.