AO-I477 N AGRITM NOASMBfO STFN S M NTIE INTO A /4

COMPCE D T SR . U A ORN AU BERKEE DEPT

01m /2 hE
0FCII ENGIER BN OUROMODE [ E 84

432N004-6
mSIFEDNEL-R
EN
homhmhI
 12.8
-21

11111.25
111
ul1111.6

MICROCOPY RESOLUTION TESTCHART

NA11ONAL BUJREAUOf STANDARDS 196 A
 gT OPP
CR 84.032

NAVAL CIVIL ENGINEERING LABORATORY

401EEPort Hueneme, California

Sponsored by
NAVAL FACILITIES ENGINEERING COMMAND)

AN ALGORITHM FOR ASSEMBLY OF STIFFNESS

MATRICES INTO A COMPACTED DATA STRUCTURE

~~ September 1984ji

LJJ

DTIG -"
.I
An Investigation Conducted by:
DEPARTMENT 01- CIVIL ENGINEI*.RIN; ~ELECTE
UNIVERSITY OF CALIFORNIA

BerkclcN, CA 95663OC2

N00014-76-C-OO1 3
N62583/83 M T265

Approved for public release; distribution in unlimited.

84 10 23 008 1
 t t
E~ I. S

z 0N CI w-O l4"
0-0 00 -. 0N
ri 0-~

E E ar

I..I rp 'U.I... II IN

S,.-
>I it
I-

N '.C

ILL 'L irO~ i biI Li

SI~U~jI~g& ~lI& l~ll6 8 19 S P LI k
I II I l ii I i ii I~aa

00

E. C 5r
 Unclassified
SC.JRITY CLASIICATION OF TH~IS PAGE (WP,.n noW. F..,...d)

REPOT DCUMNTATON
AGEREAD
REPOT DCUMNTATON
AGEBEFORE
INSTRUCTIONS
COMPLETING FORM
IREPORT NUMBER 0 3 IE CATALOG NUMBER

CR-84 .032
AILE (.d S,b-1I.)
-0"''
U1
S
7 YP
WO RPORT A PERIOD COV'EREO

An Algorithm for Assembly of Stiffness Final

Matrices into a Compacted Data -May 1980 - Dec 1983
Strucure PERFORMING ORG REPORT N~iMERI

1AUTYOR1(*, B CONTRACT OR GRANT NUNBER(AI

Bahram Nour-Omid and N00l4-76-C-0013

Robert L. Taylor N62583/83 M T265
9 PERFORMINGORGANIZATION NAME AND ADDRESS 10 PROGRAM ELEMENT paciEC T. TASK
Department of Civil Engineering AREA a WORK UNIT NUMBDERis

University of California YR023.03.01.003

Berkeley, CA 95663 ___________F
It ONTROLLING OFFICE NAME AND ADDRESS IZ REPORT DATE
Naval Civil Engineering Laboratory -September 1984
Port Hueneme, CA 93043 16NME
O AE

14MONITORING AGENCY NAME # ADDRESS,I dIlefr.I I,,, Co-,II,nll OfiIr.J IS SECURITY CLASS W~ thA.I.Po.I)

Naval Facilities Engineering Command Unclassified

200 Stovall Street IS&OECLASSIFICATION DOWNGRADING

Alexandria, VA 22332 SCHEDULE ________

IA DISTRIBUTION STATEMENT (,,l .11,. R.Pn1,l

Approved for public release; distribution is unlimited.

17 DISTRIBUTION STATEMENT (.1 111. *.st-erent-,d in Block 20. 1 d,ffe-.o I,-. R*Porf)

1B MENI ANY NOTF

EIIi' S

IS o&E BORDS ICo--~lO~IIO,,d I d.,,I IfrlIf

, block 'IO b.SOI

finite element analysis, structural engineering

20 ABSTRAC T IC-flIsn- - *Fe* nec-Aryl

aid* OIf at14 IdentlI by block n.,.h.,)

.)A data structure is described that stores only the non-zero

terms of the assembled stiffness matrix. This storage scheme
results in considerable reduction in memory demand during the
assembly phase of a finite element program. Therefore, larger
matrices can be formed in the main memory of the computer.
When secondary store must be used this approach reduces the

D IFAN"I' 1473 EDITION OF I NOV B5 IS OBSOLETE Unclassified

SECURITY CLASSIFICATION OF THIS PAGE 1
B 00 II..
 Unclsified
Sacush?Y CL*UPCATION OF THIS 10,66(Uhe owe .£...)

/O cost during the assembly stage.

1.
An algorithm is derived that starts with the element con-
nectivity information and generates the compacted data
structure. The element matrices are then assembled to form
the stiffness matrix with this storage scheme. The assembly
algorithm is described and a FORTRAN listing of the routines
are presented. The reduction in storage is demonstrated with
the aid of numerical examples.7

cos duig
I/O
Ii) easebysae
C ...
Anrali tmsdeiethtsatwihheeentc
Unnclassified -
neciviy
iforatin adgnertesUneclassifed dt
structure.
The element matrices are then assembledPto form

$SccCUP /% IT,,CtlAssiCAcoOi or twl P er oweonm

 CONTENTS

Page

INTRODUCTION ............. ........................... 2

STORAGE SCHEME ............ .......................... 3

DERIVATION OF THE ASSEMBLY PROCESS ........ ................ 4

NUMERICAL RESULT ............ ......................... 8

CONCLUSION.......................... 9

REFERENCES .......... ............................ ... 13

APPENDIX A. Program Listing ......... ................... 14

List of Figures

Figure 1. Finite element mesh of example 1 .... ............... 10

Figure 2. Nodal graph for the mesh of example 1 .. ............ 10

Figure 3. The graph of unknowns for the mesh of example 1 .... 11

List of Tables

Table 1. Connectivity sets and active degrees of freedom for the

mesh in Fig. 1 .......... ..................... 5

Table 2. The result of algorithm for establishing the row index

of the nonzero terms in K for example 1 .... ......... 6

Table 3. The row indices of the nonzero terms in the upper

triangular part of K ......... .................. 7

Table 4. Comparison of the storage demands of profile and

compacted assembly for K in examples 2 to 6 ... ....... 8

Table 5. Estimated storage needs for each scheme on regular

mesh. q is the number of equations ....... .......... 9

Table 6. Description of test examples ....... .............. 12

4!

0v V5
 -1-

Introduction

The analysis of large structures, especially in three dimensions, can result in .if:

matrices that demand an exceptionally large amount of computer storage. The storage needs of

these matrices depend to a large extent on their sparsity and the data structure that is used t.

store them. The choice of the data structure in turn depends on the method that is usc , ..

the associated system of equations. Presently, most solution schemes used in finite element com-

puter programs are based on direct methods, i.e. triangular factorization of the stiffness matrix,

K. Starting from a given mesh description, a finite element program performs the following steps:

1. determine the sparsity structure of K,

2. renumber the equations to reduce the storage demands of K,

3. reserve the required storage for K,

4. compute the element matrices,

5. assemble the element matrices into K,

6. compute the triangular factorization of K,

7. solve the associate system of equations.

In many applications the available primary memory is not sufficient to store the assembled

matrix, and therefore secondary storage is used. In this circumstance, steps 4 through 7 involve

data transfers between primary and secondary store, often referred to as 1/0. In this case, K is

partitioned into blocks and each block is assembled and stored on secondary store. The blocks

are then brought back into the main memory to form the factors of K. For large enough prob.-

lems the 1/O costs can dominate the computation costs.

A great deal of effort has been expended to develop new procedures for reordering the equili-

brium equations, thus reducing the overall storage requirements in the solution steps 5, 6 and 7

11,2,M. This ismotivated by the fact that reduction in storage translates directly into a savings

in the I/O costs. Among the many solution schemes used, the frontal method 141and the prolle

or skyline method II are probably the most popular. In 181it is shown that when the same nodal

V.l
 -2-

elimination ordering is used the profile method performs the same number of operations a., lot

frontal method; in [51 an algorithm is described that delivers a good frontal node ordering for the

profile method. The significant difference between these methods is that the frontal method otel.

combines steps 4, 5 and . In this way the I/O during the assembly step is overlapped with tht

1/O in the factorization step; thus, the frontal method results in a saving that is equal to the c.

of 1/O in step 4. Alternatively, the assembly of a profile stored matrix which is partitioned into

blocks requires a multiple pass through the elements to perform the assembly. The principle

differences in the I/O costa of the two methods during factorization may be traced to the 3bove

differences in the partitioning of the matrix into a frontal or a profile form. The principle disa4-.

vantage of the frontal method is the added overhead to retain a small front width during the tri-

angular factorization in step 6 and subsequent resolutions in step 7. For example, in resolution of

equations this added overhead may lead to CPU costs which are several times those of a profile

stored resolution.

In this paper we take a different approach. We use the following simple observation:

The assembly process I Independent of the solution procedure.

In other words, one should use the most efficient data structure for the assembly process, step 5,

and then restructure the data for ones favorite solution scheme, i.e., either frontal or profile. In

this way one can achieve the same reductions in 1/0 as the frontal method and at the same time

maintain high modularity of the program. Here we develop a data structure that stores only the

nonzero terms in the stiffness matrix in a compacted form, and present an algorithm for the

assembly of K for this storage method. This approach results in considerable reduction in the

storage needs during the assembly process. Therefore large matrices often may be fully assembled

in-core resulting in a considerable reduction in I/0.

This approach has the added advantage that the program is not built around a single equa-

tion solver. One can have many solution procedures by simply expanding the compacted strue-

ture of K into a form appropriate for each particular solution method. Furthermore, the com-

pacted structure can be used directly for iterative solution techniques such as the conjugate
 -3-

gradient type methods 16,71.

Storage Scheme

We now describe the compacted structure that is used to store K. We only consider sym-

metric matrices, although the extension to the nonsymmetric case is trivial, and store only L,-

upper triangular part of K is stored. The diagonal terms of K will be stored separately in a eta-

gle array of length n, where a is the total number of equations. The remaining off-diagonals will

be placed in a second array of length r, where r is the total number of nonzero off-diagonal terms

in the upper triangular part of K. All the elements in the same column will be placed consecu-

tively in this array, starting from the top of the column down to the diagonal (excluding the diag-

onal term). The columns are stored consecutively from the second to the last. For each entry in

this array we store its row number in a corresponding integer array of length r. The example

below demonstrates the final storage scheme.

Ezample:

Consider the matrix

9 12 15
2 10 11 13 14
3 16 17
4 19 21
K 5 1820
6 22
7 23
8

The array die# contains the diagonal terms of K as shown below.

dis= 11,2,3,4,5,6,7,81.
A real array then stores the off-diagonal terms of K and a corresponding integer array denotes the

row number of each of-diagonal term as show below.

-oi- die# - 19, 10,11,12,16,13, 17,18, 14,19,20, 15,21,22,231

irow- 11, 2, 2, 1, 3, 2, 3, 5, 2, 4, 5, 1, 4, 6, 71
In addition, a single integer array of length a is also required to point to the end of entries from a

given column. For the above example thi array is

 -4-

icol 10, 1, 2, 3,5, 8, 11, IS)

The total storage requirement is r + n real words and r + n integer words. Using a 16
5
bit integer word, -(n + r) real words (64 bit) will be sufficient. Then the largest number of

equations that can be solved this way is 215 - 1 P 32000. With a 32 bit integer word, t, .-
,c

storage required will be .- (n + r) real words.

Derivation of the Amembly Procs

In this section we give a step by step derivation of the assembly algorithm. Each step is

demonstrated with the aid of the mesh example in Fig. 1. First we introduce some notation. A

finite element mesh is denoted by M - (E , N) where E and N represent the collection of ele-

ments and nodes in the mesh. A part of the input information provided to a finite element pro-

gram is the set of nodes belonging to an element. This we denote as NCN. For example ele-

ment 4 in Fig. I has the connectivity set N 4 = {5,8,9,6). In Table I we give the complete list

of the connectivity sets N, for each element in example 1. These data are usually assembled in a

single array known as the connectivity array. The complete set ( N, V sEE ) is sufficient to =

describe the connectivity of a given mesh. Another part of the input data is the boundary condi-

tions that determine the set of the indices of all the active degrees of freedom at node p. We

denote this set as U.. In Table I we give the set ( Up VpEN) for the example in Fig. 1. In this

example, we assume that there are two degree of freedom at each node. The collection of N,

column 2, and U., column 4, given in Table I is sufficient to determine the sparsity structure of

the stiffness matrix associated with a given mesh.

b:A
 Element Connectivity Boundary Conditions

Element Set of Nodes for Node Set of Unknowns for

No. e each Element, N, No. p the Active nodes, U.

1 (1,4,5,2) 2 (1)
2 (2,5,6,3) 4 (2)
3 (4,7,8,5) 5 (3,4)
4 (5,8,9,6) 6 (5)
7 (6)
675 {(7,11,8)
10,11l, 7)
(8,11,9)
89 (7,8)
(9)

8 (11,12,9) 11 (10)

Table 1. Connectivity sets and active degrees of freedom for the mesh in Fig. 1.

Our objective here is to find the set of indices of the unknowns that are coupled with a

given degree of freedom. This is precisely the row number of each nonzero term in a given

column of the stiffness matrix, irow, that is required for the storage scheme described in the pre-

vious section.

First, we must establish the set of elements that are connected to each node. This can be

done by inspecting the element connectivity sets. Looking at the second column of Table 1, for

example node 4 appears twice, in rows 1 and 3. We then conclude that node 4 is connected to

elements 1 and 3. This process must be repeated for each node. The difficulty here is that we do

not know apriori the number of storage locations needed to identiry the set of elements for each

node. For this reason the above process iscarried out in two steps. The number of elements con-

nected to each node is determined and stored first. We refer to this as the E-degree (element

degree) of each node. In the example the -degree of node 4 is 2. The E-degree also determines

the length of the array that is required to keep the set of elements connected to each node. For

each node p we demote this set by Ep C E. Then the above process is simply to evaluate

E, -c Ip N) (E)
for each node p. This equation may be thought of as fnding the pseudo-inverse of the connec-

tivity array. The &,degreeof mode p is the number of terms inE (the cardisality of E. ). See
 -6-

columns two and three of Table 2 for the Edegree and the complete set of E, for the nodes in

example 1.

Next, we And for the set of nodes that are adjacent to each node p. We denote thib L1

ACN. This is the set of all nodes that belong to an element with p as one of its nodes IR*:vi-7
established the set of elements connected to node p, the adjacent nodes are all the 4t'

belonging to these elements. In example I node 4 is connected to elements I and 3 (see column 2

of Table 2). The set of nodes belonging to elements 1 and 3 are obtained by inspecting column 2

of Table 1; these are ( 1, 4, 5, 2 ) and ( 4, 7, 8, 5 ) respectively. Then the set of nodes adjacent

to 4 is A 4 { 1, 5, 2, 7, 8 ). The N-degree (nodal degree) of a node is the number of nodes adja-

cent to it and is the cardinality of A,. In columns 4 and 5 of Table 2 we give the N-degree and

the adjacency set of each node in the nodal graph of example 1 (Fig. 2). This step is simply to

evaluate the equation

A,= U N, - n (2)

Note that both A, and N-degree can be obtained in the same loop.

Node E N S
p degree Ep degree AP degree SP

1 1 1) 3 (4,5,2) 0
2 2 (1,2) 5 1, 4, 5,6,3) 1 {
3 1 (2) 3 (2,5,0) 1 (2)
4 2 (1,3) 5 (1,2,5,7, 8 2 (1,2)
5 4 (1,3,4,21 8 1, 4, 7, 8, 9,6, 3, 2 4 (1,4,3,2)
* 6 2 (2,4) 5 (2, 5,8,9, 3 ) 3 (2,5,3)
7 3 (3,5,6) 5 (10, 11,8,, 4) 2 (4,5)
8 4 (3,5,7,4) 6 (11, 9,6, 5, 4, 7 4 (6,5,4,7)
9 3 (4,7,8) 5 (11, 12, 6, 5,8) 3 16,5,8)
10 1 (6) 2 17,11) 1 (7)
11 4 (6,5,7,8) 5 (10,12,9,8,7) 4 (10,9,8,7)
12 1 (8) 2 (11,9) 2 (11,9)

Table 2. The result of algorithm for establishing the row index of the nonzero terms in K
for example 1.
 -7-

Since we want to store only the upper triangular part of K we need to store only a subset of

AP. This will be the set of nodes in A, with an index less than p; that is

Sp= (i iEA, and i < p ). We refer to the number of terms in S, as the S-degree (.,rc.-

degree) of a node. The set S, is only useful when the numbering of the unknowns are such tha-t

when i < j all the unknowns at node i have a smaller index than the urknowns at n,,Ie

Whenever this is not true it is necessary to use the complete set of adjacent nodes A. together

with the numbers of the unknowns for each node (e.g., see listing in Appendix A.).

Finally, for a given unknown at node p with index EUV, we find the set of the indices of all

other unknowns that are coupled with u,. This will be the set of row indices R, for nonzeros in

-i-th column of K. Then

R, U U, (3)
tIE A,

For example 1 ffj is the adjacency set of j in the graph of the unknowns in Fig. 3. Since we only

store the upper triangular part of K we scan through ifj and use the subset defined by:

R,=( i I iEf, and i <j) (4)

R) is the row index of all the nonzero terms in the i-tb column of the upper triangular part of K.

The complete set of R, for the example problem is presented in Table 3.

Ij D.O.F. Node
p R

1 2 (0)
2 4 11)
3 5 (1,2)
4 5 (3,2, 1
6 6 (4,3,1)
6 7 (3,4,2)
7 8 (6,2, 4, 5, 3)
S 8 2, 3,4,,6,,)
0, (7,4, 5, 3, 8
10 11 (9,6,7,8)

Table 3. The row indices of the nonzero terms in the upper triangular part of K.

_ r
 -8-

The listing of a FORTRAN program that performs all the steps that is described in Lbis

section is given in Appendix A. In this Appendix we also provide the subroutine that use, the

row indices to perform the assembly of K.

Numerical Result

We use the algorithm described in the previous section to assemble the stiffnems matr!,,

the problems described in Table 6. The total storage required during the assembly step is

evaluated. We compare these results to similar results obtained when the assembly is performed

directly into a profile data structure. The storage requirement of the compacted assembly is not

effected by the node ordering. For the assembly into a profile form, we numbered the nodes

across the width of the mesh to reduce the bandwidth of the stiffness matrix. Although, the

bandwidth could have been reduced further using a renumbering scheme such as 12,31, we omitted

this step for simplicity. The results for examples 2 to 6, given in Table 6, are presented in Table

4 below.

Description No. of No. of No. of Stores for Stores for

of Problem nodes elements equations profile Compacted K

Cantilever type structure 225 184 428 10204 5340

Small Cylinder structure 231 200 440 10492 5679

Large Cylinder structure 496 450 960 32542 12715

4X4X4 solid structure 125 64 300 21945 11634

8X8X8 solid structure 729 512 1944 470043 94272

Table 4. Comparison of the storage demands of profile and compacted assembly for K in
examples 2 to 6.

The reults in Table 4 is obtained based on the assumption that a real word is twice as long as an

integer word. We observe a reduction from 40% for two dimensional (2-D) problems to 80% for

3-D problems for these examples. The reductions will be more if short integer words are used. It 4
is interesting to note that the required storage for compacted structure varies linearly with the

number of equations. Therefore, the saving will be more for larger problems. In Table 5 we give
the storage counts for the two methods considered here on square mesh in 2-D and c ' .

a function of the number of equations. To obtain these estimates we assumed that there ib ,ady

one degree of freedom per node.

Dimension Halft Profile Compacted

of Problem bandwidth storage storage

1 2 5/2 n 3n

2 n/ 15/2 n
3 n 2/ n6/3 21n

Table 5. Estimated storage needs for each scheme on regular mesh. n is the number of equations.

Conclusion

The essential steps in a finite element program can be modified to make use of the com-

pacted assembly described here. Accordingly, we perform the following steps:

1. Obtain the row indices of the nonzero terms in K.

2. Assemble the matrix in compact form.

3. Choose a solution procedure and renumber the equations to reduce the storage

demands of the factors of K.

4. expand the compacted K into a data structure suitable for the solution method.

S. solve the associated system of equations.

When there is insufficient primary storage, the assembly of the matrix in compacted form

opens a number of possible avenues that one can take to reduce the 1/o cost. The expnamsloa of

the compacted form need not be done immediately after it's assembly. The matrix can be kept in

compact form and put on secondary store and expanded into a full prolle form only when a fac-

torization mut be performed. This way the number of data entries that is read (in the lpet

phse of I/O) cam be reduced considerably, which is turn results in a reductiom is the solution

Ume.
 * 10-

10 4 7

Figure 1. Finite element mesh of example 1.

Figure 2. Nodal graph for the mesh of example 1.

 I
-11-

$ -

Figure 3. The graph of umknowns for the mesh of example 1.

;. 1

I
 - 12-

., i
Example 2: Cantilever Struc-
ture, left end fixed, plane stress
elements with 2 degrees of free- .,
do. per mode.

Example 3: Small Cylinder,

both ends fixed in tangential
direction, plane strain elements
with 2 degrees of freedom per
node.

Example 4: Large Cylinder,

both ends fixed in tangential
direction, plane strain elements
with 2 degrees of freedom per
nde.
lt

Example 5: 4X4X4 Solid cube,

fixed base, solid elements with 3
degrees of freedom per node.

Example 6: 8X8X8 Solid cube,

fixed base, solid elements with 3
degrees of freedom per node.

Table 6. Description of test examples.

 -13-

Reference[

III A. George and J. W. Liu, Computer Solution of Large Sparse Positive Definite Systems,
Prentice-Hall, Englewood Cliffs, 1981.

121 E. Cuthill and J. McKee, "Reducing the Bandwidth of Sparse Symmetric Matrices," Proc.
ACM Nat. Con., New-York, 1969.

131 N. E. Gibbs, W. G. Poole, Jr., and P. K. Stockmeyer, "An Algorithm for Reducing Lbe
Bandwidth and Profile of A Sparse Matrix," SIAM 1. Num. Anal., Vol. 13, pp. 236-250,
1976.

14) B. Irons, "A Frontal Solution Program for Finite Element Analysis," Int. J. Num. Meth.
Engng., Vol. 2, pp. 5-32, 1970.

151 M. Hoit and E. L. Wilson, "An Equation Numbering Algorithm Based on a Minimum Front
Criteria," Computers and Structures,Vol. 16, No. 1-4, pp. 225-239, 1983.

161 B. Nour-Omid, B. N. Parlett and R. L. Taylor, "A Newton-Lanczos Method for Solution of
Nonlinear Finiie Element Equations," Computers and Structures, Vol. 16, No. 1-4, pp. 241-
252, 1983.

171 R. L. Taylor and B. Nour-Omid, "Solution of Finite Element Problems by Preconditioned

Conjugate Gradient and Lanczos Methods," Rep. No. UCB/SESM-8405, Department of
Civil Engineering, University of California, Berkeley, Mny 1984.

181 R. L. Taylor, E. L. Wilson and S. Sackett, "Direct Solution of Equations by Frontal and
Variable Band, Active Column Methods," in Nonlinear Finite Element Analysis in Structural
Mechanics, Proc. Europe-US Workshop, July 1976.

191 R. L. Taylor, "Computer Procedure for Finite Element Analysis," Ch. 24, The Finite Ele-
ment Method, 3-rd Ed., by 0. C. Zienkiwicz, McGraw-Hill, London, 1977.

1101 E. L. Wilson, "Solution of Sparse Stiffness Matrices for Structural Systems," Proc. of Sparse
Matrices, Ed. 1. S. Duff, SIAM, Philadelphia, 1979.

m&
 -14-

Appendix As Program Lsting

SUBROUTINE ELCNT(NUMNPNUMEL,NEN,NENI ,IX.IC)
DIMENSION IX(NENI,I),IC(l)
C .... INPUT
C NUMNP TOTAL NO. OF NODES IN THE MESH
C NUMEL TOTAL NO. OF ELEMENTS IN THE MESH
C NEN MAX. NO. OF NODES PER ELEMENT
C NENI DIMENSION OF IX ARRAY
C IX ELEMENT CONNECTIVITY ARRAY
C
C.... OUTPUT
C IC ARRAY OF LENGTH NUMNP. IT FIRST HOLDS THE ELEMENT DEGREE
C OF EACH NODE, THEN BECCMES A POINTER FOR AN ARRAY THAT
C CONTAINS THE SET OF ELEMENTS CONNECTED TO EACH NODE.
C
C.... COUNT THE NUMBER OF ELEMENTS EACH NODE BELONGS TO
C
CALL IZERO(IC,NUMNP)
DO 110 N = I,NUMEL
DO 100 J = INEN
I = IX(J,N)
IF(I.GT.0) IC(I) i C(I) + I
100 CONTINUE
110 CONTINUE
C
C.... SET UP POINTERS
C
DO 120 I --2,NUMNP
IC(I) = IC(I) + IC(I-I)
120 CONTINUE
C
RETURN
END

SUBROUTINE CASSEM(D,A,BS,P,JCOLE,IROW,LD, ID,NSTNEL,AFL,BFL)

IMPLICIT DOUBLE PRECISION (A-H,O-Z)
LOGICAL AFL,BFL
DIMENSION D(I).A(l),B(I),S(NSTI).P(1),JCOLE(I).IROW(1).LD(1)
I ,ID(1)

C ... COMPACT ASSEMBLY OF PROFILE MATRIX

c
DO 200 J = INEL
N = LD(J)
IF ( AFL AND. N .GT. I ) THEN
DO 150 I - INEL
K - LD(I)
IF ( K GT. 0 -AND. K .LT. N ) THEN
INZ = INZA(JCOLE(N.I)+I,JCOLE(N),IRC.,K)
A(INZ) - A(INZ) + S(I,J)
END IF
150 CONTINUE
END IF
C .... ASSEMBLE THE DIAGONAL
IF ( N GE I ) THEN
IF ( AFL ) D(N) - D(N) + S(J,J)
C .... ASSEMBLE THE LOAD IF NECESSARY
IF ( BFL ) B(N) -B(N) + P(J)
END IF
200 CONTINUE
RETURN
END
 - 15

SUBROUTINE CC3URO(NUh4PNUNP ELNEN.NEN NDF, IX, ID, IC, IROW. IELC.

__ JCOLE.KP)
DIMENSION IX(NENI ,I). ID(NDF, I), IC( 1), IRCN( I), IELC( 1), JCOLE( 1)
C
C FOR (NUMNP,NUMEL.NEN,NENI , IX, IC) SEE SUBROUTINE ELCNT
C INPUT
C NDF NUMBER OF UNKNOWNS AT EACH NODE
C ID ACTIVE UNKOWNMS AT EACH NODE
C. OUTPUT
C IELC HOLDS THE SET Of ELEMJENTS CONNECTED TO EACH NODE
C IROW ROW NUMBER OF EACH NONZERO IN THE STIFFNESS MATRIX
C JCOLE END OF ENTRIES IN IROW FORM A GIVEN COLUMN
C
C FIND ELENTS CONNECTED TO NODES
C
CALL IZERO (IELC.IC(NUMP))
DO 230 N = I,NMEL
DO 220 J = INEN
I = IX(J.N)
IF ( I .GT 0 ) THEN
KP IC( I)
200 IF (IELC(KP) EQ. 0 ) 0O TO 210
KP = KP I
GO TO 200
210 IELC(KP) = N
END IF
220 CONTINUE
230 CONTINUE
1C
C SET UP CO&/RESSED PROFILE POINTERS
C
IKP =
NEP =]
0

DO 360 1 = I,NUMNP
NE = IC( I)
DO 340 I = I,NDF
NEQ = ID(II,I)
IF ( NEQ GT. 0 ) THEN
JCOLE(NEQ) = KP
KPO = KP + I
IF ( NEP LE. NE ) THEN
DO 330 N = NEP,NE
NN = IELC(N)
DO 320 J I,NEN
K - IX(J,NN)
DO 310 JJ = INDF
NEQJ - ID(JJ,K)
IF (NEQJ GE NEQ OR NEQJ LT. 0) GO TO 310
C
C.. CHECK TO SEE IF NODE ALREADY IN LIST
C
IF ( KPO LE. KP ) THEN
DO 300 ICK = KPO,KP
IF( IRON(KK) EQ NEQJ ) GO TO 310
Soo CONTINUE
END IF
lCP - KP + I
IROW(KP) - NEQJ
310 CONTINUE
320 CONTINUE
330 CONTINUE
JCOLE(NEQ) =K?
END IF
END IF
340 CONTINUE
NEP - NE + 1
350 CONTINUE
RETURN
END
 -16-

INTEGER FUNCTION INZA(NI ,N2, IR.,K)

C DIMENSION IROW(1)
- C
C FIND THE TERM FOR THE ASSEMLY
C
DO oo N=Ni,N2
IF ( IRoW(N) EQ. K ) THEN
INZA = N
RETURN
END IF
100 CONTINUE
C.... ERROR IF LOOP EXITS
STOP
END

SUBROUTINE IZERO(IA,NN)
DIMENSION IA(NN)
DO l00 N = I,NN
IA(N) = 0
100 CONTINUE
RETURN
END

i[
t_

e . . .. l I