Comments
Transcript
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).