On bandsolver using skyline storage

CO45-7949/92 165.00 + 0.W C 1992 Pergamon Press Ltd

Cornpurers & Srrurturcs Vol. 44. No. 6. pp. 1187-l 196. 1992 Printed in Creak Britain.

ON BANDSOLVER

USING SKYLINE S. H.

STORAGE

Lo

Department of Civil and Structural Engineering, University of Hong Kong, Pokfulam Road, Hong Kong (Received 9 .iu!y 1991) simple out-of-core bandsolver using skyline storage is presented. Without dividing the equations into blocks, the efficiency of the working vector can be greatly improved by a dynamic

Abstract-A

storage scheme to hold as many equations as possible at any one time. Compared to the blocked-skyline method, the total number of data transfers from core to disk and vice versa is reduced and the amount of active coefficients being processed in a single data transfer is maximized. The matrix decomposition without using direct-access files is first presented to explain how triangulation is done on a one-dimensional working array using skyline profile storage. Its extension to allow for out-of-core storage is then discussed in detail. Computer codes for both versions with and without using out-of-core storage are given in the Appendix.

1. INTROD~~ION

The solution of the system Kx = f is the most important step in finite element analysis. The number of unknowns n is directly proportional to the number of nodes and the number of degrees of freedom per node. The frequency and range of application of the finite element method is limited only by the number of simultaneous linear equations that can be solved economically in today’s computer. Methods of solutions are generally divided into two broad classes [l]: (i) direct methods also called Gaussian elimination methods, and (ii) iterative methods of which the Gauss-Seidel variation is the most popular. During the earlier years of finite element analysis, iterative techniques were thought to be more effective and easier to program than the direct elimination method. However, further research into algorithmic solutions and computer programming resulted in a near-unanimous choice of elimination methods because of their near-optimum accuracy and speed for general systems. Gaussian elimination is widely used as a direct solution procedure. The most common variation of the Gaussian elimination techniques divide the operations into two steps: (i) triangulation of the matrix coefficients-a systematic series of linear transformations is applied to the system of equations in order to reduce it to a triangular system and (ii) for an upper triangular system, the last equation contains only one unknown and is readily solved. The previous equation contains two unknowns, one of which has just been calculated. By successive substitutions, the remaining unknowns can be found. The process is repeated for each line of the system in sequence from bottom to top. Gaussian elimination can be refo~ulated into a two-phase process that does not require to modify

K and f simultaneously. Such a process is referred to as matrix decomposition; matrix factorization or triangulation is widely used in finite element analysis. In studying the matrix decomposition algorithms developed, it is apparent that the amount of stiffness matrix coefficient data which is actively being processed at any time is a very small percentage of the overall coefficient data vector. When decomposing coefficients of a given stiffness equation, only those columns having coefficients which intersect the particular equation column are needed [2]. As a result, elimination methods for matrices arising from finite element analyses can be classified as: (i) methods which apply when all equations can be contained in core storage, (ii) methods which apply if the active coefficients can be contained in core, and (iii) methods which apply when there is insufficient core storage to contain the active coefficients, In this paper, matrix reduction methods belonging to categories (i) and (ii) will be discussed. The stiffness coefficients are stored in a linear array using the efficient skyline profile storage scheme. In the skyline profile storage method, the matrix coefficients are partitioned in columns and are stored column by column. The coefficients stored in a linear array are addressed indirectly through the use of a pointer vector, which indicates the position of the last element of each matrix column in the onedimensional storage form. 2. MATRIX DECOMPOSITION 2.1. Lo decomposition By direct multiplication, it can be easily verified that matrix K can be decomposed into a lower triangular matrix L with diagonal element equal to unity and a upper triangular matrix fi

1187

S.

1188

H. Lo

such that LO= K, or 1

0

l

L,,

LX, LX, l .. .. . L,,

-G, -+ L,,,- 1 1

where 2.4. Cholesky decomposition

i-l

iTi,=Kii-- 2 L,Um,, j=i,n, m-l

i=l,n

(1)

j=l,n.

(2)

If K is positive-definite, all the diagonal elements D, of diagonal matrix D are greater than zero. Let C=[,/&,& ,..., @“J, then C2=D K = LDLT=

i=j$l,n;

where D is a diagonal matrix TO,, D,, . . , D,, J and iJ is an upper triangular matrix whose diagonal elements are equal to 1, then U, = fJJD,,

j = i,n;

i = 1,n.

In a finite element analysis, the resulting stiffness matrix K is always symmetric. By Crout decomposition, the matrix K can be decomposed into a product of lower triangular matrix L, diagonal matrix D and upper triangular matrix U, in which L and U are transposes of each other. As a result, either L or U needs to be considered in the matrix decomposition process. If we concentrate on the upper half of symmetric matrix K, the matrix decomposition can be conveniently considered as a column-by-column reduction process. For instance, column j will then consist of all the elements from top of the matrix down to the diagonal, Klj, rc,,, . . . , K,_ I,i and diagonal element Di as shown in Fig. 1. With L, = U,, the equation for the reduction of column j is

(3)

2.3. Crout decomposition By matrix decomposition K = LDU

in the last section, KT= UTDLT

and

(4)

matrix K, K’= K, we have

for symmetrica

UrDLr=

LLW.

(5)

Since matrix decomposition is unique [3], UT= L and LT = U, symmetric matrix K can be expressed as K = LDL’=

i-l

fTi,= Kij-

UrDU.

c

U,,iu,,,,,

i = l,j - 1.

m= 1

(6)

,column

under

reduction _

. ' '1t

U

22

-.-

uzi

:

J

ui-,

,i

i

i.

11

L

i2'

”

Li,i-l

K K

4j

R--.

3-O”

(7)

2.5. Reduction column by column

u= DU,

Di = ui,,

L,.LJ,

where lower triangular Cholesky matrix L,. = LC.

2.2. LDU decomposition Let

LCCLr=(LC)(LC)*=

12

...

Ktj

K1~

22

...

'2j

K2r

Kii

.

=

'

'i-1

j ,

Uij

0 K

Fig. I. Reduction of upper triangular matrix column by column.

nn

(8)

A

bandsolver using skyline storage

This equation for determining oii can also be interpreted as the result of equating the dot product of column vectors i and j to coefficient Ki,. By putting i =j, the equation for the reduction of diagonal element D, is given by 1-l

c

D, = K,,-

Urn,~mm,’

(9)

m=l

The reduction for column j will be completed when U,,, are obtained from fl,,,, by dividing them with the diagonal elements D, such that Umi= u,,,,/Dmr

m = l,,j - 1.

(10)

3. SKYLINE

solution

Kx = f+LDCJ)x

=f

vector x in the is to be obtained

equation in three

steps. 2.6.1. Forward substitution. Let z = (DU)x we can write

z can

easily

1

1 L,,z,,

,=I

0

1

ZI

fi

22

.fi

z,

f;

z,,

f;,

by forward

L,

1-I

z,=f;-

>&I.

then

L2,

be solved

,&I= [K,,, h,

There are 17 off-diagonal elements within the skyline profile, and a linear array of dimension 17 is required

Lz=f

for which substitution

SCHEME

remain zero in the reduction process; however, zero terms within the skyline must be stored as these terms are not invariants in the matrix decomposition process. Considering the symmetric 8 x 8 matrix shown in Fig. 2, number of equations n equals to 8, so a linear array D of dimension 8 is required to store the diagonal elements K,, , K12, , KE8 D=[D,,D?,..

for

STORAGE

The most efficient strategy for the storage of nonzero elements of a upper triangular matrix is to store the matrix coefficients column by column of variable lengths marked off by a pointer vector on a linear array. The profile of the skyline is the envelope of the columns of variable heights. Zero terms outside the skyline envelope need not be stored as they will

2.6. Solving for x The

1189

i=2,n,

.

L,i - I 1

where

L, = 2.6.2.

L n.n I I

L,

u,,

Dividing by diagonal

elements.

Writing 4

y = Uz, we have

4

Dy=z

0

Yl

ZI

Y2

22

from which yg= z,/Di,

i = 1, n. 0

Y,,

2.6.3. Backward substitution. Finally, when vector y is known, vector x can be solved from the equation (Ix = y by backward substitution, such that

xi= yj-

i ,=ICI

u,x,,

1

...

U,,

1

i =n - 1,l.

0 1

S. H. Lo

1190

to hold all these elements as illustrated in Fig. 2, with A =[A,rA~,-..fA,71=I~~t,~*~,.r.,jY781In addition, an auxiliary pointer P of dimension n = 8 is needed to marked off each column on the linear array A, such that P, is the number of ok-diagonal element up to column i. Hence, the total memory required to store symmetric matrix K equals to

3.3. Sequence of column reduction To optimize storage space and arithmetic operations, the same storage will be used to store the matrix coefficients before and after column reduction. The diagonal and off-diagonal elements are reduced in the following sequence

Dim(D) + Dim(A) + Dim(P) = n + n + N

where N is the number of off-diagonal elements under the skyline. 3.1. Address calculation With the aid of the pointer vector P, it would be easy to calculate the corresponding address of coefficients I$ on the linear vector A. Since Pj corresponds to the address of the last element Kj_ ,j of column j, the address of element Kti is given by

where

i = h,,j - 1

W: = 0,/D,*,

(16)

where r>,, yii are initial values of the matrix coefficients, U, are intermediate values, and D:, .!J$ are reduced values. 3.4. Column reduction using linear array Initial values W, which are in fact the coefficients of the upper triangular part of stiffness matrix K, within the skyline profile, are stored in a linear array A. Their corresponding address can be calculated by

wheref=P,-i+landJ=P,--j+l.Intermsof J=P,-j+

D, the equations

(12) elements of linear arrays A and

1.

for column reduction take the form

3.2. Column red~c~~~~~ii~in the skyline prc#ie

i-i

Since zeros outside the skyline profile are ignored,

A;+i=A,,,-~~*A:,,~A;+,,

in the calculation of the sum i=h,, where the number of arithmetic operations can be reduced by adjusting the lower summation limit form. Instead of beginning the summation from m = 1, m starts from a value h given by

i= Pi--i+

h,=P,_,--P,+-i,

hi=Pi_,--Pi+j.

(141

In the reduction of column j, we now have

m=h

U,=&,iDi>

i=h,, j-l

D, = K, -

j-l

c urn, ITi,

m = b,

i j= 1,n.

J

(15)

1 and

J=P,-j+

j-l

(13)

where hi and /I, are the heights measured from the first row to the first non-zero elements in column i and column j, respectively (Fig. 2). Heights hi and hi can be calculated from pointer vector P

with

h = max(h,, h,),

for

j=l,n,

iz,= P,_ , - P, + i,

hj=Pi-,-Pi+j

D,T= Dih = max&, h,),

j-l

Af,i=A-,+ii~:,

1

1

C $+m/D:

j=

1,n.

1”= h,

i=h,,j-1

4. USlNG OUT-OF-CORE

i

(18)

ARMORY

When we are not working on a virtual memory operating system, and the size of stiffness matrix K is so large that the amount of coefficients within the skyline profile exceeds the available central core memory, the linear array A must be stored column by column on a mass storage device such as a disk drive. The number of columns (equations) to be read from disk to core depends on the size of the working vector W whose minimum length is to be determined in advance prior to the matrix reduction.

1191

A bandsolver using skyline storage

.

l,J

.

=

1

2

3

D;

A;

4

5

6

a

7

t

hi

Df

A;

.-

2

r

reduction

A-

A;

unde

Cal umn

I1

l

D:

pro file

A;

A;

AY

Skyline

D4

A;

A:0

Reduction

DS

Aft

*13/

/

element *15

i

14

=A

of

: 14

A* i 10 12

-

A* i 11 13 D7

Al7 -* Da

P=

0

1

2

5

7

11

14

17

h=

o

1

2

1

3

7

L

5

o

0

0

2

3

6

a

10

1

4

6

6

7

a

a

a

1,J

= c=

Fig. 2. Skyline profile storage scheme--column 4.1.

Length of working vector W

Vector W provides a working place for the reduction of the elements in vector A. As many equations as possible will be read from disk to vector W for reduction. Having reduced all the coefficients held in vector W, a certain number of equations may leave the working vector W to make way for the reduction of other equations. Columns leaving and entering the working area are determined by a control vector C which indicates the columns no longer required in the subsequent reduction. The elements of vector C are Ci=~y{j;hj
(19)

The significance of C, is that column i can leave the working area if column C, has entered the working area and has been reduced. Hence, working vector W must have at least a length that can accommodate columns from i to C, for any i = 1, n. The number of elements contained in columns i to C, is given by P, -Pi-,.

(20)

Thus the minimum length for working vector W is fmin=~~~{pC,-pi-Il~

(21)

4.2. Reading data from direct -access files The number of columns and amount of data to be read from the direct-access file is controlled by the

* Al4

=

A,4/D

6

7 under reduction

currently available space in the working vector W. Let L 2 C,,, be the length of the working vector W, and L, be the space occupied by the equations (coefficients) remaining in the working area, then the space available L, will be L, = L - L,

(22)

Suppose columns j,, to j, remain in the working area, and j, th column is the last column reduced, then data from disk up to column j, will be read in such that L, > P,z - P II

(23)

Now on working vector W, with the first element of W corresponding to the (P,,_, + l)th element of vector A, columns j, + 1 to j, are reduced according to eqns (18). This is always possible as columns are not allowed to leave the working area if they are still needed in the subsequent reduction of the forthcoming equations (columns). To determine which columns can leave the working area after the reduction down to column j, the Ci values of the columns are checked. Those columns with C, values less than or equal to j, can now leave the high-speed core; these reduced coefficients are transferred back to the direct-access file for later forward and backward substitutions in solving for unknown vector x Column leaving = {i = jO,jr ; C,
(24)

S. H. Lo

1192

Fig. 3. Finite element analysis of cantilever beam. or up to and include column

arguments needed in subroutine REDUCI, routine REDUC2 takes additional arguments:

j = maxji = j,, j,; C, < j2].

(25)

The cycle can be repeated by putting new values of j, and j, to .J$=j+i

and

j:=j,.

(26)

5. PROGRAMS

FOR MATRIX DECOMPOSITION AND SOLVING FOR x

The Fortran subroutines for the matrix triangulation, forward and backward substitutions with and without using direct-access files are listed in the Appendix. 5.1. Version l-without

using out-of-core

memory

SU3ROUTINE REDUCl (N, IP, D, A). N is the number of equations (columns), vectors D and A contains diagonal and off-diagonal elements of the entire stiffness matrix K. IP is the pointer vector indicating the address of the last element of each column. N and vector IP will not be altered in the reduction process; and the matrix coefficients after reduction are stored in the same memory locations as occupied by the initial values in vectors D and A. SUBROUTINE SOLVE1 (N, ZP, D, A, B). B is the ~ght-hand vector to be solved, the solution will again be put back to vector B. 5.2. Version 2-using

out-of-core

memory

SUBROUTINE REDUC2 (N, IP, IC, D, A, LAW, Nl, N2, IREC). Two direct-access files will be used to store vector A before and after matrix reduction. Vector D of diagonal elements is however always maintained in the core memory. Apart from the Table

M 5 IO 15 20 25 30 35

N

sub-

LAW = Allowable length of working vector W. N 1, N2 = channel numbers of direct-access files where A is stored before and after decomposition. iREC(I = 1, N) = pointer vector indicating the beginning record of each column in direct-access file Nl. ZC(Z = 1, N) = vector controlling the reading and writing (entering and leaving) of matrix coefficients from and to direct-access files N I and N2. Vector IC is calculated in the subroutine REDUCI itself from vector IP. SUBROUTINE SOLVE2 (N. IP, D, A, B, N2, Using SOLVE2, reduced coefficients stored in direct-access file N2 can be called upon for the solution of any right-hand vector B. IREC).

6. EXAMPLES

The finite element analysis of the bending of a beam is quoted as an example to study the performance of the bandsolver with and without using outof-core memory. The calculation was done in a micro-PC-AT with 80286 processor running at 8 MHz, and the speed of the machine is about 675% of that of a normal IBM-PC. The dimensions, boundary conditions and loading of the beam structure are shown in Fig. 3. The results of the analysis for the beam dividing progressively into M x N eight-node quadrilateral elements are summarized in Table 1, in which NE = number of elements, NN = number of node points, NEQ = number of equations, NA = number of off-diagonal terms, CPU time = time required for matrix decomposition, S = deflection at the tip

1. Finite element analysis of cantilever

ME

NN

NEQ

I

5

28

56

2

20

85

170

3 4 5 6 7

45 80 125 180 24s

172 289 436 613 820

344 578 872 1226 1640

NA 540

2,589 7,048 14,817 26,196 43,885 66,984

beam CPU time W4

6 (mm)

1.8 11.4 38.9 110.2 228.0 417.6 717.4

42.21 42.69 42.77 42.79 42.80 42.81 42.82

(value from the elementary strength of material including shear effects = 42.92 mm for Young’s modulus E = 1200 N/mm* and shear modulus G = 480 N/mm’). From (It4, N) equal to (20,4) onwards, direct access-files had to be used in the matrix reduction process, and the CPU time required was much more than the last three cases in which no backing storage was needed.

by the optimum

REFERENCES

Editeur, Paris (1984).

A simple out-of-core bandsolver using skyline profile storage is presented. The efficiency of the is augmented

files proves to be a little bit more complicate, nevertheless, the algorithm can still be coded by about 100 Fortran statements. Since the computer coding is very short, they are most suitable to be used for structural analysis programs implemented in microcomputers.

1. G. Dhatta and G. Touzot, Une PrPsentation de la mkthode des elements jinis, 2nd kin. Maloine S.A.

7. CONCLUSIONS AND DISCUSSIONS

algorithm

1193

bandsolver using skyline storage

A

2. R. L. Taylor,E. L. Wilson and S. J. Sackett,Direct

solution of equations by frontal and variable band, active column methods. In Nonlinear Finite Analysis in

and dynamic

use of the working vector so that a maximum number of matrix coefficients can be processed at each cycle of data transfer. The computer subroutines corresponding to the version without using out-of-core memory are very simpie that only 40 Fortran statements are required. The version using direct-access

Structural

Mechanics

(Edited by K. J. Bathe and

E. Stein). Springer, New York (1981). 3. T. J. R. Hughes, The Finite Element Method. Linear Static and Dynamic Finite Element Analysis. Prentice-

Hail (1987). 4. 0. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, 4th Edn. Vol. 1, Basic Formulation and Linear Problems. McGraw-Hill (1989).

APPEND%X

SUBROUTINE REDUCl (N,IP,D,A) IMPLICIT DOUBLE PRECISION (A-Ii,O-2) DIMENSION IP(*),D(*),A(*) C C

FUNCTION

C

: CROUT DECOMPOSITION A = LDU REDUCTION COLUMN BY COLUMN

C

DO 11 K=2,N KY_=K-1 LX=IP(K)-Xl KH=IP(Kl)-LK+l S=D(K)-A(LK+KH)**2/D(KH) C

C

DO 22 J=KH+l,Kl Jl=J-1 LJ=IP(J)-Jl JH=IP(Jl)-IJ+l IF (KH.GT.JH) JH=KH T=A(LK+J) DO 33 M=JH,Jl 33 T=T-A(IJ+M)*A(LK+M) A(LK+J)=T 22 S=S-T*T/D(J)

44

DO 44 J=KH,K3. L=LKiJ

A(L)=A(L)/D(J)

11 D(K)=S RETURN END

SUBROUTINE SOLVE1 (N,IP,D,A,B) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IP(*),D(*f,A(*),B(*) C C C C C

FUNCTION 1. FORWARD

: SOLVE

FOR

SUBSTITUTION

DO 11 J=2,N Jl=J-1 LJ=IP(J)-Jl JH=IP(Jl)-Wil T=B(J) DO 22 M=JH,Jl 22 T=T-A(LJ+M)*B(M) 11 B(J)=T

X, LDUx

= b

: Lz = b

where

U = transpose

of L

S. H. Lo

II94 C C C

2. DIVIDING

BY DIAGONAL

ELEMENTS

: Dy = 2

DO 33 K=l,N 33 B(K)=B(K)/D(K) C C C

BACKWARD

: UX = y

SUBSTITUTION

DO 44 K=N,2,-1 Kl=K-1 LK=IP(K)-Kl KH=IP(Kl)-LK+l T=B(K) DO 44 J=KH,Kl 44 B(Jf=B(J)-T*A(LK+J) RETURN END SUBROUTINE DIMENSION C C C C C C c C C c C C C C C C C C C C C C C C C

FUNCTION

REDUC2 (N,IP,IC,D,A,LAW,Nl,N2,IREC) IP(*),D(*),A(*),IC(*),IREC(*) : TO REDUCE A SYMMETRIC MATRIX TO LOWER AND UPPER TRIANGULAR MATRICES. OUT OF CORE BAND SOLVER CROUT DECOMPOSITION A=LDU, WHERE U=TR.ANSPOSE OF (SKY-LINE STORAGE SCHEME IS USED)

L

INPUT

: N = NUMBER OF EQUATIONS IF'(K) = ~B~R OF OFF-DIAGONAL ELEMENTS UPTO COLUMN K D(K).= Kth DIAGONAL ELEMENT A(L) = Lth OFF-DIAGONAL ELEMENT Ni I CHANNEL NUMBER OF THE DIRECT-ACCESS FILE WHERE A IS STORE BEFORE DECOMPOSITION N2 = CHANNEL NUMBER OF THE DIRECT-ACCESS FILE WHERE A IS PUT AFTER DECOMPOSITION LAW = MAXIMUM LENGTH OF WORKING VECTOR A ALIOWED IREC(K) 5 READ COLUMN K FROM THIS RECORD IN D. A. FILES

OUTPUT

: D(*),A(*)

WORKING

ARRAY

CALCULATE

: IC(K),

VECTOR

COLUMN K CAN LEAVE THE HIGH AS COLUMN IC(K) IS REDUCED

SPEEN

CORE

AS SOON

IC(N)

DO 22 K=3,N Kl=K-1 KHl=IP(Kl)-IP(K)+K+l DO 22 J=KHI,K 22 IC(J)=K c C C

ENSURE

LENGTH

OF WORKING

ARRAY

A IS NOT

LMAx=o DO 11 K=2,N IF (IC(K).LE.K) GOT0 11 L=IP(IC(K))-IP(K-1) IF (L.GT.LMAX) LMAX-L 11 CONTINUE IF (IMAX.GT.LAW) THEN WRITE (2,30) LMAX,IAW 30 FORMAT (//,/' ERROR... LENGTH OF WORKING I GREATER THAN THE EXIT STOP END IF

GREATER

THE

VECTOR A REQUIRED LENGTH ALLOWED

C LST=O KS=1 KT=l IP(N+1)=1000000 IC(N+1)=1000000 REWIND Nl C C C

Read

columns

KQ+l

2 KQ=KT LT=IP(KT)+LAW-LST 1 KT=KT+1 IF (IP(KT).LE.LT) KT=KT-1 DO 44 K=KQ+l,KT

to

KT

GOT0

from

1

direct

ALLOWABLE

access

file

Nl

LENGTH

=',I7,/ =',17)

IA

A bandsolver using skyline storage L=IP(K)-IP(K-1) IR=IREC (K) READ (Nl,REC=IR) 44 C C C

(A(I),I=LST+l,LST+L)

LST=LST+L CONTINUE REDUCING

COLUMNS

KQ+l

TO KT

DO 111 K=KQ+l,KT Kl=K-1 LK=IP(K)-Kl KH=IP(Kl)-LK+l LKl=LK-LP S=D(K)-A(~l+~)**Z/D(~) DO 222 J=KH+l,Kl Jl=J-1 LJ=IP(J)-Jl JH=IP(Jl)-LJ+l LJl=LJ-LP IF (KH.GT.JH) JH=KH T=A(LKl+J) DO 333 M=JH.Jl 333 T=T-A(LJl+Mj*A(LKlSM) A(LKl+J)=T 222 S=S-T*T/D(J) DO 444 J=KH,Kl 444 A(LKl+J)=A(LKl~J)/D(J) 111 D(K)=S C

IF (KT.EQ.N) THEN J=O DO 66 K=KS+l,N L=IP(K)-IP(K-1) WRITE (N2) (A(I),I-J+l,J+L) J=J+L 66 CONTINUE RETURN ENDIF C C C

WRITE

COLUMNS

KR+l

KR=KS 3 KS=KS+l IF (IC(KS).LE.KT)

TO KS TO

GOT0

KS=KS-1

DIRECT

ACCESS

FILE

N2

3

LP=IP(KS) LRS=LP-IP(KR) LST=IP(KT)-LP J=O DO 55 K=KR+l,KS L=IP(K)-IP(K-1) WRITE (N2) (A(I),I=Yfl,JSL) J=J+L 55 CONTINUE DO 33 I=l,LST 33 A(I)=A(LRS+I) GOT0 2 END

C C C

SUBROUTINE DIMENSION

SOLVE2 (N,IP,D,A,B,N2,IREC) IP(*),D(*),A(*),B(*),IREC(*)

SOLVE

LDUX

FOR

= b,

USING

WRITE (*I*) ' 1. FORWARD REWIND N2 DO 11 J=Z,N Jl=J-1 L-IP(J)-IP(J1) READ (N2) (A(I),I=l,L) JHl=Jl-L T=B(J) DO 22 M=l,L 22 T=T-A(M)*B(JHl+M) 11 B(J)=T

OUT-OF-CORE SUBSTITUTION

STORAGE : LZ = b'

C DIVIDING

BY DIAGONAL

ELEMENTS

WRITE (*,*) ' 3. BACKWARD DO 44 K=N,2,-1

SUBSTITUTION

: UX

WRITE (*,*) Do 33 K=l,N

' 2.

: Dy=

33 B(K)=B(K)/D(K) C = Y'

2'

1196

S. Kl=K-1 L=IP(K)-IP(K1) IR=IREC(K) READ (NZ,REC=IR) (A(I),I=l,L) KHl=Kl-L T=B(K) DO 44 J=l,L 44 B(KHl+J)=B(KHl+J)-T*A(J) RETURN END

H. Lo