Carlos A. Felippa -- Solution of linear equations+

Computers & Srructurrs. Vol. 5. pp. 13-29. Pergamon Press 1975. Printed in Great Britain
SOLUTION OF LINEAR EQUATIONS
WITH SKYLINE-STORED SYMMETRIC MATRIX
CARLOS
A.
FELIPPA
Structural Mechanics Laboratory, Lockheed Palo Alto Research Laboratory, 3251 Hanover Street, Palo Alto,
California 94304,U.S.A.
(Received in reoisedfom
26March 1974)
Abstract-Fortran IV subroutines for the in-core solution of linear algebraic systems with a sparse,
symmetric, skyline-stored coefficient matrix are presented. Such systems arise in a variety of applications,
notably the numerical discretization of conservative physical systems by finite differences or finite element
techniques. The routines can be used for processing constrained systems without need for prearranging
equations. The application to ‘superelement’ condensation of large-scale systems is discussed.
finite element methods. In the latter case, direct elimination techniques have been used to solve the resulting
linear systems Ax = b since the advent of the finite
element method in the mid-1950’s[ l] (this is in contrast to
classical finite difference discretization, where iterative
solution techniques still prevail). Commonly used algorithms include Gauss elimination, standard Cholesky
(A = LL’) and modified Cholesky (A = LDLT). The main
implementation differences, however, concern the manner in which nonzero coefficients of A are arranged in the
various storage devices (e.g. core, extended core, random
access mass storage, tapes), utilized during the solution
process. Excellent surveys on this subject are presently
available [2-4].
The routines presented here utilize only one storage
level, namely, high-speed memory. The storage arrangement of matrix A is the so-called skyline, profile,
envelope, or variable
bandwidth
method[5-61. The
symmetric matrix A is stored in a linear array as a string of
‘active’ columns of the upper triangle (or, equivalently,
rows of the lower triangle). The active portion of each
column is bounded by the diagonal and the furthest
nonzero element.
It should be mentioned that more sophisticated sparse
matrix storage schemes are available [2-4,7-91. However,
implementation of most such schemes requires considerable programming skill to ensure that the data handling
overhead (addressing, fetching, storing, etc.) is minimized. The skyline storage scheme offers a reasonable
compromise between organizational simplicity and computational efficiency. Moreover, since the associated
solution process consists largely of vector inner products,
it is a natural choice for implementation on the
‘fifth-generation’ parallel and pipeline computers [IO].
NOMENCLATURE
A symmetric coefficient matrix
A,,, A,,, A,< block partitions of A in constrained system
A* condensed A upon elimination of xf in superelement
condensation process
B A,, -A*
C(A) Euclidean condition number of A
L lower triangular matrix
D diagonal matrix
U upper triangular matrix
X see equation (19)
b right hand side vector
b<,b, prescribed and unknown portions of b, respectively,
in constrained linear system
6 modified b, in constrained system solution
b* condensed b upon elimination of xf in superelement
condensation process
r residual vector b -Ax
x solution vector satisfying Ax = b
x0 xr prescribed and unknown portions of x, respectively,
in constrained linear system
x* solution iterates produced in iterative refinement
process
A.& solution corrections in iterative refinement process
Y>= intermediate vectors in solution process, equation (3)
A linear array containing A or factorization thereof
LD array of pointers to diagonal elements of A in A
o,, uii (i, j)-th element of A, U
b,, xi j-th element of b, x
di j-th diagonal element of D
k number of iterative refinement cycles
n order of A
nf, n, number of free and constrained equations, respectively
In number of random r.h.s. vectors used for estimating
C(A)
ri length of j-th row of A excluding constrained
columns
ti tolerance for singularity test
floating point machine accuracy
B relative solution error
Il.1122-norm of matrix or vector
ll~llf Euclidean matrix norm
matrix inverse
‘;‘I.T matrix or vector transpose
ALGORM'HMDFSCRIPTION
Purpose
Subroutines SKYFAC and SKYSOL are designed to
solve the linear equation system
Ax=b
lNTBODUCTlON
where A is an n by n sparse symmetric nonsingular
matrix entirely core stored in skyline form, and x and b are
column vectors. The routines can handle the case in which
some of the components of x are prescribed.
Large sparse symmetric systems of equations generally
arise in practice as a result of numerical discretization of
self-adjoint problems by variational finite difference or
C.A.S..Vol. 5. No. I-B
(1)
13
C.
14
A.
Method
The symmetric factorization
in which entries outside the skyline template are not
shown. Matrix (6) is stored as a IS-work string A:
A=LDU
(2)
where L is a unit lower triangular matrix, D a nonsingular
diagonal matrix and U is the transpose of L, is carried out
by subroutine SKYFAC by a compact Cholesky-type
factorization algorithm similar to that described in[ll,
Contrs. I/l and I/4]. The original matrix A is overwritten
by elements of D-’ and U in the manner discussed in the
‘Organizational Details’ subsection. No pivoting is used
in the factorization process.
Once A has been factorized, subroutine SKYSOL can
be repeatedly invoked to solve equation (1) for any right
hand side b. The solution vector is obtained by the usual
three-stage process:
Lz=b
(forward substitution)
y = D-‘z (scaling)
ux=y
FELIPPA
(3a)
A: aI;, a22, alj, 0, a33, a24, aj4, ati, a55, a16,
0, 0,a46,a56, a%.
This array is complemented by a (n + 1)-word array LD
containing diagonal location pointers:
LD: 0, 1,2,5,8,9,15.
I/d,
(3c)
Constrained system
SKYFAC-SKYSOL can also be used to solve the
more general case of a constrained linear system:
where both b, and x, are prescribed (subscripts f and c
stand for free and constrained, respectively). Equations
associated with xc are called constraint equations or
simply constraints; they result from the discretization of
essential boundary conditi0ns.t A constraint equation
corresponding to a prescribed zero (nonzero) component
of x is called homogeneous (nonhomogeneous). If all
constraints are homogeneous, i.e. xc = 0, the system (4) is
said to be homogeneously constrained.
System (4) can be solved for xf upon reduction to the
form
Arrx, = bf - Af,xc = b,
(5)
where l;f = b, if the system is homogeneously constrained. I’he subvector b, may then be obtained from the
second matrix equation in (4), if desired.
It is important tb note that it is not necessary to
physically arrange a constrained system (1) to put it into
the partitioned form (4). Constrained and unconstrained
equations may be freely intermixed.*
Organizational details
The storage arrangement of the coefficient matrix A is
best illustrated by a simple example. Consider the (6 X6)
symmetric matrix:
011
A=
aI3
aZ20 a24
a33 a34
a4
symm.
[
aI6
0
0
a46
as5 as6
ah6I
is often called the constraint reaction vector.
SIf the number nc of constrained equations is significant (say
n, B 0.3 n), it may be computationally advantageous to preorder
(1)into the form (4) before presenting it to SKYFAC-SKYSOL.
tb,
u13
lldz uz ~24
l/d, ~34
Ul6
U26
U36
l/da
(9)
lids :::
l/d61 .
I
Constrained equations are flagged by a negative LD
pointer. For instance, if x3 and x5 are to be prescribed in
the linear system with coefficient matrix (6)
LD:O,1,2,-5,8,-9,15.
(10)
In this case, the decomposition (2) is performed only on
the unconstrained equations of (4), i.e. A,, is replaced by
Dff’ and U,. Constrained rows and columns are not
altered.
Decomposition failure
The factorization process (1) may be aborted if one of
the following conditions is encountered.
1. Singularity: matrix A (Aff if constrained) appears
singular to machine accuracy. The singularity test
performed at the j-th elimination stage is[ 11,Contr. I/7]:
di < tj = 8 El;
(11)
where d, is the j-th diagonal entry of D, E the smallest
positive floating-point number for which 1+ l > 1 on
computer being
and
is the Euclidean norm
(length) of the
row of
excluding constrained
columns.
2. Indefiniteness:
A
if
appears
be
and check
positive
is
by
user.
j-th
test
di
(6)
(8)
For a genera1 matrix, LD(I + 1) points to aii in array A
n, whereas LD(l) contains zero.
for i=l,...
During the factorization process, elements of the strict
upper triangle of U replace the corresponding superdiagonal entries of A, and the reciprocals of the diagonal
elements of D overwrite the diagonal elements of A. Thus
the factorization (2) of the example matrix (6) would be
arranged as follows:
(3b)
(backsubstitution).
(7)
-
(12)
The
stability the
factoriza(3) not
if
is
positive
However,
many
situations,
as
finite
analysis
a
symmetric
tion
an
matrix guaranteed
some
restrictions
the
of
equations
observed.
15
Solution of linear equations with skyline-stored symmetric matrix
Solution accuracy estimation
Two optional mechanisms for assessing the accuracy of
computed solutions are provided:
1. Matrix condition evaluation (a priori estimate). A
lower bound estimate of the Euclidean condition number
of A[13, Sec. 3.111:
2. Problems involving local nonlinearities, if these can
be confined to the bottom of the discretization operators.
Then the (linear) upper portion of A can be factorized
once and for all during the nonlinear iterative solution
process (a similar situation occurs in the implicit time
integration of locally nonlinear dynamic systems).
USAGE
C(A) = ]IA~I&‘!E
(13)
where ]]jE denotes the Euclidean matrix norm, is returned
by SKYFAC on request. The estimation is based on
performing one cycle of iterative refinement on a specified
number m of ‘random’ right hand side vectors b, and
applying equation (13.3)of chapter 3 of [ 131.The larger m,
the better the estimate, but in practice m = 2 is often
sufficient. If C(A)E & 1, log10C(A) may be used as a gross
estimate of the number of significant digits expected to be
lost in the solution process.
2. Iterative refinement of solution (a posteriori estimate). Given the initial computed solution x0 of (I), a
k-cycle refinement process is defined by the algorithm [ 11,
Contr. I/2]:
ri = b-Axi_,
i=l,...k
LDUAxi = ri
xi = xi-1 t Ax, 1
The routines presented here include the equation solver
proper (SKYFAC-SKYSOL) and a test package (SKYTEST). Relationships among the various routines are
illustrated in Fig. 1.
r---------
(14)
where only the calculation of the residual ri has to be
carried out in higher precision arithmetic. Given k >O,
subroutine SKYSOL returns xk and the relative error
Sk= llhxkllzlllxt11z.
If S, 4 1, 8, provides an estimate of the accuracy of xkml.
Note: In practice it is pointless to specify k > 1 unless the
entries of A and b are exact within their single precision
representation (as it happens, for instance, in some finite
difference schemes on regular grids), since the sequence
xe will not generally converge to that solution of system
(1) associated with infinite precision representation of A
and b.
Partial factorization
Subroutine SKYFAC allows the user to specify that the
decomposition of A is to be performed from equation
(NBEG t 1) through NEND, where NBEG, NEND are
input parameters. The usual complete factorization of A is
achieved by setting NBEG = 0 and NEND = matrix
order. A partial factorization may be requested by
suitably constraining that range. For instance, if matrix (6)
is submitted to SKYFAC with NBEG = 0 and NEND =
3, the output is
I/d,
Ml3
lldt u23 a24
l/d, a34
au
ala
0
0
a46
(16)
ass a56
I
aba: .
A subsequent call with NBEG = 3 and NEND = 6 would
produce (9).
This capability is useful in the following cases:
1. The calling program may wish to examine factorization of principal minors (e.g. in local stability analysis of
physical systems).
(
TEST
PACKAGE
L____-_____i
(W
Fig. 1. Subroutine flowchart illustrating relationships among
program modules.
Documentation for using subroutines SKYFAC and
SKYSOL is embodied in the source listing. Both
subroutines utilize an external, user-supplied vector inner
product function whose name is transferred in the formal
argument list. Fortran samples of such function are
provided in the listing, but for maximum efficiency the
function should be handcoded in assembly language
(many computer installations have such a procedure in the
utility library).
The test program SKYTEST is intended to illustrate the
use of SKYFAC-SKYSOL, and to certify their implementation on specific computers. SKYTEST uses four
auxiliary routines: SKYLDU, SKYOPC, SKYMAP and
SKYPRT. The last two are likely to be useful on their own
for matrix display purposes. Guidelines for setting up and
running the verification program are provided in the
comment block heading SKYTEST. A test example
produced by the test program on the UNIVAC 1108 is
also given.
Program modifications
Only one machine-dependent constant is utilized in this
package: the floating-point precision constant E denoted
by EPSMAC in subroutine SKYFAC. Appropriate values
of l for various computers are indicated there.
Conversion of SKYFAC-SKYSOL to operate in full
double precision is straightforward: all REAL declarations should be changed to DOUBLE PRECISION (or
REAL*8 on IBM 360/370 models), or appropriate
IMPLICIT statements inserted. If this is done, the
optional iterative refinement mechanism loses meaning
C. A. FELIPPA
16
SUBRO”,INE
SKYFAO
...
PRLIGRHM
07,
NBEG,
. PDCHEK,
ICOND,
IRFILE,
c
C..*.............*..**.*..*.*...**.*......**.....***.*.............*..**
LISTING
NEND,
FlCLiND,
. . .
LD.
Nr
DETCF.
DOTPRD.
VI
IDETEX,
SINGABs
NEGEIG,
IFFaIL>
_______-__-------
t*
c*
C.
.
.
*
-PU*POSE_______________--
.
C.
C+
TO
PERFORN
SYNMETRIC
IN-CORE
FRCTOR1ZRTION
NHTR~X
OF
.
.
.
A SPHRSE,
SKYLINE-STORED,
SYMMETRIC
C.
C*
C........*..,.*.*.........*..**........*....*........*...*..*.....*.....
C*
C*
.
PROGRRMNED
-
CRRLOS
A.
FELIPPA,
MRHCH
l
1767.
l
::
C.
CO
C.
CO
C*
LlPDRTE
-
JClNE
.
1974.
l
LRNGUAGE
-
FORTRAN
tI+J>-TYPE
IV
<ANSI
ESCEFT
SUBPCRIPTS)
FOR
“RCHINE
EPSMAC,
INDEPENDENT
(EXCEPT
SEE DISCUSSION)
.
.
OCCASIONAL
.
EBUIPNENT
-
FIIH
l
TOLERRNCE
.
.
::
C*.r....***...*****.*.**.........*....*...........***...*.........*.....
----------------------_
- D 1 ; L I, 2 S I 0
N
.
.
.
-
________-
l
THE UPPER
TRIRNGLE
OF THE INPUT
ZYNNETRIC
MHTPIX
CR, IS 3TORED
IN COLUPIN-WISE
‘ShYLINE’
FRSHION
,N THE FIRST
NWA WORDS OF
ARRRY R. THE STORRGE
SCHEME
IS BEST
ILLUSTRRTED
BY THE
(6 X 64
MRTRIX
<DOTS DENOTE
ZEROS
OUTSIDE
SKYLINE,
:
l
.
l
.
l
.
.
.
l
l
.
Au’
I
l
I
l
THIS
IS
PRESENTED
TO
SKYFAC
AS
H
I’+WORD
PTRING
IN ARRRY H.
THE POSITION
OF THE DIRGONRL
ELEMENTS
A tN+1,-WORD
INTEGER
ARRAY LD,
SlUCH THAT
LDcl>
=
LD<I+I,
POINTS
TO R<I,I,
IN RRPRY R FOR I = IclaN.
FOR THE EXAMPLE
MRTRIX
(1)
:
LD
=
Or
I,
2,
5,
3,
9,
.
.
.
:
IS NARKED
U HND
BY
l
.
.
.
.
15
(2R>
.
l
THE
SYMMETRIC
.
.
DECOMPOSITION
CH,
=
[Ll
CD1
GUI
13
.
l
WHERE
WI
15 THE TFRNZP[I:E
OF CL,,
WHICH
11 il iUNIT
LOWER
TRIHNGULW
NRTR,::,
“ND
CD,
I? F, DIHGOllRL
MATRIX,
I;
PERFORMED.
THE STRICT
UPPEY
TRIANGLE
OF WI
OVEFUFITES
THE OFF-DlH6ONHL
THE RECIPCOCHL;
OF THE ELE”ENTS
OF
ENTRIES
IIF
[RI I LIHERE‘,:
OF [D,
DEPLACE
THE DIHGONHL
ELEMENT:
OF 0%.
THE DUTPClT
;TORHGE
ISHEME
FOP THE EXRNPLE
MHTPI::
a 11 I I THEN:
.
THE DECONPO;ITION
PPOCEIC
FOLLOWING
CONDITICiM:
:
.
THE NHTPI.‘:
:ING,,LHF.
“H‘i
BE
ABORTED
UNDEF
EITHER
OF
.
.
.
.
.
.
.
THE
l
rm~*FlEn
[HI.
THE IINGULWITY
ABIID<.,b,
B’i THE R~WNDING
ERRORS.
TE?T
WED
HT THE J-TH
LE
t3PPmRS
3TEP
IS
3.iP’NHC.PIJ~
35,
WHEYE P\,’
DENOTE’I
THE ELlCLlDEHN
NOPM OF THE I-TH
POW OF
[R,,
HND EP5NW:
II
THE :MHLLELT
PO‘lITli’E
FLDRTING-POINT
NUABEF
%I:,,
THRT
1. O+EPt”ili
6,
1.0
OH THE CONPUTEP
THIS
51N61,~HPI,r
TE?T
li
UNL’I’ ENFORCE11
IF
BEING
U:ED.
INPUT
FLHG
JIW?HB
I:
. TRUE.
c:EE
INPUT
HPGUMENTT,.
.
+
.
.
.
.
.
l
.
.
l
l
THE
NHTRIX
CH,
I:
INDEFINITE.
l
I.E.9
l
DC,*
LT
.
S*EP?MRC.RIJ’
l
THE POSITIVE
DEFINITENES?
TE:T
I:
ONL’i ENFORCED
IF
INPUT
FLRG PDCHEh
15 . TFIJE.
a5EE
INPUT
RPGUNENT:I.
IT ZH‘WLD
BE
ZTRESSED
THHT
THE NUNEPli
iTHBILITY
OF THE DECOMPO‘;ITIDN
(3,
15 NOT GURPRNTEED
IF NHTPb
[H,
IT
INDEFINITE.
l
l
l
.
l
.
THIS
SUBFOLlTlNE
MHY RL:O
BE UTILI:ED
TO HHNDLE
THE NORE
GENEPRL
tRCE
OF R CONSTWINED
LINEAP
‘I’rTTEM,
IN WHICH
SOME
THE X-ENTRIEC
WE
PRESCRIBED.
H CONITPHINED
EYTTE”
NW
BE
FOR”kLY
WRITTEN
IN THE PRRTITIDNED
MHTRIX
FORM :
l
BF
l
.
.
+
17
Solution of linear equations with skyline-stored symmetric matrix
c*
C.
C.
C.
C.
CO
C.
C.
C.
C.
CO
,: .
;:
1:l
C.
.
.
WHERE
CBFI
fiND
IXCI
HRE GIVEN
AND CXFl
RND CBC, RRE TO BE
DETERt‘INED
<SUFFICES
‘F’
RND ‘C’
STAND
FOR ‘FREE’
AND
‘CONSTRAINED’,
RESPECTIVELY.
HOWEVER.
IT “L&T
BE EMFHHSIZED
THRT
IA3 NEED NOT BE PHYSICRLLY
RPPRNGED
TO CONFOPM
TO (7,.
.
.
l
IF X(1>
IS TO BE PPESCPIBED.
THE I-TH
POCl RND COLUMN ilF
IA3
ARE S,3ID
TO BE CONSTPRINED.
SUCH CONDITION
I5
NHPIED
BY H
NEGATIVE
VALUE
OF LDaI+l.I.
FOR INTTANCE.
IF RObIS 3 RND 6 OF
EXHNPLE
MHTPIX
(1)
RRE CONSTRRINED.
THE LD HPRRY MUST BE:
LD
=
0,
I,
2,
-5,
8,
‘3,
l
.
THE
.
.
.
(7RI
-15
l
.
.
.
.
.
.
.
CON?TRHINED
PI,@
tiND COLURNJ
HFE IGNORE11
IN THE FHCTOPIZHTION
OF R CONITPAINED
MHTPIX
HND VHLUEr.
THEREIN
PETUPN
LINHLTEYED.
FOR INCTRNCE,
THE OUTPUT
[email protected] ?kYFAC
FOP THE EXHNPLE
NHTPIX
(11
lJITH
THE L” HF‘YAY
‘7R,
r\tOULD BE :
c*
C.
C.
0:.
I
I
I
I
1
I
c*
c.
i.
c.
cl
I ‘DI
I
HI3
HZ3
H3’
112
,_I24
,_,1:4
1, DJ
1,D5
HI&
I
0
n
H46
Fl5c;
“66
I
I
I
I
I
l
IS>
:
.
.
l
.
.
+
C.
c*
c+
THhT
III
OHLY THE ,uNRE~PP~~NGED~
EFFECTIVELY
FHCTORIZED
INTO
CLFF,
FFEE-FREE
POPTION
CDFF,
CUFF,.
IRFFI
IS
l
l
c”:
C.
HDDITIONHL
FEHTlJPEl
OF
THIS
KIITINE
.
.
.
.
.
INCLUDE:
c.
C.
THE DETEF”INRNT
BE CONPUTED
ON
1.
c+
OF
CHJ
,OF
IHFF,
,F
CONSTRHINED1
“AY
FEOUEST.
.
0.
C.
2.
FbN ESTINHTE
OF THE EUCLIDEHN
CHFFI
IF CONSTRAINED,
MHY BE
c.
CONDITION
OBTHINED.
NUMBER
OF
IHI
.
iOR
l
C.
.
c.
r:*
3.
HCCUFIIJLHTION
OF INNER
PRODUCT’ i IN DIFFERENT
PRECISION
MOIlE5
IT FHCILITHTED
BY INCLUDING
THE INNERPRODUCT
FUNCTION
IN THE FORNHL
RPGUNENT
LIST.
VHRIOLIS
ENTRY
POINT
NRMES
MHY THUS BE UIED
WITHOLIT
NEED FOR EDITING4ECOMPILING.
C.
C.
2:
C*
:
t
l
l
4.
R CHPABILITY
FOR PF1CTIHL
FACTORIZRTION
IS
INCLUDED.
T,,E
FHCTORIZHTION
‘3’
I:
HCTUALLY
CRPRIED
OUT FROM EQURTION
*NBEG+Ia
THROUGH
rlENDr
WHERE NBEG,
NEND RRE TWO INPUT
PRRHMETERS.
NORIIALLI’
NBEG=O
RND NEND=N.
INATRIX
ORDEk.
HOWEVEF.
IN CERTHIN
CHZES
(E.G..
LOCHL
STRBILITY
ANRLYSIS,
LOCHL
NONLINEHKITIES;~
I H PRPTIAL
FRCTOPIZRTION
MRY BE
UCEFUL.
c*
c*
C.
$1.
c*
c*
.
.
.
.
.
.
.
.
C...**..+**...*.**....**.....*.**.***.*..**.......**............*....*..
C.
__----_____--_
l
c+
-
.
.
.
.
c*
C.
C....
1.
INFUT
HPGCIMENTS
IJ ;
R 15 E -
-
c.
C.
l
NBEG, NEND
c.
c.
C.
C.
c..
CO
C.
C.
C.
c.
TWO PRRHMETERS
SPECIFYING
FRCTORIZRTION
RRNGE
HI EXTENDING
FRON EQUHTION
INBEG+l>
THROUGH
NEND.
c.CEE DI’CU:‘:
ION) . CONPLETE
FRCTOPIZRTION
OF CA1
I;
qBTRINED
WITH
NBEG=O
RND NEND=N
<PIRTRIX
ORDER,.
IRDEP
LD
,N+Is-UURD
HPRRY
ISEE
DICClJ;SION,
Y
H
N
DOTPRD
FUNCTION
RETURNING
INNER
(DOT,
SCHLRP)
PRODUCT
UF TWO VECTOYS.
THIS
FUNCTION
IS
INVOKED
BY
DOTPRD(P,Q.N,
AND SHOULD
EVRLURTE
THE INNER
PRODUCT
OF THE TWO N-WORD
VECTORS
P ,?ND Q.
SINGRB
LOGICRL
FLRG CONTROLLING
SINGULRRITY
TEST
(5)
IS
MeTHI?:
CR,
l
C+
OF
DIAGONRL
SCRRTCH
FLOATING-POINT
IF
ICOND
= 0, 4.N
IF
LOCRTION
POINTERS
RPRRY OF LENGTH
ICOND
IS POSITIVE.
AT
LERST
C.
C.
C+
CO
c.
C.
C.
C.
CO
SINGRB
SINGHB
=
=
.TRUE.
.FRLSE.
c.
ACTION
VERIFIED
TO
:
BE
TRKEN
RBORT
SETTING
IFAIL
= J
SET D(J)
TO Y.EPSNRC.R<A
PROCEED.
.
.
.
.
.
.
IF
THE
co
C*
C*
C.
c.
.
FIND
FOR NORMAL USE OF SKYFRC
IN EQURTION
SOLVING,
ONE
SHOULD
ALWAYS SET SINGRB
= .TRUE..
THE .FRLSE.
SETTING
IS RECOMMENDED
ONLY IF THE DECOMPOSITION
(3)
IS TO BE USED
IN THE CRLCULATION
OF EIGENVECTORS
OF [RFFI
BY INVEPSE
POWER ITERATION.
.
.
.
.
+
LOGICAL
FLFlG
CONSTRAINED)
DEFINITENESS,
.
.
INDICRTING
WHETHER
CR1 <OR CRFF,
IS TO BE CHECKED
FOR POSITIVE
I.E.,
CONDITION
16)
TESTED
:
IF
C*
c*
C.
.
.
.
.
l
PDCHEK
C*
C.
C.
C.
C.
C.
.
l
C.
CO
C.
.
.
.
.
l
c*
CO
.
.
N
OF
.
.
.
.
PDCHEK
PDCHEK
ICOND
=
=
.TRUE.
FRLSE.
ABORT
DON’T
SETTING
IFAIL
= -J
TEST
FOR INDEFINITENESS
IF POSITIVE,
A LOWER BOUND ESTXMRTE
OF THE EUCLIDEFIN CONDITION
NUMBER C(A)
OF CR3 <OR CAFF,
IF
CONSTRRINED,
WILL
BE CALCULATED
BY PERFORMING
ONE
ITERATIVE
REFINEMENT
CYCLE
ON ICOND
‘RRNDO,,’
RHS
VECTORS.
THE LRRGER
ICOND,
THE BETTER
THE ESTIIIRTE..
HOWEVER,
FOR PRRCTICAL
R PRIORI
RSSESSMENT
OF
SOLUTION
RCCURRCY.
ICOND
- 2 IS OFTEN
SUFFICIENT.
l
.
l
.
:
.
.
.
.
.
.
C. A. FELIPPA
CO
C.
C.
IAFILE
IF NONZERO,
COPIES
OF IA,
IJILL
BE WRITTEN
IN LOGICRL
‘OUTPUT
FILES’
%ELOWj.
IF
WRITTEN.
H NONZERO
IRFILE
FOLLOWING
TWO CASES
:
1:.
C.
C.
c*
iH>
,B,
E:
fh!HRNING
PHPTIRL
c*
.
l
.
.
l
.
.
.
.
.
ICOND
IS POSITIVE
<iiF,,
REQUESTED,
ITERRTIVE
PEFINENENT
OF SOLUTION
VECTOR:
WILL
BE REQUESTED
WHEN CRLLING
6lrYSOL.
C.
c*
RND ITS
FHCTDRIZATION
UNIT
IAFILE
(SEE
ZERO.
NO COPIES
ARE
IS REQUIRED
IN THE
- DO NOT USE
FRCTORIZATION
FERTUREC
OPTION
IS
<RI 01;
BEING
iB>
IF THE
EXERCISED.
.
.
.
CO
c+.rr
‘2.
INPUT-OUTPUT
HPGUNENTS
l
-
CO
c.
C.
C.
C.
C.
Cl...
C.
C.
H
3.
HPPRY
CONTRININ‘
INPUT
MATRIX
lJNFACTORIZED
PORTION
THEREOF
OUTPUT
FRCTOPIZHTION-DEFINITIDN
DESCRIBED
IN THE DICCUSZION.
OUTPUT
RPGUMENTS
TO BE FACTORIZED
IF NBEG ST 0,
ELENENTC,
HZ
l
.
.
l
-
.
RCOND
EPTIMRTE
C.
OF
CIA>
IF
ICOND
GT
0
(RND
IFRIL
=
0’
.
DETi=
IS PREi;T-TO
I:<1 &;NPUT.
IF DETLF
I;
ZERO
ON INPUT,
THE DETEPMINHNT
IS NOT CHLCULHTED.
IF OUTPUT
IFAIL
= +,I
DETCF
IZ
ZET TO ZEPO.
IF OClTPlJT
IFAIL
= -JI
<J 6T O>,
THE DETERNINHNT
a,lbi
I’
THRT OF THE UPPER
*J r: i’
PPINCIPHL
MINOR.
i.
C.
C.
C.
NE6E I5
COCiNT OF NESATIVE
[HI
I;
CONSTPHINED,
EIGENVALCIES
IF
IFRIL
=
OF [R,
m.OR tF,FF,
0 ON OCITPUT.
IF
c.
1: l
FACTOPIZHTION
IFHIL
1:.
C.
C.
1:.
C.
FRILUPE
IFRIL=
IFRIL=
IFRIL=-J
5.
CON”ON
INPUT-OUTPlUT5
NONE.
c....
n.
COMNON
OUTPLlTl
NOllE.
0
J
FHCTOPIZRTION
CUCCESSFULLY
EXECUTED
:INGVLFIPITY
CONDITION
r51
DETECTED
INDEFINlTE
CONDITION
(6,
DETECTED
ICOMNON
.
l
-
.
l
.
.
7.
1PlFlJT
FILES
.
NONE.
l
l
:3.
OUTFIUT
FILES
.
-
l
IHFILE
IF NONZERO.
LOGICAL
UNIT
IAFILE
RECEIVES
COPIES
OF THE ORIGINHL
AND OUTPUT
CONTENTS
OF RRRRY A,
SR’VED F,; 2 FORTRAN
BINRPY
RECORDS
OF SIZE
LDTN+l).
‘3.
XRRTCH
FROGRHM
l
:UBPOUTINES
DOTFPD.
.
XYRD”,
l
I~Y?OL
.
l
:UBPOLITINE;
:ORT,
FOPTRHN
BINRRY
.
1’0
l
T
ir’E
DI”EN:ION
HND
:
T H
T
E ”
NOTE - RPPROPIATE
‘JHLClEl
IIF
CDC 6000,‘i’“O”
SERIES
IB”
36W320
EERIES
IBN
360,370
SERIES
UNIYHC
1108~1110,
IBFl 7094
DRTR
EPCMHC
c
c
i
,‘I.
DETCF,
DOTFPD,
EPPNRC
FOR VARIOUS
:.llE-15
‘TINGLE
9.54E-07
+EAL.4
2.2X-16
cRERL.8
,.49E-08
CIINGLE
SAVE
=
=
=
‘;
4?E-8.’
0
0
0
qPIGINHL
MATRIX
NldH = IRBS<LDcN+l,
IF
I IRFILE.EO.
0,
IF
‘:NBEG.NE.
0,
REWIND
IRFILE
WRITE
(I-FILE)
IF
IRFILE
NE
0 AND
NBEG
I
GO TO
GO TO
cR~J,rJ=I.NWFO
200
200
EF5MRC
COMPUTERS
PRECISION)
PRECISION)
PRECISION,
PRECISIONr
INITIRLIZHTION
IFAIL
IDETE:I
NEGEIG
E N T
11
V
fIbI>,
HIJ?,
DI
LDf.1)
ZINGRB
PDCHEK,
IHIJ~,II,DMAII
REAL
REHL
INTEGER
LOGICRL
EOUIVRLENCE
c
c
c
.
l
IO.
i
c
.
.
.
.
FILES
HONE.
LIRRHRY
C....
11.
c*
c.
C,,r...+...rr.*....*.*........****+.*..
c
.
.
l
c*
$.
c
c
c
.
.
INPUTS
NONE.
c.
c
c
.
l
INDICRTOP:
l
C.
C.
C.
C.1..
C.
C.
.
.
.
.
.
.
4.
I:*.*.
.
l
.
l
C....
c*
C.
C***.
C.
,- l
C.
C.
C....
C.
C.
C....
C.
.
.
.
r,,P
=
0
RRE
19
Solution of linear equations with skyline-stored symmetric matrix
I;
c
COMPUTE
200
SQUHRED
LENGTHS
OF
UNCONSTRRINED
NBEGPl
= NBEG + 1
00 IO”0
I = NEEGPlrN
II
= LOiI+l>
IF
<III
V<I>
= HrIIa**2
m=
II
-I
K = MHXO <NBEGPl.IRBS~LO<I>~-M+1)
L = “IN0
INENO,I)
- ,
IF
<t-L>
400
00
500
600
8 00
1000
800
J = IrL
IF
<LD iJ+l)
)
HIJ2
= HsM*J1..2
V’I,
= ‘$‘I>
+
fiIJ2
V<Jl
FzIJ2
=
V<J,
+
RDWS
NBEG+l
THR”
1000,1000,400
500.500,1000
800~800~600
CONTINUE
tONlINUE
c
FHCTOPIZHTIUN
SECTION
1:
00
c
C
c
4000
J
=
NbEGPlrNEND
COMPUTE
CU SUPEROIRGONF1L
IF UNCON‘:TRRINEO
ENTRIES
OF
J-T,,
COLU,,,,
OF
U,,
c
1200
JJ = LDfJ+Ia
IF
<JJ)
D = HB’J,,
4000.4000.1200
JMJ = IHB:rLD<J,
>
Jli = JJ - JMJ
hU = JK 1
IF
tCU.EQ.
0)
DO 2000
E = 1,y.u
I = J - JY + K
V ,‘b, > = 0.0
GO TO
II
= LD’I+l)
IF
(III
,I = “IN0
III-IHBS,LDrI,a,k,
I,
= JMJ + c
V’IO
= HOI,>
- DOTPRD
H<IJ,
=
V*k.>.RtII>
CONTINUE
LB00
2000
-
2200
2000~2000.1800
,
IH~II-M>,V’C-llj,fl)
c
c
c
COMPUTE
0
=
OIRGONRL
D TJ)
ELEMENT
D -
OOTPPD
SINGULRPITY
TEST
<AtJMJ+I),“.KU,
c
c
c
TOLROW = 8. +EPSMFIC*SQRl
IF
CABS (0).
GT. TOLROW)
IF
<SINGfiB>
0 = TOLROW
R<JJ,
= 1.0/D
2200
2500
r” <J> >
GO TO 2500
GO TO 6000
C
c
UPDRTE
DETERIIINANT
IF
OETCF
IS
NONZERO
C
-3000
3200
3400
C
C
C
0. 0)
IF
<OETCF.EQ.
OETCF
= DETCF*O
IF
<ABS (OETCF)
.LT.
1.0)
DETCF
= DETCF*.
0625
IDETEX
= IDETEX
+ 4
GO TO 3200
IF
(FIBS (DETCF)
. GE.
.062?5)
DETCF
= DETCF.16.
IOETEX
= IDETEX
- 4
GO TO 3400
POSITIVE
DEFINITNESS
IF
<Q.GT.
0.0)
NEGEIG
= NEGEIG
IF
(POCHEK>
CONTI~~UE
3500
4000
C
C
r
SRVE
C
C
C
3500
GO TO
3400
GO TO
3500
POCHEK
IF
IRFILE
NE
07 (J,
GO TO
6500
0 AND
SET
A RANDOM
I
UP
CRLL
=
N -
NEND
GO TO
GO TO
5000
5000
GO TO
5000
ESTIHRTION
CONDITION
0.0
.TRUE.)
4000
,J==,,NWR>
0)
K=O
nwax
=
DO 4200
=
GO TO
1
0)
MATRIX
C
C
C
+
FFICTOPIZRTION
IF
fIRFILE.EQ.
IF
<NEND.NE.NI
WRITE
<IAFILE>
IF
(ICONO.EQ.
(IF
CHECK
GO TO
1,ICOND
SKYRD”
DHS
<VI
VECTOR
NWA.
0.0,
IN
V
K)
C
c
SOLVE
OBTRIN
C
C
.
4200
c”
C
CALL
SKYSOL
1,
IRFILE,
DMRX = AMRXl
CONTINUE
C
RETURN
=
DMAX/
CYCLE
OF
SOLUTIOW
(A,
Nr LD,
DOTPRD,
V<2*N+l>,
DELTA)
IDELTA,OMRX>
NOW ESTIMRTE’C<@)
RCONO
5000
FOR V AND DC, ONE
2-NOR(1
RELRTIVE
FROM
< (1.
LARGEST
O+D"AX,
.EPSMRC)
ITERATIVE
ERROR
IN
01
OELTR
-1,
RND
REFINEMENT
DELTA
VI
TO
V<N+l>r
MRCHINE
PRECISION
N
C. A. FELIPPA
20
1: .
C.
1: .
I‘.
I: .
z.
4: l
4: .
C.
c*
C.
C.
c*
c.
,: .
1: .
c.
c.
c.
WHERE
[BF,
RND [BC,
RYE PPE:CRIBED
FtND THE FREE-FFEE
POPTION
b3FFl
HH:
BEEN PPEFHCTOPED
INTO
CLFFI
[DFFI
[IJFFI
81’ ;kYFHC.
IIrKE
THE CONiTFHINEI!
CH’E
ENBOI1IEI
THE UNCONSTEHINED
ONE.
THE
FROZE::
DESCFlPTION
BELUbl A:‘,LI,,ED
H i,,NiTRHINED
‘IITEN.
.
l
1: .
C.
C.
C.
CO
THE
IOLUTION
PROCESS
IN
ZC’iSOL
CON51STS
OF
UP
TO
:I):
STHGE?
-
FHS NODIFICATION
:
tBF1
IF REPLRCED
THIS
STEP
IS SKIPPED
IF THE LINEFIR
OR HONOGENEOUSLY
CONSTRAINED
<[XC1
BY [BFI
CRFCI
[XC,.
SY?TEM
IS UNCONSTRAINED
= CO,,.
2.
FOYWRPU PEDUCTION:
I? COLVEP
FOP
[ZFl
C.
3.
SCRLING
4.
BACKSUBSTITUTIOfl
IS SOLVED
FOR
THE
TRIFlNGULRP
SYSTEN
[LFF,
[ZF,
=
CBF,
.
.
.
l
l
:
THE
DIAGONRL
SYSTEN
tDFF1
CYF,
=
KZF,
IS
SOLVED
.
.
C.
:
THE
TDIANGULRR
SYSTEM
LUFF1
C%FI
=
CYFl
CAFI
CO
CO
C.
C.
C.
.
l
C.
C.
C.
c.
c.
.
l
1.
c+
c.
c+
l
l
C.
C.
C+
l
.
t
.
.
l
5.
ITERRTIVE
REFINEMENT
: IF HRGUMENT
IFEF
GT 0,
THE RESlDUHL
[RF1
= CBF,
[HFFI
tXF1
[HFCI
[XC,
IS COMPUTED
USING
HIGHEF
PRECISION
HPITHMETIC,
FlNIt STRGES
2-4
FEPERTED
TO
%ILVE
[AFF,
CDXF,
= CRFl
lNOTE
THAT
CRC1 = tOI,.
THE
COFRECTIQN
CDXFI
I;
HDnER
TO [XFI
RND THE PROCESS
FEPEHTED
IREF
TIMES.
.
c*
C.
6.
CONSTRAINED
FHS RECOVERY:
CBCI = CACFI
CXFI
+ [KC,
[XC,
C.
IS CRLL‘ULATED
IF REBUESTED
<SEE INPUT
APGUNENT
1810
CO
C....+.*..+.*rr............*...*........*.*...*..*.*.*.....**...........
-__-__---_---C.
CO
CI;HGE----_-_--_-_-C.
c.
C....
1.
INPUT
RR‘UNENTS
C.
c*
R
OUTPUT
MRTRIX
FACTORS
PRODUCED
BY SKYFRC
C.
C.
N
SYSTEM
ORDER
C.
C.
IN+I,-WORD
RFDRY OF DIRGONAL
LOCRTION
POINTERS
LD
C.
ENTRY
POINT
OF INNERPRODUCT
FUNCTION
KF,
SCYFHC’
C.
DOTPPD
C.
OPERRTION
CONTRDL
FLAG
:
C.
IOP
C.
c*
IOP = 0
: NOR”RL
SOLUTION
PROCESS
C.
IOP =-I
: EXECUTE
STRGES
1-2
RND EXIT
(USEFUL
c*
IN CERTRIN
SUBSTRUCTURING
APPLICATIONS~
C.
IOP = 1
: EXECUTE
STFlGES 3-4
AND EXIT
(USEFUL
C.
ON STARTING
RN INVERSE
ITERATION
PROCESS,.
C.
c+
FLHG CONTROLLING
(IUTPUT
CONFIGURRTION
OF ARRRY B.
IBX
CO
C.
IBX
= 0 : B RND X RRE IDENTIFIED
IN THE CALLING
C.
PROGRRN.
THE SOLUTION
VECTOR
X THEN OVERWRITES
8.
.
l
.
l
l
.
.
.
l
l
.
.
+
.
l
.
.
l
l
.
.
l
L
.
.
+
Solution of linear equations with skyline-stored symmetric matrix
c*
C*
C*
c*
C*
c*
C*
C*
21
THIS OPTION MRY BE USED ONLY IF
0% THE SYSTEM IS UNCONSTRRINED.
OR HOMOGENEOUSLY
CONSTRAINED WITW NO RECOVERY OF CBCI DESIRED
<S> NO ITEPRTIVE
PEFINERENT
REQUESTED <IREF=U>
.
.
.
IN TWE
IBX = -1 : B RND X RRE NOT IDENTIFIED
CRLLINCI PROGRRM BUT CONFUTRTION
OF EBCI IS ND1
DESIRED.
l
l
l
c+
LB)! = 1 : B AND X FlPE NOT IDENTIFIED
CHLLING PROGRAIl RND RECOVERY OF CBCI
c,
c*
IN
IS
TWE
REQUIRED.
l
l
.
l
l
.
C.
C1
NOTE - IF THE INPUT SYSTEM CONTRINS NONHOK16WEOUS
I.ONITPAINTC,
IBX NUST RE -1 OR +lr
I.E.,
ARRRYS
Et AND ,< MAY NDT SHFIRE THE SANE STORRSE SPACE.
C.
c.
c.
CC
CI
c*
c.
c*
c*
IF GT 0, PERFDR” IREF CYCLES OF ITERATXVE
REFINE”ENT
OF THE SOLUTION VECTOR.
IN PRRCTICE
!,NE SHOULD NOT REQUEST IREF GT 1 Uht_ESS THE
ENTRCE; OF ‘RI HRPPEN TO BE EXRCT NUNBERS
IPEF
.
.
.
.
.
l
.
l
l
;:
i.
C+
CC
f.
C....
IHFILE
LOGICAL UNIT
WCTORI’HTION
v
.;[email protected] FLOHTIN6-POINT
ARRAY OF LENGTH 2*N
TO BE U-ED IN ITERHTIVE
REFINERENT
PROCESS.
NRY BE A DUMMY RRSUMENT IF IREF = 0.
2.
INPiJT-OUTPLIT
CONTRINING
COPIES OF CR1 AND ITS
<SEE TKYFACI.
REQUIRED IF IREF GT 0
AFlGUMENT’
-
.
.
.
.
.
.
.
l
c.
l
c+
ON INPUT,
HRRHY E CONTAINS PRESCRIBED COMPONENTS
.
OF CB1 HND IX,,
I.E.,
IF THE I-T”
EWURTION IS FREE l
rCONSTRRINED> I EtI,
‘ZTORES THE RHS <LNS) VRLUE.
THE OUTPUT CONFIGURRTION
OF B IS CONTROLLED BY
INPUT FLHS 16): <SEE INPUT ARGUMENTS>.
B
c.
C+
C.
1:.
c*
;:,.+
3.
.
.
.
.
.
OUTPUT RRGUMENTS -
l
c.
x
CDMPUTED X?LUTION
C.
C+
VECTOR
DELTFI
IF IREF 67 0, RRTIQ OF LRST CURRECTION LENGTH
TO FII,AL
XLUTION
LENGTH. NOT UiED IF fRfF = 0.
l
l
0.
1%.
C....
C.
1:l
Ct...
c+
.
l
.
4.
COMFION INPUTS NONE.
5.
COMNON INPUT-OUTPUTS
NONE.
.
.
.
c*
-
.
l
.
C...+
C-+
c+
C....
c+
c*
c.
C.I..
C.
6.
CONNON [IUT~UTS NONE.
7.
INPUT
l
l
.
FILE- ., -
l
.
IHFILE
3EE IUBROUTINE
;Fr’FRC.
l
.
8.
JUTPUT
z:...?.
C.
FILES
NONE.
.
.
.
SCRATCH FILES
NOI+E.
.
.
C.
l
C.1..
CC
10.
PRO&RR” SUBROUTINEi
DOTPRD, X’YMUL
II.
LIBRHR’i’
.
l
l
C.
Cr...
c.
c+
JUHROUTINE6
:DPT,
FORTRAN BINHRY
.
.
f
c
.
I/O
T
7
P E
D I
HND
N E N ? I
I
J N
T H T E M E ‘l T ”
C
INTEGER
YEHL
DEHL
EQUIVALENCE
C
C
C
LD’l)
Bll>r
RI11 I
‘i < 1 > ,
PNOWN,
rBI.:xI,:wORN~
BXFRC,
Xl,
XI,
x*.1,,
INITIHLI~ATION
150
100
HS”IGN
3100 TO NEXT
a
hREF = 0
B.eHC = *.*
IF IIBX.EQ.
Ot
BXFAf = 1.0
DO 150
I = ItN
x <I ) = B < 1 )
IF iIOP.GT.0,
1F <IBX.EQ. 0)
C
C
C
RHP
300
500
6(30
=
TO 280
130 TO 1800
50 TO 1100
F I
C H T I
DO 1000
I = 1.N
II = LDfI+l>
IF III>
BI = BfI,
IF <DI.EQ. 0.03
II
400
N 0 D I
al
3*0*
0 N
1OOOt 1000
GO TO 1000
-11
K = I - II + IABS<LD<I,>
DO 900
J = KIN
JJ = LD(J+li
IF <JJ>
M=J-I
IF tN>
x1,,
= xc.0
- A~II+M>*BI
6D TO 900
IJ = JJ - N
+ 1
900.Y0~,400
50”,600,d00
DELTR
XNORM
22
C. A. BLIPPA
IF
‘IJ-IHBSrLD(JI~>
X<J,
= xc,>
-
8UO
900
900,900,*00
ATIJ>.BI
CONTINUE
CONTINUE
1000
C
C
C
FORWHRD
1100
1200
1300
1500
5
U
E
1500
1 = IrN
II
= Lncx+1,
IF
III,
XTI,
= 0.0
GO TO 1500
IMI
= IRBS(LD<I,,
M = II
IMI
- 1
XII>
= Xii,
- DOTPRD
CONTINUE
IF
cIOP.NE.0)
5 T
I
taoo,
no
2000
2000
II
=
X<I>
c
,;
:
GO TI1
II I =
GO
IF <“.EO.
0,
2500
2800
2’500
:(.I-J’
J
=
Z T
I
T 1, T
I
J
PHZ3
22uo.2’BOo.240,:,
I Ii
I
-
= t,n
<*.I-,,
-
Gil TO
L2800
GO TO
4000
+?,,*-_,,.:?:,I,
CONTINUE
CONTINUE
IF
sIB:X*IHFILE.EQ.
0,
i31OI>.
3500,
GO TO NEXT,
c
I
I ‘I
T E P H T
c
3100
E,
I =I - I
3000
c
SOOCl
P-SC
1, B
I = PI
DO 3000
i = IrN
II
= ia*r+u
1F ,,I,,
.A
BXFAC.B
TD 8:300
M = II
- IHBSacI,,
no
1200.1300
I = IrN
IfiBxrLntt+i,,
= A<II,*YcI.I
B H 1:’ /,
2400
PHI’
N
~RII~I+l,rX,,-~‘,M)
TCHLING
220n
1 0
DO
C
C
C
1800
T CI T
E:REF = lPEF
IF
I~REF-IR~F,
+
P E F
E
I
II
E M E PI T
;
E C
T
I
0
N
1
32,,0,3500.4u0l3
i
c
CALCULATE
RESIDUAL
WJ
RETURNS
IN .?r
I‘
‘VECTOR
HNn OLD
[RI
[Xl
= [Bl
IN V
-
[HI
[Xl
LI:INl:
SKYMUL.
c
3200
C
c
C
,.
REWIND
IHFILE
NW
= IFIBS rLn n,N+l>
REHn
<InFILE)
CRLL
iKY”UL
‘.A,
N,
)
Ln,
X.
~,H’Jl,J=,rNWHI
‘<, -I,
8,
50LVE
FOP CORRECTlO,,
CI1:*:,
OLD :,,LUTION
RNn EVHLIJHTE
I
3500
RERn
IHFILE)
B%FAC = 0.0
HSEI‘N
3500
TO NEXT
GO TO 1100
XNORN = 0.0
PNORN = 0.0
no 3800
RNORM + i~.Il.XcI,
lirl>
= V,‘I)
+ .<<I>
XNOP”
= SNOP”
+ .S<I,.X~ll
CONTINUE
DELTA
= 0.0
IF
<XMORM.NE.
0.0)
IF
<KREF-41
sUHICH
2-NDR”
‘“‘)
FlPPEHRZ
RELHTIVE
IN Y’r
ERROF
CORRECT
UELTA
*‘HIJ~,J=~,NWH)
I = 1,N
PrlOPM
=
3300
C
C
,-
CaN;TPHINEn
-4000
4200
IF
no
PH3
<IBY.LE.O~
4800
I = 1.N
II
= LncI+i>
IF
~11)
IMI
= IRBs~Ln~I)~
M = -II
- 1111 DCITPRD
,HllMI+,brX,r-n,,n,
DO 4600
J = IIN
4400
4600
4soo
5000
C
R E C 0
V E R Y
GO TLI 5000
4200.4200r4S00
1
B<I, =
C
C
C
C
C
C
C
C
C
C
DELTA
= CORT <RNORM, XNORPI)
3100r3100~4000
I,
= IRBS’LD~J+~;~
+
IF
IIJ-IRB~TLD~JI,,
B<,>
= BrIl
+ H<IJl.ll(J>
CONTINUE
CONTINUE
RETURN
END
I -
SUBROUTINE
Ln,
Sk.‘fWL
<HI
N,
J
4600.460074400
:“:,
RX,
IOP,
B,
R)
THIS
SUBPOUTINE
POSTWLTIPLIES
SKYLINE-STORED
MATRIX
[RI
BY
VECTOR
[Xl.
THE RESULT
[WI
IS PRODUCED
IN DOUBLE
PRECISION.
FURTHER
OPERRTIONS
RRE CONTROLLED
BY INPUT
FLAG
IOP
:
IDP0
EXIT
AFTER
OBTRINING
[RX,
FORM RESIDURL
VECTOR
[RI
= LB,
[RX,
IN R
IoF=
I
W,VE
[Xl
TO RI RND FORM LB,
[AX,
,N X.
IOP = -1
HRRRYS
AX
DOUBLE
INTEGER
RERL
PRECISION
(n.P.,
RNn
Lnci)
HII,
I
R MHY
BE
EQUIVRLENT
IPI
THE
HXil>,
HIJl
FIX,
*cl,,
XII>,
i?(i)
CALLING
PROGRR”
23
Solution of linear equations with skyline-stored symmetric matrix
DO
x =
~000
1-N
II
= IHB’ILD‘I+lI,
i
HXtI>
= DBLEcH~II,,*DRLEiXrI,I
tl = II
- IABScLDrI>>
IF
(ll.EU.
0)
no
is00
GO TO
2000
K = lrrl
J=I-k
AIJ
= R\IX-K.I
,%<*.I,
= AX(I)
R%<J,
= RXrJ,
CONTINUE
CONTINUE
1500
~000
1
+
+
RIJ*DBLE0ZIJ>)
HIJ*DBLEiX*I>,
r
TF
ZOO
Wl
!?-‘I,
XCI,
300
3500
3300
5000
2500.5000.3500
**[lp>
1 = 1rN
DO 2800
iixtr,
=
=
=
%<I,
m_E<B<l~~
IF
,LDII+l)
CONTINUE
GO TO 5000
-
.LE.
I = f,N
DC) 3300
P%Il
= *BLE<Brl>,
IF
~LDII+I>.‘E.O~
CCINTINUE
KETUFN
END
SUBRCOJTINE
HXI
0)
-
SCYRDA
‘VI
x,1>
= 0.
R<I>
=
Fixc1>
V(IERN,
Ns
a.
IffUR)
c
c
c
c
c
c
THIS
RaUIINE
FILLS
OUT THE N-WORD
FLEIRTING-POINT
RRRHY V WXTH
PSEUDO-RRNDOM
NUPlBERS UNIFORMLY
DISTRIBUTED
IN THE RANGE
FLRG
IRDM
NUST
BE
SET
TO
ZERO
BEFORE
THE
IVMERN-.S.VNERN+.S).
FlRST
CRLL
TO SCYRD,,
TO INITIRLIZE
THE RFlNDO,, NUHRER
SEQUENCE.
E
:
,,OTE - RNY ‘RERSONABLE’
RRNDOR NUNBER
GENERRTOR
WILL
DO FOR
THIS
RPPLICATION.
THE ONLY REQUIREPIENT
IS
THRT THE FIPST
TWO SIGRlFlCRNT
DIGITS
PRSS
THE EYEBALL
TEST.
R RRTHER
SFIALL
,,ODULE
120.16
= 65536)
IS USED HERE 50 AS TO PRODUCE
IDENTICFIL
SEQUENCES
IN ALL COflPUTERS
WITH
32-B,T
OR LRRGER
WORD SIZE
C
C
300
500
REAL
V(N)
C = VMEAN - 0.500
IF
<IRDS.NE.
0)
ix
= 26983
IRwl
= 1
DO 500
I = 1,N
IX = MOD <97181+IX+13849,6S536>
x = IX
‘)(I>
= .15258?89E-4*X
+ C
CONTINUE
RETURN
END
GO TD
300
C...*.**.*.*.....*r,+rr+r,r**....*.******..**~...***...*..*..***.***.****..*
c+
TEST PRDGRR”
SKYTEST
Cr+*.*.+*.**.**rr..*..****.**..*.***.**~*...*.*~*********.*.*..***..***.
FEiP
SIYFRC
C*
C.
C*
A
C.
ISKY
[AI-SKYLINE
TYPE
INDEX
? = FIXED-BRND,
3 FULL
c*
N
3YSTEN
CO
C.
LEVC
RPPROXIMATE
NUMBER
OF CONSTRRINED
INTEGER
PERCENT
<NO CUNSTRRXNTS
TEST
Ci
CRSE
CR3
[Xl
=
tB1
15
DEFINED
BY
HND
SIX
*
SKYSOL
INTEGER
: I = ‘RHNDOfl’
TRIANGLE
PRRRflETEPS:
SKYLINE,
l
C.
c.
c.
C*
CO
C+
C‘
CO
ICOND
C.
c.
ORDER
:
2
THPU
NflH:<
(SEE
DHTR
IF
STHTEnENTSl
EQURTIONS
IN
LEVC = 0,
A POSITIVE
ICOND
REQUESTS
COMPUTHTION
OF RN ESTIHATE
OF THE SPECTRAL
CLlNDITION
NUMBER OF CR> USING
ItDND
‘RRNDDM’
VECTZIRS
<SEE SUBRnUTINE
SKYFAC>
IDU
R POSITIVE
REFINEMENT
IF NUNZERO,
CA1 WILL
BE
IREF
WHEN
PEQUESTS
SOLVING
THE EXRCT
PRINTED.
IREF
CYCLES
FOR
[Xl
(SEE
RND
COMPLITEB
OF ITERRTIVE
SUBROUTINE
FACTOPtZRTION
IHE
SIX
PARRtlETERS
IISKY,N,LEVC.iCOND,IREFIIIiUi
CARD INPUT
FILE
5 IN
<6X4>
FOR,,RT.
ONE CRPD PER
DATR DECK SHOULD
BE TERMINATED
BY R BLRNK
CARD.
:CYSO,.I
OF
c.
C.
C.
C*
C,
c.
c.
CO
C.
C*
C
C
C
C
C
E
C
C
l
t
*
*
.
l
*
*
ARE RERD FPllM
TEST.
THE TEST
C.
c.
*
.
.
.
*
*
l
IPEF
c.
c.
C.
CI
C*
C*
c
*
.
*
.
SKYTEST
CRLLS
FOUR RUXILIRRY
ROUTINES:
;KYLDU,
SKYMWP.
SkYOPC
RND SKYPRT,
WHICH
RRE LISTED
AFTER
SKYTEST.
SAWPLE
INNEIPRDDUCT
ILlUTINES
WHDSE ENTRY
POINT
NAHES
<VIP%.
VIP%,>
RPE FRSSED
IN
ARGUMENT
LISTS,
RRE FiLSD LISTED.
SKYTEST
ALSO
INVOKES
THE TIMING
FUNCTION
CPTIME..
X = CPTIMEiO2
INITIRLIZES
THE CPU CLOCK.
X = CPTINE<I)
RETURNS
CPU TI”E
IN
FLORTlNG:POINT
SECONDS
ELAPSED
SINCE
THE LF,ST CPTIME<O)
CHLL.
THE USER SHOULD
INSERT
APPRUPIRTE
REFERENCES
TO NIS
SYPTEN
CLOCK RWTINES
IF HE DESIRES
TO TIllE
SKYFRC.
l
.
t
*
*
+
*
l
.
.
l
*
*
l
*
IF TESTING
ON CDC 6OOOf7000
SERIES,
INSERT
THE FOLLOWIPIG
PRUGRRM SKYTEST
<INPUT.OUTPUT~TFIPES=INPUTITAPE6=OUTPUT,T~PE2~
CRRD
DI%NSIONS
OF TEST
RRRRYS
RRE R(NWANAX>,
DUiNWAMHX,,
B<NMRXk,
BEX <NHRX> , DBEX <Ni,RX> , LD <NNRX+ 1) I X <Ni,RX) , XEX rN,,(RX> I Y <J.NllFtXI
WHERE N”RX
= LRRGEST
N, NWARRX = LRRGEST
CR, STURRrjE
ALLDCRTIDN
CONCERTING
NATRIX
DISPLRY
PARRHETERS
IPPT.
IUPLOW,
SEE SKYPRT
FOR FACT[fRIZRTION
ARGURENTS
IAFILE,
PDCHEK.
SINGRB,
SEE SKYFRC
C. A. FELIPPA
24
DOUBLE
PRECISION
DBEX<,OO)
EXTERNRL
VIPSS,
VIPSD
INTEGER
BLRNK,
CMHRK,
LD<lOl>
LOGICAL
PDCHEC,
SINGHB
RERL
R(2000,.
DU<2000,
RERL
B<lOO>,
BEX’lOO>.
‘V<JOO~,
DATA
BLANK
,‘1H ,I
C”RRK
~IH.,,
DRTR
IPRT
,2x.
IUPLOW
2’2,
DATA
NMAX /100/r
NWRNRX ~2000~
DATR
PDCHEK
..FRLSE.,,
SINGAB
X<lOOl,
IRFILE
J2,
XEX<lOO,
,‘. TRUE.
/
C
C
C
RERD
100
TEST
PARRMETERS
RERD
<5r10>
IF
(N.LE.
0)
WRITE
t6,20>
IF
<N.GT.NNAX>
IRDN
= 0
FROM
CARD
FILE
ISK’l,
N,
LEVC.
ISKY,
N,
LEVC,
<UNIT
5)
ICOND,
STOP
ICOND,
GO TO
IREF,
IDU
IREF.
5200
IDU
S Y S
T E M
C
C
C
C
TEST
BUILD
CONSTRUCT
DIAGONHL
LOCATION
POINTER
RRRRY
IN
LD
c
C
400
450
C
500
550
C
600
650
1000
C
C
C
C
LD(,>
= 0
600,500.400
,F
(ISKY-2)
.
SYMMETRIC
MATRIX
ISKY=
: FULL
DO 450
I = LIN
LD<I+l>
= LD<I)
+ I
GO TO 1000
ISKY=
: FIXED-BRND
NATRIX
(HALFWIDTH
M = MAXO<BI
M = MAX0<3,2.N15)
DO 550
I = 1111
LD<I+,,
= LD<I,
+ NINO<I,M)
50 TO 1000
ISKY
= 1 : ‘Rf,NDOM’
SKYLINE
CALL SKYRDM
<VI N,
0.5,
IRDM>
C = O.G*FLORT(N)
DO 650
1 = 1,N
K = RNRXI (C.VII>,
1.000*>
LD<I+l>
= LD<Ij
+ NINO<I,KI
NW,? = IRBSrLD<N+l))
GO TO 5400
IF
<NWA.GT.NWAMRX)
MHRK HPPROXIMATELY
‘RANDOIILY’
ON LD
1200
C
C
C
IF
<LEVC.EQ.
0)
CRLL SKYRDM
(VI
Nr
DO 1200
I = 1rN
IF
<V(I).LT.O.O>
CONTINUE
SENEPRTE
1500
1600
IF
TEZT
CD1
LEVC.NllOO
LEVC
IS
CONSTRAINTS
NONZERO
GO TO
Ol*FLORT(LEVC)
0.5-O.
1500
I IRDFO
LD<I+,)
HND
[ill
O.J*N))
IN
RRRRY
CRLL ‘ShYHDM
<DC,, NUH,
0.0,
IRUM,
DO 1600
I = l,N
II
= IRBS<LDcI+l)r
DLlrIIl
= ,C..DU<II,
IF
<IDU.EQ.
0)
WRITE
(6725,
CALL SYYPRT
!DU,
N, LDr
IJPLOW,
=
-LD(I+l>
DU
GO TO
2000
CHF12.6,,
IPPT)
C
RSSEMBLE
MHTRlX
CR,
=
CL1
CD1
GUI
r
2000
c
c
CRLL
WRITE
CRLL
WRITE
CRLL
CRLL
SKYLDU
(6.30)
SKYPRT
c6,32.,
SKYMRP
SKYOPC
GENERATE
[BEX,
=
C
C
CRLL
CRLL
rDU,
R,
N,
LD,
Vr
‘VIPSD)
<A.
N.
LD,
IUPLOW,
6HF12.6.r
tR,
(0,
Np
NI
LD,
LD,
IUPLOW,
KOUNT,
IPRT,
RRNDOM LHS VECTOR
IN XEX HND
IF!,
[XEX,
IN DOUBLE
PRECISION
PKYRDN
SKYNUL
<XEX,
Nt
0.0,
IRDM)
<RI N, LD,
XEX,
DBEX,
0,
IPRT)
RSSOCIHTE
DUNMY,
RHS
VECTOR
DUMMY)
C
C
S V ‘i
F
H C
TElT
r
DETCF
= 1.0
DUNNY = CPTI”EtOl
CALL SYlFAC
CR, 0, N, N, LD,
V,
VIPSS,
SINGRB,
PDCHEK,
.
ICOND.
IHFILE,
RCOND,
DETCF,
IDETEX,
NEGEIG.
IFRIL,
TOP = 1. E6rCPTIME
111 /FLORT
<KOUNT)
WRITE
q6.50,
SINSRB,
PDCHEK,
IFRIL.
DETCF.
IDETEX,
NEGEIG
IF
<IFHIL.NE.
01
GO TO 100
WRITE
r6,54,
KOUNT,
TOP
IF
<ICOND.NE.
0)
WRITE
~6.56)
HCOND
IF
<IDU.EQ.
0)
50 TO 2200
lWPITE
(6,60>
CRLL ‘KYPRT
(H.
N, LDI
IUPLOW.
6HFl2.69,
IPRTl
c
c
SKY:OL
C
C
C
3ET
2200
2500
C
UP
SINGLE
PRECISION
DO 2500
I = 1,N
BEXrI,
= DBEX(I>
B<I,
= BEX(I>
IF
rLD<I+l>.LE.O)
CONTINUE
TEST
RHS
VECTOR
B(I)
I:
=
XEX<I>
Solution of linear equations with skyline-stored symmetric matrix
SOLVELINERR
c
c
CALL
SYSTEM
SKYSOLw,,
N,
WRITE
(6170,
ENllPM
= 0.0
XNORM = 0.0
DO 3200
I = 1.N
XNOR”
3 XNORM +
XERR = XEX*I>
-
RND
LD,
PRINT
25
RESULTS
VIPSS,
0,
I,
B,
X,
IREF,
IAFILE,
V,
DELTR)
Xll,.X<I,
,-:<I>
ENORM = ENODM + XERP.XEPR
BERR = BEXcI)
- BtI’
MhRk = BLRNK
IF
<LDrI+l,.LE.O>
MRRK = CMhPh
IWRITE
(6.72,
IrMHRbr
BEXrI)rPcI*.BERP,
XEXrI~,S,I~,AERP
CONTINUE
:XERP = SQRTcENOP”+XNOPM,
WPITE
t6.741
.<ERR
IF
rIREF.EQ.
1)
ldRITE
t6r76,
IlELTR
GO TO 100
3200
C
C
C
DIHGNOSTICS
WRITE
GO TO
WRITE
GO TO
5200
5400
C
1:
C
10
FOR
ILLEGFlL
t6,82>
IO0
(6.86)
100
TEST
PI,
NUH,
FORflAT
SThTEMENrS
FORMAT
(614,
PHHAMETERS
NMAX
NulAMhX
20
.
.
llH.
C-LEVEL
= 13.
22H
O-O,
lCOND,!REF,ID”
= I2rlH,.I,,,H,,I,~
5x 67<1H-)
)
25 FOPPIRT
<, ‘5XZ”TEST
CDINV,
HND CL11 : I
30 FORnAT
!,‘%3lHTEZT
MHTRIX
CR, = CL,
CD,
[“,
: >
32 FORllRT
~I’~;(~OHSC*LINE
MHP OF CRI
:I
50 FOHRRT
i,’ ‘SA27HSKYFRC
EXECUTED
,JITH
IINGRB
= ~3. Ilk,,
PDCHEb
= L3,
. SXPHIFRIL
= 14,16H,
DETERMINANT
= FlS.:.dH
. 2..
14,
5X22HNEGHTIVE
EIGENVRLUES
= 14,
54*FORllhT
~SX17HOPERRTION
UNlTS
=I9/
. 5X20HhVG
TI”E,OP
l,NIT
=
F8.2r9H
MICROSEC)
56 FORFlhT
t’5A 34HF(hTRIX
CONDITION
NUNBER
ESTIMRTE
= lPG12.3)
i0
FORMRT
rl4XlSHCOMPUTEP
FRCTOPS
:)
70 FORMAT
r/15:<43HShYSOL
TEST
RESULTS
‘0 = CONZTRHINED
EQS)
: ‘,
. 2X 2HEQ BXTHESRCT
B
5X7HCUNF.
B 3:-: 7HB-ERROR
4:i 7HEXhCT
N
.
5?: 7HCOMP.
.X
3X TN’<-ERROR
a
72 FOPMAT
~‘l:~rI3rRlr2F12.3,,PE~O.l~~~PFll.~,Fl2.7,~~E,O.1~
74 FORMRT
< ‘5A3SHHCTUHL
RELRTIVE
ERROR OF COMPUTED
[X3
1PE12.2,
i6
FORMRT
< 5X3YHESTIAATED
REL.
ERROR qIF UNREFINEII
[Xl
lPE12.2)
82 FORMRT
(,14H
. . . ORDER N = I4r15H
EXCEED’
NNRX =
14,
86 FOR,,hT
r25H
. . . H-IIRTRIX
WORDS NWH = 15r17H
EXCEEDS
NYRFlRX = 15)
END
;VBPOUTINE
c
C
C
C
C
SCYLDLl
#HFRC,
R,
N,
Ln,
‘i’r
DOTPRD,
GIVEN
MHTRIX
FRCTORS
CDINVI
RND [Lll
‘TORED
IN hRPRY
BFRC,
THIZ
iUBROUTINE
CONPCITE:
[HI
= [L,
CD1 CIJlr
WHERE
CL1 IS
THE TRHNSFOSE
OF [Ul
REHL
INTEGER
RF,%<,,,
LDO,l,
Htllr
‘~*I>
c
DO
4000
J = IrN
JFlJ = IHB~~.LU<JI>
JJ = LDnJ+I,
IF
(JJ.LE.
0)
JK = JJ - JMJ
V~JKI
= I.OIHFHC<JJ~
TO
?OOO
DO
2000
<so0
3000
3200
4000
2500
K = 1,JK
IJ = J”J
+ C
I = J + IJ
- JJ
II
= LD(I+I>
IF
(II.LE.0)
IF
<I.NE.
J,
,, = “IN0
<II-IHB”<LDiI)>,K)
F,<IJ>
= V<K>
+ DOTPRD
GO TO 2500
hcIJ,
= FIFt3CrIJ>
‘v’<K , = 0.0
CONTINUE
GO TO 4000
JJ = -JJ
II
= JMJ + 1
DO 3200
IJ
= II.JJ
ha:IJ,
- RFRC<IJ,
CONTINUE
RETURN
END
GO
SUBROUTINE
C
C
C
SKYOPC
<NBEG,
GO TO 2OOO
VIK)
= RFhC(IJ)/tiFRC(II>
1
rhFRC~~I-n,,V(K-fl,rn)
NEND,
LDr
COUNT)
RETURNS
NUMBER OF qPERRTION
UNITS
<KOUNT>
REQUIRED
FRCTOR
h SKYLINE-STORED
SYMMETRIC
MhTRIX
FROM ROW
(NBEG+i,
THRU NEND
(EHCH
N-WORD
INNERPRODUCT
IS
COUNTED
AS M OPEPRTION
UNITS>
c
INTEGER
LD(l>
KOUNT
= 0
NBEGPI
= NBEG + 1
DO 4000
J = NBEGPl,NEND
JJ = LD<J+L)
IF
(JJI
TO
26
C. A. FELIPPA
IuBPOUrINE
c
c
4:
c
c
C
c
c
Ir\NHP
*h,
IGENEFHTE”
OF’DEP N.
h TENPLHTE
:‘>‘“1I:DLi +s-,
IllFLOG
IUPLOLI
ItE’.;‘ICE
DEVICE
1
2
1
2
=
=
=
=
REAL
INTEGER
INTEGER
DHTA
:
:
:
:
N.
LD,
FIHP’
OF
0 DENOTE
PPlNT
“HP OF
PRINT
FlHP OF
FOPMHT
PPINT
FOPMHT
PRINT
IUPLOl,,r
DEVICE>
PIRTPIX
CA1 OF
ENTRIES
R LD-1‘k.‘fLlNE-:TOYED
PO’ITIVE,NEGhTIVE,:EfO
CIPPEF TPIRNGLE
LOGlEP TRIHNGLE
FOP 72-CHHPHCTEP
TELETYPE
F[lR
128+
CHHRRCTEP
LINE
PPINTEF
ht1,
RLhNb 7 CFLHG,
C”AEk,
PlJU160~
La’l!,
BLHNk,CmARC,POSIT,tiEGRTI:EPO
DEVICE,
/lH
NEGRT,
, I!,..
POSIT,
lH+r
I,+r
ZERO
IHO,
c
KCMHX = MRX0~30,‘3OrDEVICE>
JREF
- 0
JBEG = JPEF
+ 1
JENR - HINO
<JPEF+YC,,AX,N>
kC = JEMD - JREF
JREFPS
= JREF + 5
WRITE
<6.30,
(J. J=JREFPSr
,END,S>
1,7X12,10)
30 FQRRRT
DO 1300
J = JBEG.JEriD
I = J - JREF
ROW\Il
= BLRNC
IF
(LD(J+l).LE.O,
ROWrl)
= CMARK
COPlTINUE
1300
WRITE
(6,401
tROW<J,rJ=,.bC,
<8X 60R2,
40 FORMHT
IBEG
= 1
IEND
= JEND
IF
<IUPLOW.EQ.
1)
GO TO ,500
IBEG
= JBEG
IEND
= N
1500 DO 4000
I = IBEG.IEND
DO 1600
J = lrKC
1600
ROW tJ>
= BLHNC
CFLRG
= BLhNh
II
= LDiI+l>
IF
(II.LE.0)
CFLHG
= CMHYI
J1 = I
J2 = JENIl
IF
<IUPLOW.EQ.
11
GO TO ,800
II
= IhBS\II>
J1 = I + IhBS(LD(I>’
- I,
+ 1
Jl
= MAX0
rJl.JBEG>
J2 = MIPlO
<l,JEMD)
IF
lJ.?-J,,
4000,
*Boo,
lSO0
DO 3000
J = JlrJ2
1800
k = J - JREF
M=J-I
I,
= II
+ II
IF
<IUPLOW.NE.
1)
GO TO 2000
IJ = IRBS<LD(J+l)a
- M
IF
!IJ.LE.IABSrLD~J>>>
GO TO 3000
IF
(h<lJ,>
2000
2200.2400,&00
ROW<K,
= NEGRT
2200
GO TO 3000
ROW ib.)
= ZERO
2400
GO TO 3000
YOWtkl
= POSIT
3000
CUNT 1 NLIE
WRITE
t6.60)
I,
CFLHG,
<ROW (h> I K=, I KC>
4000
CONTINUE
60 FORMRT
<lX14.‘,3.60R.?>
JREF
= JEND
IF
(JREF-II)
1100.5000,5000
5000
RETURN
1000
1100
?W.FOIJTINE
IhYPPT
‘H.
II,
LD,
IsIPLOU.
F”T,
IiEVlCEl
R OF
OPDEF
c
PRINT
LD-:b,‘LINE-ITOPEU
iPHPSE
“RTPIX
N.
:
c
c
c
ILlPLObl
IUPLOCl
DE?‘ICE
“E’rVICE
L
=
=
=
=
1
2
1
2
:
:
:
:
PRINT
IIPPEP
FRINT
LDIJER
FOPb,AT PRINT
FOPFlRT PRINT
TRIANGLE
‘C’ILINE
TRIRNOGLE 5k YLINE
TELETYPE
(TI”E:HRFE
W[IRb )
FOP :2-CHHP
FOR 128+
CHRP LINE
PRINTER
‘BHTCH
RUN,
c
c
c
c
MATPI:<
ENTPIE:
HPE
OF I,IIUTH
I2
(E.G.
INTEnSER
INTEGEF’
INTEGER
DHTH
DhTa
IlhT”
.
c
PPlNTED
Ei2.4,
H‘COPDINS
TO E,F
OF
Ph;FEU
THPU fiRGI_I”ENT
6 FOR”HT
FMT
Solution
P:Cr!HX
lo00
C
c
C
IlOO
1200
1300
40
c
C
C
=
SET
COL
c
symmetric
matrix
21
15, S.DEVICE,
PRINT
PHNGE
ROW PRINT
IBEG
with skyline-stored
0
idBEG
JBEG = JREF
+ f
JEND = MINO
<[email protected]+KCMHX,NI
kC = JEND - JF’EF
_I = JBE‘,JEND
DO 1300
I = J - JPEF
IDCONS< I,
= BLANC
IF
rLD<J+i)>
IDCONS
= CMHYC
CONTINUE
<C(ILLRB,J,
WRITE
<6,40,
FOYM(RT
),‘14X
10’43?14rR1,4X>)
;ET
1500
mwo
=
JREF
of linear equations
RRNGE
TO
JEND)
RHD
LIST
,200rl2OO,
IDCONS<J-JREF>
(IBEG
THF’U
COL
IDENTIFIERS
1300
, J=JBEG,
JEND)
IEND,
= 1
IEND
= JEND
IF
IIUPLOW.EQ.
L>
IBEG
= JBEG
IEND
= N
Do 3000
IBEGrlEND
KR = 0
IF
~IUPLOW.EQ.2I
DO 1600
c = f.KC
PPEPHRE
I-TH
ROlrl OF UPPER
TRIRtiGLE
RUWF,,TrK+3,
= LKIPI~)
J = JREF
+ C
M=J-I
IF
cM.LT.0,
IJ
= IfiBS’LD’J+l))
- M
IF
IIJ.LE.IHBS<LD<J,,,
POWF,,T*h+3,
= FMT(1,
kR = hP +
ROUcKR,
= Ro:I.J,
CONTINUE
60
TO
1500
50
TO
2000
GO TO
lliO0
GO TO
1600
I =
1
,600
5
IF
‘kH,
PREPRRE
I-TH
ROW OF LOWER
II
= IHBS<LDrI+,>I
M =
II
- IAB?<LD~~I~~
,I
= MRAO <I-WrJBEG,
J2 = “IN0
#:I, JEND)
L.R = J2 - J1 + 1
IF
II;P>
r = 11 - .lBEG
2 100
ROWFNT~JI
= TTIPth+l’
POWFMTtS,
= 4H
20
ROWFMTa.G>
= FpITc,)
IJ
= II
Jl
DO 2200
K = lrrR
P,,W<r,
= RaI.J>
3000,3000.2500
TPIHNGLE
2000
1
300o,3000,2l00
I +
2200
i
c
I?
LIST
=
IJ
I-TH
1
+
POW
C
;‘s*i,
TFLHG
=
RLHNC
iF-‘;D~x;il.Ll.O>
3000
4000
5000
POWLHB,I,CFLHG,
WF’ITE
IG,ROWFMT,
CONTINUE
JPEF
= JEND
IF
(JPEF-PI,
PETURN
END
RERL
FUNCTION
SFLAG
= CUARK
<ROW<J,rJ=1.CY,
,10*,5000.5000
VIPPE
<HI
Bl
Na
c
C
c
C
C
:HMPLE
SINGLE
SHOULD
200
500
5000
PEAL
‘VIPS:
=
IF
tN)
DO 500
VIPSS
ON
fi
GIVEN
SKYSOL.
MRCHINE)
B<il
fir,,,
0.0
5000,5000,200
I
=
= l.N
‘VIPSC
FUNCTION
SA,,PLE
DOUBLE
SHOULD
200
500
1000
5000
SERVING
SKYFAC
RND
(RSSEMBLY
LRMGURGE
+
A~I>.B~I~
PETURN
ENU
RERL
C
c
C
C
h
INNERPRODUCT
FUHCTION
PROCEDUPE
PPECISION
ACCUflClLATION
VERSION
BE U?ED
IN HCTUHL
IMPLEME~TRTION
VIPTD
CR.
B,
PO
INNERPPODUCT
FUNCTION
PROCEDURE
PRECISION
RCCU,,ULHTION
VERSION
BE USED
IN ACTURL
IMPLEMENTATION
RERL
A<l,,
B(l)
DOUBLE
PRECISION
UM
SU”
= O.ODO
IF
(PI>
DO 500
I = IrN
SUM = SU”
+ DBLE(AtIl)rDBLEtB(I))
VIPSD
= sun
RETURN
END
SERVING
SKYFAC
RND
(ASSEMBLY
LANGUAGE
ON 1 GIVEN
,,RCHlNE,
1000.1000.200
SKYSOL.
28
C. A. FELIPPA
...
SKY-TYPE
TEST
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
1,
NATRIX
CR1
N
=
=
CL,
COL
044013
-:061730
-.052174
-.326797
.720432
.053372
.319443
-.008196
COL
6
.174217
.029194
-.031293
-.359695
-.412109
-.039510
.085142
16.
17
18
19*
20
SKYLINE
+--++----+-------+-+++
13
14
15
16.
17
18
19,
20
COL
COL
COL
.007034
.410675
COL
16*
-2.709717
-. 282257
-. 184418
-.248322
-.236160
[A,
COL
9.
5.821777
166977
:375198
-. 000900
.238645
- .044449
-. 149384
115341
-: 287247
13
COL
1.048912
- .34 0589
-. 042077
.480682
.250496
17
.689583
-. 157687
-.416275
.048782
14
COL
lOI
-3.020508
.478058
-.258881
-.318405
047058
-:182816
.300110
-.441345
COL
15
-.230179
.373713
018601
-:361099
-. 053287
-. 005294
-.279541
CDL
COL
13
.269398
.439880
-. 058580
19*
-6.405029
.254974
-.050121
172791
:016881
052201
:006973
-.0x4154
COL
20
-.17989,
:
‘0
15
1 I>
. *
.
l
+
+_+__-_-+
+-++--++-+---+++-++---++++---+-+-+
-+--+
+--_-++--+--+__+_+-
SYYSOL
19.
20
COL
-. 090900
-. 156782
.417618
-.200029
OF
1.032168
110855
-: 443243
-.417660
-.411774
-. 013899
.040987
-. 277634
.446950
-.302399
12
-. 061659
.043387
101807
-: 098354
MRP
3
.230813
-.042492
-.089312
-.016849
-.240814
.208542
-.093933
-.054065
.473480
-.011119
.494019
;;;;tC=EXECUTED
WITH
CIN‘AB
01
DETERMINHNT
=
NEGATIVE
EIGENVALUES
=
5
226
OPERATION
UNITS
=
34.50
AVG
TIME.jOP
UNIT
=
13
14
15
16.
17
18
-. 382484
-.015320
-. A46509
7
5
__+____++
12
8
.611516
.212633
COL
+
+--__+
-+++-+-++---++
-+--+++
10.
9.
LO.
11
12
4
COL
-. 151699
067160
-: 036890
.105057
.362091
-. 145353
-. 077434
COL
O.l,O
+-+-
11
:
6
7
3
COL
.213440
.026407
.365494
11
.079424
.144593
=
+
++-+
4
5
6
3
3
9.
EQ
1
2
3
ICOND.IREF,IDU
O/O,
-. 077390
119363
: 182684
-. 298889
-. 196823
-. 209992
-.122696
-.023942
5
1
2
3
2
-.320265
082017
:329712
-.333603
.272368
-. 035472
.025108
COL
ROW
ROW
ROW
ROW
ROW
30
[Lb, :
-1.770384
-. 177671
-. 135824
11
12
I3
14
15
16+
17
18
19.
20
=
.040587
9r
lo*
11
12
14
15
160
17
19.
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
CD,
...
EXRMPLE
C-LEVEL
COL
1
.125876
018694
:017722
.013656
1
2
3
4
5
6
7
8
94
10.
11
12
16,
17
6
7
8
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ROW
ORDER
TEST
20,
TEST
RESULTS
EXHCT
B
-.0368776
-.3354552
-.0437254
-.3132310
1016244
-:0281424
.3474827
-.3593653
-2.2627583
.3993893
-.4837782
-.0035529
.2332043
-.1591169
.0594882
1132414
-:3659491
.1939173
-.9196062
-.0221439
HCTUR!_
RELRTlVE
ESTI"ATED
REL.
(.
=
ERROR
ERROR
T, PDCHEK
-.1393258
=
.
F
2..
-24
MICRDSEC
CONSTRAINED
CmlP.
I:
-.0368776
-.3354552
-.0437254
-.3192.310
.I016244
-.0281424
.3474*17
3593653
-2.2627583
.9993393
-.4837782
-.0035529
.>33.1049
-.1591169
.0594882
.1132414
-.3659491
.1939178
-.9196062
-.0221439
-.
=
B-ERROR
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
3.0-08
-7.5-09
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.5-03
0.0
OF COMPUTED
OF UNPEFINED
EP;I
:
EXifCT :i
-.010'3253
.2503308
-.2339832
.4102020
-.0374023
-.4717565
-.4175110
-.1548004
-.4187622
-.1642303
.2662653
-.4160319
4168701
:1573742
1092224
-:0142975
-.4036304
.P2?2023
1582336
:1593579
@.I
[Xl
t:ONP. X
-.oio9253
.2503308
-.2389832
.4102020
-.0874023
-.4717865
-.4175110
-.1548004
-.4187622
-.1642303
e66e659
-.4160919
.41CJ701
.1573742
.I092224
-. 0142975
-.4086304
.2292023
1.24-03
1.'37_07
.1582336
. 1593573
X-ERROR
-7.2-09
0.0
0.0
3.7-09
3.7-09
7.5-09
3.7-09
-5.6-09
0.0
0.0
a. 0
-3.7-03
0.0
-5.6-09
4.7-09
0.0
0.0
1.9-09
0.0
1.9-09
Solution of linear equations with skyline-stored symmetric matrix
and should be removed, unless extended double precision
arithmetic (e.g. REAL*16 on IBM 370 computers) is
available for the residual calculation in subroutine
SKYMUL.
A slightly more efficient implementation of these
procedures is possible if array LD is used to pack the
following information:
case, the routines presented here can be made use of as
follows:
1. Use SKYFAC to factor Ag upon marking the tat
equations associated with the external freedoms x, as
constrained.
2. Solve the homogen~usly constrained block system
(a) diagonal location pointer
(b) local skyline width excluding diagonal
(c) constraint flag
(20)
into 3 bit fields. Such a version is not presented here since
partial word manipulation with Fortran statements is
highly machine dependent.
APPLICATIONS
General
The equation solving procedures presented here may be
useful in the following application areas:
1. As fast linear equation solvers for small-scale
systems (normally not exceeding 1000equations). This is
typical of most finite element codes written in universities
and research institutions, as well as pilot pro~ms
developed in industry for special applications or feasibility studies of nonlinear or time-dependent problems.
2. As freedom “condensation” processors for iirstlevel superelement analysis. This application is elaborated
in the following subsection.
3. As a conceptual aid in the design and coding stages
of the development of large-scale equation solvers based
on the skyline storage approach.
Superelement condensation
Superelement te~hniques[I4] (also called ~subst~cturing’ or ‘tearing’ by structural engineers and ‘dissection’
or ‘blocking’ by numerical analysts[9]) are modelpa~itioning procedures intended for a multistage analysis
of complex finite element idealizations. A block of finite
elements modeling a simply connected portion of a
complex system is called a level one superelement. A
connected set of level one superelements makes a level two
superelement and so forth until the complete model is
realized. The fundamental characteristic of the associated
matrix assembly process is that ‘internal’ degrees of
freedom of all superelements of given level (those
freedoms not shared by two or more superelements) are
reduced out or ‘condensed’ before proceeding to the next
level.
The governing lineart equations for a particular
superelement may be represented by equation (4), in
which subvector xf collects all internal freedoms whereas
xc includes external or connected freedoms. Elimination
of x, yields:
A*x, = b”
(17)
where the (symmetric) condensed superelement matrix
A* and right hand side vector b* are given by
A* = A, - A$A;’ A,
(18)
b* = b, - Af:Ai’bf
(19)
provided Aff is nonsingular. The matrix condensation
process (18) can usually be carried out in core for level
one (and sometimes level two) su~reiements. If such is the
~Superelementtechniquesare mosteffectivein linearproblems.
id..
Vol. 5;No.
I-C
29
by submitting to SKYSOL (with argument IBX = 1) nc
right hand sides consisting of the n, columns of A&
completed with zeros on the constraint equation positions.
3. The condensed matrix A* is
A”=A,,-A;X=A,,-3
(21)
so that one should simply substract each column of B (as
they are being put out by SKYSOL) from the approp~ate
column of A,, to produce A* in situ (overwriting A,,).
Note that X is not used.
The reader may verify that the entire process can be
carried out without prearranging A, and that only two
auxiliary moating-lint vectors of length n itre required in
addition to the A-definition arrays A and LD. The
reduction (19) of the load vector can be executed in a
similar manner.
REFERENCES
1. C. A. Felippa, Finite element and finite difference energy
techniaues for the numerical solution of martial differential
equatidns, Proc. 1973 Summer Comp~fer-S~m~~~~onCffnf.,
Montreal, Canada, 1-14 (July 1973).
2. R. A. Willouahbv, A survev of snarse matrix technoloev. in
Tech. Report A&DL-TR-jl-79, ‘Air Force Flight Dyn&&s
Lab., W~ght-Pa~te~on AFB, Ohio, 65-171 (June 1971).
3. .I. K. Reid (ed.), Large Sparse Sets of Linear Equations. pp.
17-24, 41-56, 97-104, 255-277, Academic Press, London,
(1971).
4. A. George, A survey of sparse matrix methods in the direct
solution of linear equations, Proc. 1973 Summer Computer
Sjmu~uf~o~Co& Montreal, Canada, U-20 (July 1973).
5. A. Jennings, A compact storage scheme for the solution of
symmetric simultaneous equations, Comput. J. 9, 281-285
(Sept. 1969).
6. R. J. Melosh and R. M. Bamford, Efficient solution of
load-deflection equations, J. Struct. Div. ASCE, Proc. Paper
I\io. 6510, 661-676 (April 1969).
7. B. M. Irons, A frontal solution program for finite element
analysis, Int. J. Num. Math. Engrg. 2, 5-30 (January 1970).
8. R. P. Tewarson, Sparse ~airjces. Academic Press, New York
(1973).
9. 0. Birkhoff and A. George, Elimination by nested
dissection, in Complexity of Sequential and Parallel Nume~caf Algotithms. J. F. Traub (Ed.), pp. 221-270, Academic
Press, New York, (1973).
10. J. 0. Owens, The influence of machine organization on
algorithms, in vol. quoted in[9], 11l-130 (1973).
11. J. H. Wilkinson and C. Reinsch (eds.), Handbook for
Automatic Computatjoa Voi. IX-Linear Algebra Springer,
Berlin (1971).
12. C. Johnson, On the convergence of a mixed finite element for
analysis of plate bending, Num. bats. 21, 43-42 (1973).
13. 5. H. Wilkinson, Rounding Errors in Algebraic Processes.
Prentice-Hall, London, 1973.
14. P. 0. Araldsen, the application of the superelement method in
analysis and design of ship structures and machinery
components, presented at the National Symposium on
Computerized Anaiysis and Design, George Washington
University, Washington, D.C. (May 1972).