Home

SPARSPAK: Waterloo Sparse Matrix Package User's Guide for

1. Program 2 C SPARSPAK A ANSI FORTRAN RELEASE III NAME EXSB o C UNIVERSITY OF WATERLOO JANUARY 1984 CRA REE REET b kk k t kk bkd bkd db k dakk RARER RRA RARA EE CT E EE TAKTA ek MA 1 NL IN E PROGRA M TELE DR COFRE AAA bb b bkd k AAA RA AAA d dd k dd tva RRE AA RARA RRA CARRERA ERA AAA REA ii C o FILE REQUIREMENT z o UNIT 3 FOR SAVEA AND RSTRTA SHOULD NOT BE o DESTROYED C CREE RARER k bk bb tk bb C o INTEGER I TERRA IPRNTE IPRNTS MAXINT MAXSA e e MSGLVA NVARS REAL MCHEPS RATIOL RATIOS TIME REAL S 999 REAL ERROR FOUR ONE TWO ZERO i TERRA RARA November 1984 1 il 49 SPARSPAK A c COMMON SPAUSR MSGLVA IERRA MAXSA NVARS COMMON SPKSVS IPRNTE IPRNTS MAXINT RATIOS RATIOL MCHEPS TIME C HRARRRARRECRDCORARDARAORARORODURRDORHDAENANAODAODHAADANSORHSBHANA BANABAS QQQQ qQaQQ QQQ Q QQQ A 00002 AAAQ AAAA 100 200 300 400 a e z e e s e e e s ea a e a a e a e e vw s e e a we CALL SPRSPK MAXSA 999 CALL IJBEGN o DO 100 I 2 200 CALL INIJ I 1 1 CONTINUE CALL IJEND IF IERRA EQ 0 GO TO 00 CALL STATSA STOP CONTINUE w a 4 m e we e e e e e e e e e e e e e gw e we s e CALL ORDRAS S IF IERRA EQ GO TO 300 CALL SAVEA 8 S CALL STATSA STOP CONTINUE ZERO 0 0E0
2. 94SPK UNIVERSITY OF WATERLOO 4 4 SPARSE MATRIX PACKAGE RARA AA SPARSPAK ERA AAA RELEASE 3 EXRHRAERRRASR C JANUARY 1984 eee eeeeeee ANSI FORTRAN SINGLE PRECISION eeeeeeeeee LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES OUTPUT UNIT FOR STATISTICS _IJBEGN BEGIN STRUCTURE INPUT INIJ INPUT OF ADJACENCY PAIRS IJEND END OF STRUCTURE INPUT eet ORDRAS NESTED DISSECTION ORDERING EMSGA SYSTEM A ERROR ORDRXI X A B AND Im1 2 3 4 5 6 ERROR NUMBER I NSUFF STORAGE FOR ORDERING MAXSA MUST AT LEAST BE SAVEA SAVE STORAGE VECTOR STATSA SYSTEM A STATISTICS SIZE OF STORAGE ARRAY MAXSA NUMBER OF EQUATIONS NUMBER OF OFF DIAGONAL NONZEROS TOTAL TIME REQUIRED MAXIMUM STORAGE REQUIRED Program 3 C C C UNIVERSITY OF WATERLOO 124 1400 999 200 998 SPARSPAK A ANSI FORTRAN RELEASE III NAME EXSC 1SPK JANUARY 1984 2SPK CASARSARARRARAREARRRRARAAARAMRA KARA EE SSPK November 1984 SPARSPAK A CEPT OH OOOOH EHO TO SEH ANA NANA NNA rr rr rs EE Crrerer eres MAINLINE PROGRAM Oco CHOOSE lod OOH AAA AAA AAA AAA AAA AAA AAA AA ONES OREM EEN CIPNAPNPFDRDAODGCDODRDROROBKEBORDBOSDORERODBDARARODAENNAREREAKBARHBHHBBAR DSANN c C FILE REQUIRMENT c UNIT 3 FOR SAVEA AND RSTRTA SHOULD NOT BE C DESTROYED CroerssrsarRON dr rr rro rn res SP SR Pos P PPnoroenan cr pS K XSn prrs n INTEGER I IE
3. l CALL IJBEGN DO 100 I 2 10 CALL INIJ 1 1 1 S 100 CONTINUE CALL IJEND S C o Corr iii gt wee eee eee ee eee ee C DETERMINE SYMMETRIC ORDERING NA RA CALL ORDRA2 S c C onic Se a et cer Ste Se Seca eee C INPUT NUMERICAL VALUES CG O sania a a an a enee man emm e e e IF I GT 1 CALL INAIJ2 I 1 1 ONE S IF 1 LT 10 CALL INAIJ2 I I 1 ONE S CALL INAIJ I I FOUR Mil ii November 1984 ST Te own 35 1SPK SPK 3SPK 4SPK gt 5SPK 6SPK 7SPK gt 8SPK 9SPK 10SPK gt 11SPK 12SPK 18SPK 14SPK 15SPK 16SPK 17SPK 18SPK 19SPK 20SPK 21SPK 22SPK 23SPK 24 SPK 25SPK 2 6 SPK 2 7SPK 28SPK 2 9SPK SOSPK 1SPK S2SPK 33SPK 34 SPK 5SPK 36SPK SISPK 38SPK 99SPK 40SPK 41SPK 4 2SPK 49SPK 44SPK 45SPK 46SPK gt 47SPK 48SPK 49SPK SOSPK 51SPK 52 SPK SPARSPAK A 200 09999 QQQQ QQQQQ QQAQAQ CALL INBI I FOUR CONTINUE CALL INBI 1 ONE CALL INBI 10 ONE S o e e e 46 e e 46 wg oe e e e we e we we e e e e we oe e e we we we O e e e e e e we we oe e e e we e e COMPUTE AN ESTIMATE OF THE RELATIVE ERROR IN THE COMPUTED SOLUTION e e a a e S e DC e e e e e e e m a e e e e e an v e e e ww e m m Ma e e AM 2 O e User s Guide e a e e e e we oe e we we e e we e we e we e e e ew e e e e e e e e e e e we e e e e e e
4. next table ORDRxi ORDR S S Matrix input INADI INAIJi 1 J VALUE S INROWi I NIR IR VALUES S INMAT NIJ II JJ VALUES S INBI I VALUE S INBIBI NI II VALUES S INRHS RHS S oa Factorization and or Solution Solution SOLVEt SOLVER S S BREST RELERR S saa Save and Restart the computation and Restart the computation SAVEA AALS K S F RSTRTA K S Ordering Choices Right hand side input Reverse Cuthill McKee ordering 7 symmetric A Reverse Cuthill McKee ordering 7 unsymmetric A One way Dissection ordering 2 symmetric A One way Dissection ordering 2 unsymmetric A Refined quotient tree ordering 3 symmetric A Refined quotient tree ordering 3 unsymmetric A Nested Dissection ordering 4 symmetric A Nested Dissection ordering 4 unsymmetric A Minimum Degree ordering 6 symmetric A Minimum Degree ordering 6 unsymmetric A November 1984 w 29 SPARSPAK A input the structure of A order and allocate storage input numerical values for a numerical relative error estimation values ge k for b au SPRSPK IJBEGN INIJ INROW INIJIJ INCLQ IJEND ORDER where 2 A B 1 1 2 3 4 5 6 INROWI INALJi INMATi INRHS INBI INBIBI SOLVE1i EREST Sketch of possible execution paths through SPARSPAK A modules November 1984 Use
5. REAL MCHEPS RATIOL RATIOS TIME REAL S 250 REAL ERROR FOUR ONE RELERR TWO 1 ZERO 8 C ISIT III TITIS TITIS SIS TISI YET YET TIT YTITTIYYTIYYYTTTYTITTTSTT C COMMON SPAUSR MSGLVA IERRA MAXSA NVARS COMMON SPKSVS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME C CtrrrrsrrrrrrssSSrrarnrrrrrnSenSSSVprrSSssrrSSnsss p Mi e EE C INITIALIZE SPARSPAK A CO di bad a kek oes CALL SPRSPK MAXSA 250 S CE madrinas C INPUT STRUCTURE O 4 EE CALL IJBEGN DO 100 I e 10 CALL INIJ I 1 1 S 100 CONTINUE CALL IJEND S ge a u m ew ew wm e m e we e e e e e e e e e a e ao mp e we e e e e we we ew e e we e e e e e v e e e e ww e s e ZERO 0 0E0 ONE 1 0E0 TWO 2 0E0 0990 anana November 1984 wa 38 1SPK 2 SPK 8 SPK 4SPK SSPK 6SPK 7SPK 8SPK 9SPK LOSPK 11SPK 12SPK 1 3SPK 14 SPK 15SPK 16SPK 17SPK 18SPK 7 19SPK 20SPK 21SPK 2 8SPK 23SPK gt 24SPK 2 5SPK 26SPK 2 7SPK 28SPK 29SPK 90SPK 1SPK 92SPK 9 3SPK 34 SPK 3 5SPK S6SPK gt 97SPK 98 SPK SISPK 40SPK 41SPK 42SPK 49SPK 445PK 45SPK 46SPK 4 7SPK 48SPK 49SPK SPARSPAK A aaaaa aaaa ARAAA AAAA 900 200 500 11 400 2994 AAGA anan FOUR 4 0E0 DO 200 Im 1 10 User s Guide IF I GT 1 CALL INAIJS I l 1 ONE S CALL INAIJ I I FOUR S CALL INBI I TWO S CONT
6. hopefully successful run is recorded In a batch oriented environment IPRNTS and IPRNTE are usually the same Note that the user and or the computer installation must ensure that the files associated with IPENTS and IPRNTE are available to the user s program before execution begins 74 Message level MSGLVA The first variable MSGLVA in the common block SPAUSR stands for message level and governs the amount of information printed by the interface modules Its default value is two and for this value a relatively small amount of summary information is printed indicating the initiation of each phase When MSGLVA is set to one by the user only fatal error messages are printed this option could be useful if SPARSPAK A is being used in the inner loop of a large computation where even summary information would generate excessive output Increasing the value of MSGLVA up to 4 provides increasingly detailed information about the computation Note that the module SPRSPK sets MSGLVA to its default value if the user wishes MSGLVA to be different from two he must reset it after SPRSPK has been called In many circumstances SPARSPAK A will be embedded in still Boies super package which models phenomena producing sparse matrix problems Messages printed by SPARSPAK A may be useless or even confusing to the ultimate users of the super package or the super package may wish to field the error conditions and perhaps take some correctiv
7. interface routines Before the actual execution of each interface routines the variable STAGE is used to check that all previous interface modules have been successfullv completed This avoids producing erroneous results due to an improper processing sequence or accidental omission of steps The content of the variable STAGE i is ae changed after a SH has been successfully executed When an error occurs during the execution of the phase the variable STAGE remains unchanged This prevents the execution of all the subsequent phases even if they are invoked by the user The variable STAGE is also used by the modules to determine whether some initialization is necessary in a module or whether part of the module has already successfully executed during a BESSER call to it November 1984 SM 62 SPARSPAK A M NA User s Guide 10 6 Storage allocation of integer and floating point arrays The ANSI FORTRAN standard specifies that the number of bits used to represent integers and floating point numbers are the same However some vendors provide the user with the option of specifying short integers either explicitly in the declarations such as INTEGER 2 or via a parameter to the FORTRAN processor which automatically represents all integers using fewer bits than used for floating point numbers Since a significant portion of the storage used in sparse matrix computations involves integer data for pointers subscripts etc it is
8. the storage pointers in the storage map block as well as the storage vector S In this way the state of the computation can be re established by executing the module RSTRTA which restores the control block and the storage map block and the storage vector S The variable MXUSED in the control block is used to avoid saving irrelevant data from S After the successful completion of each phase MXUSED is set to the maximum number of storage locations in S used thus far It is then only necessary to save the first MXUSED locations of S whenever the routine SAVEA is called November 1984 Ta A 61 SPARSPAK A User s Guide Some operating svstems allow a program to change the space it occupies in main storage during execution Thus in some installations the user of SPARSPAK A may be able to 7 _ dynamically increase or decrease the size of the working storage S He can determine what the value of MAXSA should be by declaring the labelled common block SPACON in his mainline program and examining the value of MXREQD At the end of each successfullv executed phase of the computation MXREOD is set to the minimum value of MAX SA required to successfullv execute the next phase of the computation It is often the case that when this dynamic growing of program space is provided the effect is to increase the space allocated to the unlabelled COMMON which is usually assigned the highest memory locations in the user s program a
9. 1 ERROR ABS S 1 ONE 300 CONTINUE WRITE IPRNTS 11 ERROR 11 FORMAT 10X 9 5HMAX IMUM ERROR STOP END Output re r UNIVERSITY OF WATERLOO SPARSE MATRIX PACKAGE eeeeeeee eae ee ee REE ETTIZTZZZENI November 1984 RELEASE SPARSPAK 3 33 35SPK 96SPK SISPK 88SPK S9SPK 40SPK 41SPK 42SPK 49SPK 44SPK gt 45SPK 46SPK 4 7SPK 48SPK 49SPK 50SPK 51SPK 52SPK 5SSPK 54SPK 55SPK 56SPK 57SPK 58 SPK 59SPK 60SPK 61SPK 62SPK 63SPK 64SPK 65SPK 66SPK 67SPK 68 SPK 69SPK 70SPK 71SPK 7 2SPK 78SPK 74 SPK 75SPK 76SPK 77SPK 78SPK 79SPK 8 0SPK 81SPK 82SPK 83SPK 84 SPK 85SPK 86SPK gt 87SPK SPARSPAK A k C JANUARY 1984 er rerste ANSI FORTRAN SINGLE PRECISION ee aes LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES OUTPUT UNIT FOR STATISTICS IJBEGN BEGIN STRUCTURE INPUT INIJ INPUT OF ADJACENCV PAIRS IJEND END OF STRUCTURE INPUT ORDRA1 RCM ORDERING INAIJ1 INPUT OF MATRIX COMPONENT S INBI INPUT OF RIGHT HAND SIDE SOLVE1 ENVELOPE SOLVE EREST1 ERROR ESTIMATOR STATSA SYSTEM A STATISTICS November 1984 SIZE OF STORAGE ARRAY MAXSA NUMBER OF EQUATIONS NUMBER OF OFF DIAGONAL NONZEROS TIME FOR ORDERING STORAGE FOR ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION STORAGE FOR SOLUTION TIME FOR FACTORIZATION TIME FOR SOLUTION OPER
10. 2924 0 217 2924 SIZE OF STORAGE ARRAY MAX SA NUMBER OF EQUATIONS NUMBER OF OFF DIAGONAL NONZEROS TIME FOR ORDERING STORAGE FOR ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION TOTAL TIME REQUIRED MAXIMUM STORAGE REGUIRED 53 TOSPK TISPK 72 SPK T9SPK 74SPK 75SPK 716SPK TISPK 18 SPK 79SPK 80SPK 81SPK 82SFK 83SPK 84SPK 85SPK 86SPK 87SPK SPARSPAK A Program 4 c SPARSPAK A ANSI FORTRAN RELEASE III NAME EXSD C C UNIVERSITY OF WATERLOO JANUARY 1984 Podkova idi a vd zl v EE Si a a d EE EE SET va A A ono a id o i lo aa Cees errr MA I N L I N E P R O G R A M AA E O OSO AA EE EE E EE E EE ERS LOS See eee od odolal li dla E LCL TS Eee lon d ol alt C o FILE REQUIREMENT l o UNIT 3 FOR SAVEA AND RSTRTA C TA ok RRR RRR AAA AAA c INTEGER I IERRA IPRNTE IPRNTS MAXINT MAXSA E MSGLVA NVARS REAL MCHEPS RATIOL RATIOS TIME REAL l s 82509 REAL ERROR FOUR ONE TWO ZERO C l ostases rr rr rr rro rr rr rar rro rr rr ra rra rr rro rr rr rro o COMMON SPAUSR MSGLVA IERRA MAXSA NVARS COMMON SPKSYS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME rressasesapesonpasonssnanonosasansnnanznazanosa ansansansa anaaaa QQQQ tu ch t wl dny PA ba bong N ta ta y gt a ta gt e e a e e e we oe e e e e e e s e e e we e CALL SPRSPK MAXSA m 2509 w5 a e e
11. 5 Solving many problems having the same structure sessen 19 8 Solving many problems which differ only in their right hand side 20 7 Output from SPARSPAK A assesseer i at a el 22 7 1 Message level MSGLVA coocoiconicnccicconnonos EIE E pab u tie Ma eka ee ay 22 7 2 Statistics gathering STATSA e TEETS T Ee SA A OT 23 7 3 Error messages IERRA ENEE ea Na O wi i 24 7 3 1 Save and restart routines az eeh s man a 2 73 2 Input of the matrix structure EENS 25 7 3 3 Ordering and storage allocation routines NENNEN 26 7 3 4 Input of the numerical values een EE 26 November 1984 ES io SPARSPAK A o User s Guide 7 3 5 Factorization and solution ssseseeseszzzenssnzznznnsrezzonnezzznnenzzznzsnarazonnonzz 27 7 3 6 Relative error estimation ussesssssessmsnznnnnzznnnzznnsnzzzzznnenaan nanniet 28 8 Summary listing of interface routines ses ennnnnnnzznnnnnnnnnenzznznnennnnezn 29 9 Examples EA O Stan ko 31 Example 1 0i rrrrrrreerererreorereseroreaen EE Example 2 a A e A eat BEE 35 Example 3 ENN Gees icon r in 38 Example A v n Ss Nese E CaN Ueda as te Mega ke 42 Example bin ata on s no A 47 Example 6 da A oa 57 10 Appendix implementation OVE VIEW L sessevirzazzonnnzznsenzonnsnneon esnensnensnnenensansnusuang 60 10 1 User module communication r tu a ee 60 10 2 Module module communication ssssmm
12. 6 11 S I I 1 10 11 FORMAT 10H SOLUTION 5F12 5 GAN a lw siw m m Ga m kw w m a m aa m a a a a Wa a mw m aa a a m a a e e c PRINT STATISTICS GATHERED BV SPARSPAK A e O EN dn CALL STATSA c STOP END Note If the SPARSPAK A available to you is a double precision version the REAL declaration in this example should be changed to DOUBLE PRECISION The module SPRSPK must be called before any part of the package is used Its role is to initialize some system parameters e g the logical unit numbers for output files and to set default values for options e g the message level indicator and to initialize the timing routine The routine needs only to be called once in the user program and the FORTRAN statement is simply CALL SPRSPK Note that the only variable in the common block SPAUSR that must be explicitly do a value by the user is MAX SA It is assumed that the ros which comprise SPARSPAK A have been compiled into a library and that the user can reference them from his FORTRAN program just as he references the standard FORTRAN library subroutines such as SIN COS etc Normally a user will use only a small fraction of the subroutines provided in SPARSPAK A Warning The modules of SPARSPAK A communicate ith each other through labelled common blocks whose names are SPKSYS SPAUSR SPACON SPAMAP and SPADTA Thus the user must not use labelled common blocks with these names in his program Nov
13. SPARSPAK A correspond to different El alias for choose P along with appropriate storage methods and whether or not A is symmetric When A is symmetric U is replaced by LT in the above EES and of course only one of L and L7 is stored The user chooses a particular method by calling the appropriate subroutines in Steps 2 3 4 6 and di The methods are distinguished by a numerical digit 7 1 lt lt 6 which is the last character of the subroutine names The subroutines used in Steps 1 and 5 apply to all the methods The best method to use depends very much on the particular problem and the context in which it is being solved so we cannot provide rigid rules as to which method to use Some guidelines and considerations regarding the choice of method are given in Section 3 Restrictions and assumptions 1 SPARSPAK A assumes that the nonzero structure of A is symmetric If this is not the case the package will still work but if A has highly unsymmetric structure this may lead to some inefficiencies because on matrix will be treated as Re its structure is that of GE November 1984 y 3 SPARSPAK A a User s Guide 2 SPARSPAK A assumes that for any permutation matrix P Gaussian elimination applied to PAP without row or column interchanges yields an acceptably accurate factorization LU In other words the package assumes that A can be symmetrically permuted without regard for numerical stability This is true for e
14. and then restart the calculation at exactly that point some time later To save the results of the computation done thus far the user executes the statement CALL SAVEA K S where K is the FORTRAN logical unit on which the results are to be written along with other information needed to restart the computation If execution is then terminated the state of the computation can be re established by executing the following statement CALL RSTRTA K S Examples 4 5 and 6 provided in Section 9 illustrate the use of SAVEA Sch RSTRTA Note that executing SAVEA does not destroy any information the coniput li i can proceed just as if SAVE A were not executed When errors occur in a module the routines SAVEA and RSTRTA are useful in saving the results of previously successfully executed modules see Section 7 3 and Example 5 i in Section 9 Another potential use of the SAVEA and RSTRTA modules is to make the working storage gt array S available to the user in the middle of a sparse matrix computation After SAVEA has been executed the working storage array S can be used by some other computation Finally the SAVE A and RSTRTA modules allow the user to segment the computation into several distinct phases and thereby reduce the amount of program that must be resident in storage at any given time Important Notes a In the subroutines SAVEA and RSTRTA information is either written on or read from the FORTRAN logical unit K u
15. d or e will occur upon re execution If c occurs the minimum value of MAX SA needed for d and e is provided When c or d occurs after executing SAVEA adjusting the storage declaration then executing RSTRTA one must again call ORDRzxi However the interface will detect that the gt ordering and or storage allocation have already been performed and will skip that part of the computation Note that if a user is simply using SPARSPAK A to select a particular EES c mar be an acceptable termination state See Example 6 in Section 9 2 See Section 4 for details on how to use SAVEA and RSTRTA and Examples 4 5 and 6 in Section 9 November 1984 e 1 o a NE 1 SPARSPAK A wa e User s Guide 2 4 Modules for inputting numerical values The modules in this group are similar to those for inputting the matrix structure They provide a means of transmitting the actual numerical values of the matrix problem to SPARSPAK A Since the data structures for different storage methods are different the package must have a different matrix input subroutine for each method For the user s convenience SPARSPAK A uses the same set of subroutine names for all the methods except for the last digit which distinguishes the method and the parameter lists for all the methods are the same Important Note The elements of A and b transmitted to SPARSPAK A by these routines are elther single or double precision floating
16. desirable to try to exploit these short integer features whenever it makes sense to do so SPARSPAK A contains parameters RATIOS and RATIOL set in the module SPRSPK which specify the ratios of the number of bits used for floating point numbers to the number used for short and long integers For example in a double precision IBM version of the package which exploits short integers RATIOS is 4 and RATIOL is 2 Let 1 be the smallest integer such that rx gt r The package then uses RATIOS RATIOL to allocate only p RATIOS p RATIOL elements of S for short long integer arrays of length p SPARSPAK A assumes that the declaration of S that the user makes in his program is of the same type as that used for floating point computation We also make the reasonable assumption that RATIOSPI and RATIOL gt 1 10 7 Statistics gathering SPARSPAK A contains a labelled common block called SPADTA which appears below These variables are used to provide the output described i in Section 72 COMMON JSPADTA ORDTIM ALOCTM FCTIME SLVTIM ERRTIM FCTOPS SLVOPS ERROPS ORDSTR ALOSTR SLVSTR ERRSTR OVERHD ANORM RCONDA ERRFCT RELEST SVPAD 88 In order to supply timing information SPARSPAK A assumes the existence of a real function DTIME which returns the processor execution time that has elapsed since DTIME was last referenced Thus the DTIME function is also installation dependent 4
17. e G e e e a m e e e a Ga S PV a s a a 990 anaa e a e ss nu e e a S e e e e e e e e CALL ORDRAS S IF IERRA EQ 0 GO TO 100 CALL SAVEA 3 S CALL STATSA STOP 100 CONTINUE e e e we e e a e a e e se ns e 099 a m a e e e M a oe 4 e e e e e e e s DO 200 I 1 200 S IF I GT 1 CALL INAIJS I l 1 ONE S CALL INAIJS I I FOUR S November 1984 User s Guide 1SPK SPK 9SPK 4SPK 5SPK 6SPK 7SPK 8SPK 9SPK 10SPK 11SPK 12SPK 13SPK 14SPK 15SPK 16SPK 17SPK 18SPK gt 19SPK 20SPK 2 1SPK 225PK 23SPK 24SPK 25SPK 2 6SPK 2 7SPK 28 SPK 29SPK 90SPK 31SPK 32SPK 33SPK 34 SPK 95SPK 36SPK 37SPK 38SPK 39SPK 40SPK 4 1SPK 4 2 SPK 43SPK 44SPK 45SPK 4 6SPK 47SPK 48SPK 49SPK 50SPK 51SPK 52SPK 53SPK 54SPK 55SPK 56SPK 57SPK 58SPK SPARSPAK A 200 999A AAAA an a 300 11 Output DE E E EEE EEEE KEE INBI STATSA November 1984 User s Guide CALL CONTINUE CALL INBI 1 ONE CALL INBI 200 ONE S INBI I Two S mm gt e e e e e e e we e we e e e e e we e e oe e e m a a e a a m a e e e e e e e e u m w ew e e e m m a e m e e e e e a a e e e sm a gt we e sm we we e a we a e we e m ep a e e e e we e we we e e e me e we e e e e e e m e e a e e e a ge we e e e w o m w w s m COMPUTE TH
18. e e e e e e e fF F F FF F e se e e e e e e we we COMPUTE THE ACTUAL RELATIVE ERROR IN THE COMPUTED SOLUTION SINCE THE TRUE SOLUTION IS KNOWN ERROR ZERO DO 500 I 1 200 ERROR AMAX1 ERROR ABS S I ONE CONTINUE WRITE IPRNTS 11 FORMAT 10X ERROR 3 5HMAXIMUM ERROR STOP END November 1984 1PE15 5 User s Guide 48 SISPK 32SPK 33SPK 34 SPK 35SPK 36SPK 37SPK 38SPK 99SPK 40SPK 41SPK 4 2SPK 48SPK 44SPK 45 SPK 4 6SPK 4 7SPK 48SPK 4 9SPK 50SPK 51SPK 52SPK 59SPK 54SPK 55SPK 56SPK 57SPK 58SPK 59SPK 60SPK 61SPK 6 2SPK 63SPK 64SPK 65SPK 66SPK 67SPK 68SPK 69SPK 70SPK 71SPK 72SPK 73SPK T4SPK 715 SPK 76SPK 77SPK 78SPK 79SPK 80SPK 81SPK 82SPK 8 9SPK 84SPK 85SPK 86SPK 87SPK 88SPK 89SPK 90SPK 91SPK 92SPK 93SPK 94SPK SPARSPAK A va al User s Guide Output te 444 UNIVERSITY OF WATERLOO t htH S SPARSE MATRIX PACKAGE A f G PARSPA K RARA AAA RELEASE 3 FHERRERIES C JANUARY 1984 waeeeee ete ANSI FORTRAN eeeeeeetee SINGLE PRECISION sib dal LAST UPDATE JANUARV 1984 OUTPUT UNIT FOR ERROR MESSAGES 6 OUTPUT UNIT FOR STATISTICS e IJBEGN BEGIN STRUCTURE INPUT INIJ INPUT OF ADJACENCY PAIRS IJEND END OF STRUCTURE INPUT EMSGA SYSTEM A ERROR IJEND ERROR NUMBER 116 INSUFF STORAGE FOR ADJ STRUCTURE MAXSA MUST AT LEAST BE po 999 STATSA SYSTEM A STATISTICS NO STATISTICS AVAILABLE
19. e e oe a we we e a e e e e COMPUTE THE ACTUAL RELATIVE ERROR IN THE COMPUTED SOLUT ION SINCE THE TRUE SOLUTION IS KNOWN ew we e e m e eee e e e e e e we e a we e we a we a e e we e e e e e we e rt ve e e e e e e e e ERROR ZERO DO 00 Tei 10 ERROR AMAX1 ERROR ABS S 1 ONE CONTINUE WRITE IPRNTS 11 ERROR FORMAT 10X S5HMAXIMUM ERROR STOP END 1PE15 5 s s es ests UNIVERSITY OF WATERLOO eeeeeeeee SPARSE MATRIX PACKAGE KRAK EEO SE i GP AR S PA K tornos RELEASE 3 sesssesess C JANUARY 1984 soosesssss ANSI FORTRAN 999999994 SINGLE PRECISION eeeeeeeee LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES OUTPUT UNIT FOR STATISTICS IJBEGN BEGIN STRUCTURE INPUT INIJ INPUT OF ADJACENCY PAIRS IJEND END OF STRUCTURE INPUT ORDRA2 RCM ORDERING INAIJ2 INPUT OF MATRIX COMPONENTS November 1984 53SPK 54SPK 55SPK 5 6 SPK 57SPK 58SPK 59SPK 60SPK 61SPK 62SPK 63SPK 64SPK 65SPK 66SPK 67SPK 68SPK 69SPK 70SPK 71SPK 72SPK 73SPK 74 SPK 75SPK 76SPK 7TISPK 78SPK 79SPK 80SPK 8 1SPK 8 2SPK 83SPK 84SPK 85SPK 86SPK SPARSPAK A User s Guide INBI INPUT OF RIGHT HAND SIDE SOLVE2 ENVELOPE SOLVE EREST2 ERROR ESTIMATOR STATSA SYSTEM A STATISTICS SIZE OF STORAGE ARRAY MAXSA 250 NUMBER OF EQUATIONS 10 NUMBER OF OFF DIAGONAL NONZEROS 18 TIME FOR ORDERING 0 017 S
20. format used by the later phases The user does not have to concern himself with this representation or transformation Ni Fr Notes a SPARSPAK A assumes that the value of NEQNS id m number of equations is equal to the maximum column or row index supplied by the routines which transmit the 7 7 pairs to the package Thus it is imperative that the user supplies at least one 7 7 pair for which 7 or j is equal to NEQNS The routine IJEND assigns the value of NEQNS found by the package to the corresponding variable in the common block SPAUSR b SPARSPAK A assumes that the diagonal elements of the coefficient matrix are nonzero Common Errors The most common cause of error during matrix structure input is insufficient working storage If we denote the number of off diagonal nonzeroes in the matrix by OFFDA then the minimum amount of storage necessary to EE input the structure is given by OFFDA 2XNEQNS 1 Of course sometimes the user does not know the value of OFFDA and may guess too low SPARSPAK A will still accept and count the i j pairs even after running out of storage and the user can obtain an upper bound for OFFDA by calling the module STATSA described in Section 7 after all pairs have been input The number reported may be unnecessarily large because duplicate input pairs may not now be detected and thus mav be counted more than once by the package For a complete list of errors which may be e genera
21. nei o vor 0 000008002 0008008080800000908000008000009000980900000900 O ERRAEBEBA wi M A I N L I N E P RO G R A M 44440090 SS A E A EE ll AS E lu Cees eee sr r lt lt a rrsrrnrs s SssssrsrrrVMR9I 9O ass rare nn9 lt 4M n ans pa C C FILE REQUIREMENT d C UNIT 3 FOR SAVEA SHOULD NOT BE DESTROYED C l CICOHRRADEOOCUEDBDROCCORADRCOAOORDODSODONHGDCOC ARDSDBHHADBDKDOREOBREDNDRH B6 C INTEGER I IERRA MAXSA MSGLVA NVARS REAL S 250 REAL FOUR ONE RELERR PE COMMON SPAUSR MSGLVA IERRA MAXSA NVARS C l A A AE A A ALA LE AAA A A A Ad dada Aa HAND anons dn Z UN wel CN gt te bonn N ta la gt a Le gt ns e s a e e aa we MAXSA 250 QQQQ 2 u G h3 4 S Q S s ta CALL IJBEGN DO 100 I 8 10 CALL INIJ I 1 1 S 100 CONTINUE CALL IJEND S ose a we a e e e e a we s e s e QQQQ QQQQ ONE 1 0E0 FOUR 4 0E0 DO 200 It 10 IF 1 GT 1 CALL INAIJ1 I l 1 ONE November 1984 o i 42 SPARSPAK A 200 CONTINUE a oooo QQQQ anana AAAA Output CALL mm e m e w m Mw m w e e e e e oe a m a m a m a e DA Aa e a a a w a e mm e e e e e e e e e e wm e a ep e e em e o o e e e e e e e we e SOLVE1 S CALL ma e e e e G ww e e m a am m e E e e w w e m e m am w m ep ep w m wm e e u COMPUTE AN EST
22. of the right hand side and ERESTi has to be called only once regardless of the number of right hand sides that are to be solved Furthermore ERESTi should be called only after the factorization has been do that is after SOLVE has been executed successfully November 1984 w dE 15 SPARSPAK A User s Guide 3 Some guidelines on selecting a method We mentioned in Section 1 that there are six basic methods distinguished by a numerical digit satisfying 18186 These six methods can be viewed as grouped into three odd even pairs the only distinction between method 7 odd and method 1 is that method assumes the coefficient matrix A is symmetric and method i 1 assumes A is unsymmetric Thus we really only provide three essentially distinct methods with each one having a symmetric and unsymmetric version Hence in this section we will largely confine our remarks to methods 1 3 and 5 comparative remarks about them will also apply to their unsymmetric analogues methods 2 4 and 6 The basic methods are as follows the remarks comparing them and the advice provided _ should be regarded as at best tentative Characteristics of sparse matrices vary a great deal Method 1 2 3 4 5 6 November 19 4 Basic Strategy and References The objective of these methods is to reorder A so it has a small bandwidth or profile 7 The well known reverse Cuthill McKee algorithm is used For relatively small prob
23. point numbers depending on the version of SPARSPAK A being used The examples in san manual assume a ange precision version of the pack eee is being used There are three ways of passing the numerical values to SPARSPAK A In all of them indices passed to the package always refer to those of the original given problem The user need not be concerned about the various permutations to the problem which may have occurred during the ordering step idi o When any of the three numerical input routines is first called the storage used for storing the numerical values is initialized to zero a Input of a single nonzero component The subroutine INAJJi is provided for this purpose and its calling sequence is CALL INAIJi I J VALUE where I and J are the row and column indices and VALUE is the numerical value The subroutine INAIJi adds the guantity VALUE to the appropriate current value in storage rather than making an assignment This is helpful in situations e g in some finite element M pasa don where the numerical values are obtained in an incremental fashion For example the execution of CALL INAIJ2 8 9 9 5 CALL INAIJ2 4 effectively assigns 5 5 to the 3 4 th component of A November 1984 NY 12 SPARSPAK A wa NEE User s Guide b Input of a row of nonzeroes The routine INROWi can be used to input the numerical values of a row or part of a row in the matrix Its calli
24. the factorization has been performed in a previous call to SOLVEi SPARSPAK A will automatically skip the factorization step and perform the solution step directly The solution vector is returned in the first NEQNS locations of the storage vector S If SOLVEi is called before any right hand side values are input only the factorization will be performed The solution returned will be all zeroes See Examples 3 and 4 in Section 9 2 6 Modules for estimating the relative error in the computed solution An estimate of the relative error using the infinity norm IL in the computed solution can be obtained CR executing the Tollowing FORTRAN statement CALL ERESTi RELERR S Here S is the working storage array for SPARSPAK A and RELERR is a variable which will contain the relative error estimate after the subroutine is invoked successfullv The last digit 2 O TE November 1984 1 14 SPARSPAK A User s Guide is used to distinguish between subroutines for different storage methods If the problem is too ilkconditioned with respect to the precision of the machine or unacceptable rounding error in the factorization has occurred a message will be printed depending on the value of the message level indicator MSGLVA and RELERR will be set to 1 0 The mki l based on estimates of the condition number of the coefficient oi and the error incurred in its factorization 1 Thus the estimate is independent
25. was executed and the structure was successfully input however the ORDRA5 module could not execute because MAXSA was less than 1400 The module SAVEA was then executed and the run terminated The third run had MAXSA 2500 and the ordering and storage allocation were successfully performed However ORDRAS terminated with an error because it detected that too little storage was available for the numerical computation SOLVES so SAVEA was again executed Finally the last run was executed with MAXSA set to 2509 the maximum value printed in the third run and the solution to the problem was obtained Nole The following examples were run using a single precision version of SPARSPAK A on a DEC VAX 11 780 computer The working storage required will therefore be oo if a different version of SPARSPAK A or a different computer is used Hewes 1 C SPARSPAK A ANSI FORTRAN RELEASE III NAME EXSA C C UNIVERSITV OF WATERLOO JANUARV 1984 CI k TTT CREE kk kk kk koko k kb bk bb kb CTRA MAINLINE P R OGRAM akk oo ak A ii i i i a i a a il E EE dd EE a dda dd c C FILE REQUIRMENT c UNIT 3 FOR SAVEA AND RSTRTA SHOULD NOT BE c DESTROYED c prepr C l INTEGER I IERRA IPRNTE IPRNTS MAXINT MAXSA 1 MSGLVA NVARS REAL MCHEPS RATIOL RATIOS TIME REAL S 900 REAL ERROR FOUR ONE TWO ZERO C to ika aka aa aa ka a a nl C COMMON SPAUSR MSGLVA IERRA MAXSA NVARS COMMON SPKSYS IPR
26. 017 STORAGE FOR ALLOCATION 60 STORAGE FOR SOLUTION TO TIME FOR FACTORIZATION 9 017 TIME FOR SOLUTION 0 OPERATIONS IN FACTORIZATION 13 OPERATIONS IN SOLUTION 0 TIME FOR ESTIMATING RELATIVE ERROR 9 093 OPERATIONS IN ESTIMATING REL ERROR 160 STORAGE FOR ESTIMATING REL ERROR 90 ESTIMATE OF RELATIVE ERROR 1 785e 07 TOTAL TIME REQUIRED 0 067 MAXIMUM STORAGE REQU I RED 90 SAVEA SAVE STORAGE VECTOR Program 2 C SPARSPAK A ANSI FORTRAN RELEASE III NAME EX4B OG C UNIVERSITY OF WATERLOO JANUARY 1984 CHT SHS SSS c 0000000090000 0860606000 0000090000 000900000000000089000 Cee soneo61od 6046090066680 8060 8000806060096 0606040096060686900 90060680009 a U S MA I N L I N E P R O G R A M bab dd Co svo o do to dod 0000000032080 80800080600909006080600000808000060090680000800 C 9 O gt o gt oO BOBO Boo BOB Bo dd sl 0800800060680 0008090000000 0 C EC FILE REQUIREMENT C UNIT 3 FOR RSTRTA A Mia AAA AAA AAA AAA AAA AAA AAA AAA AAA AAA EEE AAA EA E RAEE S INTEGER I TERRA IPRNTE IPRNTS MAXINT MAXSA 1 MSGLVA NVARS REAL MCHEPS RATIOL RATIOS TIME REAL S 250 REAL e ERROR ONE TWO ZERO C Ctra cr sernnreesrr Sn nenrVrrarno9SSSSRnnnessnnoo C i A COMMON SPAUSR MSGLVA IERRA MAXSA NVARS COMMON SPKSYS IPRNTE IPRNTS MAXINT RATIOS RATIOL QQAQQQQ 1 i MCHEPS TIME O I RR RN INITIALIZE SPARSPAK A CALL SPRSPK MAXSA m 250 C Bw Se
27. 1 and the entries i in the right hand side vector b are all ones REAL 250 FOUR ONE INTEGER I IERRA MAXSA MSGLVA NEGNS COMMON SPAUSR MSGLVA IERRA MAXSA NEQNS G CALL SPRSPK MAXSA 250 INPUT THE MATRIX STRUCTURE THE DIAGONAL IS ALWAYS ASSUMED TO BE NONZERO AND SINCE THE MATRIX IS SYMMETRIC SPARSPAK A ONLY NEEDS TO a KNOW THAT THE SUBDIAGONAL ELEMENTS ARE NONZERO CALL IJBEGN DO 100 I 2 10 CALL INIJ I l 1 100 CONTINUE WW CALL IJEND S Qaaaaa w e a ms e e e ee e e e e we a e oe a e e e w oe e n em we a e e e mm w ep w a we am e w e e w e eo INPUT THE NUMERICAL VALUES FOR A AND B SINCE THE MATRIX IS SYMMETRIC ONLY THE LOWER TRIANGLE AND THE DIAGONAL ARE I NPUT 099a AAA 1 Declared either REAL or DOUBLE PRECISION depending on the version of SPARSPAK A that is available The exam ples in this manual assume a single precision version is being used November 1984 a 5 SPARSPAK A User s Guide erate AR A A a a a A FOUR 4 0E0 ONE 1 0E0 DO 200 I i 10 IF 1 GT 1 CALL INAIJ1 1 l 1 ONE CALL INAIJ1 I I FOUR CALL INBI I ONE 200 CONTINUE CG vao EE se c SOLVE THE SYSTEM g lt lt V V T T T T CALL SOLVEL E A AAA A E c PRINT THE SOLUTION FOUND IN THE FIRST TEN O LOCATIONS OF THE WORKING STORAGE ARRAY S C A AD WRITE
28. 66SPK 67SPK 68SPK 69SPK gt TOSPK gt TISPK 72SPK 78SPK 74SPK 75SPK 16SPK T1SPK 78SPK 79SPK 80SPK 81SPK 82SPK 83SPK 84SPK 85SPK 86SPK gt 87SPK 88 SPK 89SPK 90SPK 91SPK 9 2 SPK 98SPK 94 SPK 95SPK 96SPK 97SPK 98SPK 99SPK 100SPK 101SPK gt 108SPK 103SPK 104SPK 105SPK 106SPK 107SPK 108SPK 109SPK 110SPK 111SPK 112SPK 118SPK 114SPK 115SPK SPARSPAK A 1 User s Guide Output 464 444 UNIVERSITY OF WATERLOO eeeeeeeees SPARSE MATRIX PACKAGE SPARS PAK seeeeenense RELEASE 8 eanaeeeeee C JANUARY 1984 4 ANSI FORTRAN SINGLE PRECISION 9 LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES 6 OUTPUT UNIT FOR STATISTICS s 6 IJBEGN BEGIN STRUCTURE INPUT INIJ vo INPUT OF ADJACENCY PAIRS IJEND END OF STRUCTURE INPUT ORDRAS ONE WAY DISSECTION ORDERING INAIJS INPUT OF MATRIX COMPONENTS INBI INPUT OF RIGHT HAND SIDE SOLVES IMPLICIT BLOCK SOLVE ERESTS ERROR ESTIMATOR STATSA SYSTEM A STATISTICS SIZE OF STORAGE ARRAY MAXSA 250 NUMBER OF EQUATIONS 10 NUMBER OF OFF DIAGONAL NONZEROS 18 TIME FOR ORDERING 0 033 STORAGE FOR ORDERING 91 TIME FOR ALLOCATION 0 017 STORAGE FOR ALLOCATION 94 s STORAGE FOR SOLUTION 94 TIME FOR FACTORIZATION E 7 TIME FOR SOLUTION 0 017 OPERATIONS IN FACTORIZATION 18 OPERATIONS IN SOLUTION S TIME FOR ESTIMATING RELATIVE ER
29. 9SPK 10SPK 11SPK 12SPK 19SPK 14SPK 15SPK 16SPK 17SPK 18SPK 19SPK 20SPK 21SPK 22SPK 9SPK 24SPK 25SPK 26SPK 27SPK 28SPK 29SPK 30SPK 31SPK 32SPK 33SPK 34SPK 95SPK 96SPK 37SPK 38SPK 99SPK 40SPK 41SPK 42SPK 48SPK A4 SPK 4 5SPK 4 6SPK 4 71 SPK 48SPK 49SPK 50SPK 51SPK 52SPK 59SPK 54SPK 55SPK 56SPK 57SPK 58SPK 59SPK 60SPK 61SPK 62SPK 63SPK 64SPK 65SPK 66SPK 6 7SPK 68SPK 69SPK SPARSPAK A 09999 AAA 900 11 User s Guide COMPUTE THE ACTUAL RELATIVE ERROR IN THE COMPUTED SOLUTION SINCE THE TRUE SOLUTION IS KNOWN ERROR ZERO DO 300 I 1 00 ERROR AMAX1 ERROR ABS S I ONE CONTINUE WRITE IPRNTS 11 FORMAT 10X ERROR l 5HMAX IMUM ERROR a 1PE15 5 STOP END Output November 1984 ee se ee ee RSTRTA ORDRAS EMSGA SAVEA STATSA UNIVERSITY OF WATERLOO SPARSE MATRIX PACKAGE SPARSPAK RELEASE 8 C JANUARY 1984 ANSI FORTRAN SINGLE PRECISION LAST UPDATE JANUARY 1984 RARA AAA EXA ERA AAA AAA EEA EE RR ER EARAARRAR AA RARA OA PA A A A AO OUTPUT UNIT FOR ERROR MESSAGES 6 OUTPUT UNIT FOR STATISTICS 6 RESTART SYSTEM A NESTED DISSECTION ORDERING S YSTEM A ERROR ORDRXI X A B AND I 1 2 8 4 5 6 ERROR NUMBER INSUFF STORAGE FOR SOLVEI MAXSA MUST AT LEAST BE 127 2509 SAVE STORAGE VECTOR SYSTEM A STATISTICS 2500 200 398 0 167 1400 0 050
30. ATIONS IN FACTORIZATION OPERATIONS IN SOLUTION TIME FOR ESTIMATING RELATIVE ERROR OPERATIONS IN ESTIMATING REL ERROR STORAGE FOR ESTIMATING REL ERROR ESTIMATE OF RELATIVE ERROR TOTAL TIME REQUIRED MAXIMUM STORAGE REQUIRED MAXIMUM ERROR S 18 38 0 050 160 990 1 785e 07 0 100 90 1 19209e 07 User s Guide 34 SPARSPAK A l all IMMA User s Guide Example 2 This is the same as Example 1 except that the matrix A is unsymmetric The diagonal elements of A are all 4 the superdiagonal elements are all 1 and the subdiagonal elements are all 1 The right hand side vector b is chosen so that the entries of the solution vector z are all ones Program C SPARSPAK A ANSI FORTRAN RELEASE III NAME EX2 C C UNIVERSITY OF WATERLOO JANUARY 1984 CARER REE EEE kk k k EERE odb Lekt MA NL I NE P ROGRAM L EE CARRERA ARA AAA AAA AAA ARA AAA AAA AER AAA AAA RA ARA AAA AAA CEE ERE EE RRA EERE AAA RARA ERA ARA AAA AAA ARA RARA AAA AAA AAA AAA C INTEGER I IERRA IPRNTE IPRNTS MAXINT MAXSA 1 MSGLVA NVARS REAL MCHEPS RATIOL RATIOS TIME REAL S 250 REAL ERROR FOUR ONE RELERR ZERO HRA RRR EERE HS o COMMON SPAUSR MSGLVA IERRA MAXSA NVARS COMMON SPKSTS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME ZA V r EE C a re C INITIALIZE SPARSPAK A Co EC CALL SPRSPK MAXSA 250 C z sisse em ieS Se fe INPUT STRUCTURE C A
31. AXSA NUMBER OF EQUATIONS NUMBER OF OFF DIAGONAL NONZEROS TIME FOR ORDERING STORAGE FOR ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION STORAGE FOR SOLUTION TOTAL TIME REGUIRED MAXIMUM STORAGE REGUIRED November 1984 User s Guide 0 017 2396 2398 0 288 2396 9500 1194 0 217 2696 0 083 3599 8801 3599 59 SPARSPAK A Users Guide 10 Appendix implementation overview In this section we describe briefly the use of labelled common blocks in the internal implementation of SPARSPAK A and the various methods of communication between modules 10 1 User module communication As described in previous sections of this user guide the user supplies a one dimensional floating point array S from which all array storage is allocated In particular the interface allocates the first NEQNS storage locations in S for the solution vector of the linear system of eguations After all the interface modules for a particular method have been successfully executed the user can retrieve the solution from these NEQNS locations There is one labelled common block SPAUSR that the user must provide having four variables l COMMON FALSE MSGLVA IERRA MAXSA NEQNS The variable MAXSA is the declared size of the aie inenzi ijal floating point arrav S and it must be set by user at the beginning of his program For each module in the interface that allocates storage e g INIJ IJEND ORDRzi
32. E ACTUAL RELATIVE ERROR IN THE COMPUTED SOLUTION SINCE THE TRUE SOLUTION IS KNOWN em a n a e e e e e e a au e e e e w we e n m a a w m e e e we e e e e o e e e wm e e e e we w e w e a m a a o a we a nm ERROR ZERO DO 00 I 1 200 ERROR AMAX1 ERROR ABS S 1 ONE CONTINUE WRITE IPRNTS 11 ERROR ESA FORMAT 10X 85HMAXIMUM ERROR 1PE15 5 STOP END KARA AA ARA AAA ERA RARA UNIVERSITY OF WATERLOO SPARSE MATRIX PACKAGE SPARSPAK RELEASE 9 C JANUARY 1984 ANSI FORTRAN SINGLE PRECISION LAST UPDATE JANUARV 1984 ERA k ktkhk hkiktkkha DT OUTPUT UNIT FOR ERROR MESSAGES 6 OUTPUT UNIT FOR STATISTICS 6 RSTRTA RESTART SVSTEM A ORDRA5 NESTED DISSECTION ORDERING INAIJS INPUT OF MATRIX COMPONENTS INPUT OF RIGHT HAND SIDE SOLVES GENERAL SPARSE SOLVE SYSTEM A STATISTICS 2509 200 398 0 167 SIZE OF STORAGE ARRAY MAXSA NUMBER OF EQUATIONS NUMBER OF OFF DIAGONAL NONZEROS TIME FOR ORDERING 55 59SPK GOSPK 61SPK G2SPK 63SPK 64 SPK 65SPK 66SPK 67SPK 68 SPK 69SPK 70SPK 71SPK T2SPK 73SPK gt T SPK gt 75SPK TESPK 77SPK 78SPK 79SPK 80SPK 81SPK 82SPK 83 SPK 84SPK 8 5SPK 86SPK SPARSPAK A STORAGE FOR ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION STORAGE FOR SOLUTION TIME FOR FACTORIZATION TIME FOR SOLUTION OPERATIONS IN FACTORIZATION OPERATIONS IN SOLUTION TOTAL TIME REGUIRED MAXIMUM STORAGE REQ
33. IMATE OF THE RELATIVE ERROR IN THE COMPUTED SOLUTION w a e v e e e we m W m m e A sw e w e e e ep e CALL mp m e e e we G m e e e e w m w e e e e e e e e au we e e e nm e m w a e e e e e e e e e a e e au e a w e a s s e e e e e e e e e e s CALL STATSA SAVEA 8 STOP END eeeeeeee UNIVERSITY OF WATERLOO eee ee eee te SPARSE MATRIX PACKAGE z 7 SPARS PAK eek ik RELEASE 8 RRE C JANUARY 1984 Z it ANSI FORTRAN weeeeeeeete SINGLE PRECISION trr LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES OUTPUT UNIT FOR STATISTICS IJBEGN INIT I JEND ORDRA1 INAIJ1 SOLVE1 BEG IN STRUC TURE INPUT INPUT OF ADJACENCY PAIRS END OF 2 h CTURE INPUT RCM ORDER ING INPUT OF MATRIX COMPONENTS ENVELOPE SOLVE NO RIGHT HAND SIDE PROVIDED SOLUTION WILL BE ALL ZEROS EREST1 STATSA November 1984 ERROR ESTIMATOR SYSTEM A STATISTICS INAIJ1 I I FOUR S User s Guide 2 SISPK 52SPK 5 8SPK 54SPK 55SPK 56SPK 57SPK 58SPK 59SPK 60SPK 61SPK 6 2 SPK 63SPK 64SPK 65SPK 66SPK 67SPK 68SPK 69SPK TOSPK TISPK T2SPK 18SPK 74 SPK 75 SPK 76SPK SPARSPAK A User s Guide SIZE OF STORAGE ARRAY MAXSA 250 NUMBER OF EQUATIONS 10 NUMBER OF OFF DIAGONAL NONZEROS ES 18 TIME FOR ORDERING 0 STORAGE FOR ORDERING 60 TIME FOR ALLOCATION 0
34. INCLQ i is useful Its calling seguence is a CALL INCLQ NCLQ CLQ S November 1984 8 SPARSPAK A User s Guide where NCLQ is the size of the submatrix and Gier is an array containing the column or row indices of the submatrix Thus to inform the package that the submatrix corresponding to indices 1 3 5 and 6 is full we execute CLQ 1 1 CLQ 2 8 CLQ S CLQ 4 6 CALL INCLQ 4 CLA The type of structure input routine to use depends on how the user obtains the matrix structure Anyway the one or ones that best suit the application can be selected SPARSPAK A allows mized use of the routines in inputting a matrix structure The package automatically removes duplications so the user does not have to worry about inputting duplicated index pairs CALL INROW 4 3 IR CALL INIJ 1 3 S The code above would input the matrix structure k ix x into the package Note that the diagonal elements of the input matrix are assumed to be nonzero see notes below When all pairs have been input using one or a combination of chee input routines the 1 user is required to tell SPARSPAK A explicitlv that structure EE is iii by calling the routine IJEND The statement to use is November 1984 SPARSPAK A ko a E User s Guide CALL IJEND S and its purpose is to transform the data from the format used during the recording phase to the standard
35. INUE SR CALL INBI 1 ONE CALL INBI 10 ONE S m e e e o e e we a e e m e e e e we we e e e e e e e e oe e e e me an e we e e e e e e e e e e e e e e e ge e e e we e e e e e e e e we we e e e we a a s e e ee e e a we we we e e we me e e e CALL SOLVES S a e e e e e s s ERER w au w a wm m s e me S e S S e S S S e e e ann COMPUTE AN ESTIMATE OF THE RELATIVE ERROR IN THE COMPUTED SOLUTION a we e e e e e e e e e e e e e e e e e e am ge au e am e o a o a e e e e e e e e e e e CALL STATSA e m e e w w w e w u we w e e e m e e e S S e e e o e we e wm e s om s e ew s e s e e s em e e COMPUTE THE ACTUAL RELATIVE ERROR IN THE SOLUTION SINCE THE TRUE SOLUTION IS KNOWN ERROR ZERO DO 00 I 1 10 ERROR AMAX1 ERROR ABS S 1 ONE CONTINUE WRITE IPRNTS 11 ERROR FORMAT 10X S5HMAXIMUM ERROR a e a m w m m ewe a e e e e e e e e e e e e we DO 400 I 1 10 CALL INBI I FOUR CONTINUE CALL INBI 1 TWO CALL INBI 10 TWO w m w e ew w we e w e e wg w w e m w e w e w a CALL STATSA m m e a m m e e w w e w o e e W a m a m e mm we e e n ees a ERROR ZERO DO 500 I 1 10 ERROR AMAX1 ERROR ABS S 1 TWO CONTINUE WRITE IPRNTS 11 ERROR STOP END November 1984 1PE15 5 39 50SPK SISPK 52SPK SSSPK 54SPK 55SPK 56SPK 57SPK 58SPK 59SPK _ 60SPK 61SPK 62SPK 68SPK 64SPK 65SPK
36. MAXSA is used to make sure that en is at storage to carry out the particular phase 10 2 Module module SE There are several labelled common blocks used for communication among modules within the interface Two important ones are the control block SPACON and the storage map Block SPAMAP COMMON SPACON STAGE MXUSED MXREQD NVARS NEDGES METHOD and other method related control variables COMMON SPAMAP PERM INVP RHS and other method related data structure pointers i The control block has fifty integer variables and contains control information about the specific problem being solved There are fifty variables in the storage map block which keep the locations origins in S of the various arrays used in the particular storage scheme These storage schemes differ in complexity across the methods so the same storage map block must be used in the corresponding routines ORDRxi INAIJi INROWi INMATi SOLVEi and ERESTi An example is given below November 1984 _ A Ni G ia 60 SPARSPAK A d Se User s Guide RHS ane right hand side vector PERM gt permutation vector INVP ia inverse permutation vector XENV E index to E SE of L DIAG gt diagonal of the matrix factor L envelope of the matrix factor L Storage allocation for the symmetric envelope method ORDRA1 10 3 Save and restart implementation The SAVEA routine saves the control information in the control block
37. NTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME c CHAR KRKA RAKA c A A A ete C INITIALIZE SPARSPAK A November 1984 wa S 47 ISPK 2SPK SSPK 4 SPK 5SPK 6SPK 7SPK 8 SPK 9SPK 10SPK 11SPK 12SPK 18SPK 14SPK 15 SPK 16SPK 17SPK 18SPK 19SPK 20SPK 21SPK 22SPK 23SPK 24SPK 25SPK 2 6SPK 2 7SPK 2 8SPK 29SPK 30SPK QQQQ QAQ A QQQ Q 100 200 300 400 0000 AAAA AAA 500 11 a a a e s s e e S S S e S a s e CALL SPRSPK MAXSA 900 w a e v e e e e s e we we INPUT STRUCTURE CALL IJBEGN DO 100 Te 2 00 CALL INIJ l CONTINUE CALL IJEND S IF IERRA EQ 0 CALL STATSA STOP I 1 S Go TO 200 CONTINUE CALL ORDRAS S IF IERRA EQ 0 CALL CALL STOP GO TO 300 SAVEA 9 STATSA CONTINUE e an s e 4 e e e we 0E0 I 1 00 IF 1 GT 1 INBI L 1 ONE S 200 ONE S 0 i 2 0E0 4 CALL TWO S INAIJ5 I 1 1 ONE S CONTINUE CALL INBI CALL INBI e e e w e e ge am e e e e e e e we e e e e e e e e we e e e e e e e oe e e e e e e m e e e a a G e 6 s e v m am e we s gw we NM a we e s e e we e we ww e a CALL SOLVES S m e a G a e e e e u e e e we e e e m a e e e a e e e a e a e e 5 e e e e e e e e e e e e e e
38. ONE 1 0E0 TWO 2 0E0 FOUR 4 0E0 DO 400 I 1 200 f IF I GT 1 CALL INAIJ5 I 1 1 ONE S CALL INBI 1 TWO S CONTINUE CALL INBI 1 ONE CALL INBI 200 ONE em e e e e e e e e e e e e we e e a a e we e e 5 e e e e e e e e e e e a e we e e e e CALL SOLVES S CALL STATSA e e e e e e e e e e s a gw e e e e e e e e e e e e s e e e e e e e e a n a a e e a a e COMPUTE THE ACTUAL RELATIVE ERROR IN THE COMPUTED SOLUTION SINCE THE TRUE SOLUTION IS KNOWN e e a e e e e e m w e e e e e e e e e e e e e e s e e o e we e e e e e e e e a e e e a e w m u w m e e we a e e ERROR ZERO DO 500 I 1 200 November 1984 User s Guide 50 2 SPK 23SPK 24SPK 25 SPK 26SPK 27SPK 28SPK 29SPK 30SPK 31SPK 32SPK 33SPK 34SPK 35SPK 36SPK 3 7SPK 38SPK 39SPK 40SPK 41SPK 42SPK 49SPK 44SPK 45SPK 46SPK 4 7SPK 48SPK 49SPK SOSPK 51SPK 52SPK 53SPK 54SPK 55SPK 56SPK 5 7SPK 58SPK 59SPK 60SPK 61SPK 62 SPK GISPK 64 SPK 6 5SPK 6 6SPK 67SPK 68SPK 69SPK 70SPK 71SPK 72SPK 73SPK 74 SPK 75SPK 76SPK 77SPK 78SPK 79SPK 80SPK 81SPK 8 2SPK 8 9SPK 84 SPK 85 SPK 86SPK 8 7SPK SPARSPAK A C ERROR AMAX1 ERROR ABS S 1 ONE 500 CONTINUE WRITE IPRNTS 11 ERROR 11 FORMAT 10X S5HMAXIMUM ERROR STOP END User s Guide 88SPK 89SPK 90SPK 1PE15 5 91SPK 92SPK 93SPK
39. ROR 0 067 OPERATIONS IN ESTIMATING REL ERROR 157 STORAGE FOR ESTIMATING REL ERROR 144 ESTIMATE OF RELATIVE ERROR 1 785e 07 TOTAL TIME REQUIRED MN 0 133 MAXIMUM STORAGE REQUIRED 144 MAX IMUM ERROR 1 19209e 07 INBI INPUT OF RIGHT HAND SIDE SOLVES IMPLICIT BLOCK SOLVE FACTORIZATION ALREADY DONE STATSA SYSTEM A STATISTICS November 1984 S A0 SPARSPAK A SIZE OF STORAGE ARRAV MAXSA NUMBER OF EQUATIONS NUMBER OF OFF DIAGONAL NONZEROS TIME FOR ORDERING STORAGE FOR ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION STORAGE FOR SOLUTION TIME FOR FACTORIZATION TIME FOR SOLUTION OPERATIONS IN FACTORIZATION OPERATIONS IN SOLUTION TOTAL TIME REQUIRED MAXIMUM STORAGE REQUIRED MAXIMUM ERROR User s Guide 18 88 0 050 94 2 98419e 07 November 1984 41 SPARSPAK A A User s Guide Example 4 This example illustrates the use of the save sala feature of SPARSPAK A After the factorization is computed SAVEA is executed which writes the current state of the computation on FORTRAN logical unit 3 In a second program the module RSTRTA is execute to read the information from unit 3 and the computation resumes at the point at which SAVEA was invoked Program 1 C SPARSPAK A ANSI FORTRAN RELEASE III NAME EX4A C C UNIVERSITY OF WATERLOO JANUARY 1984 NoidAd d d d d ddd d d ddddd ddddddd addddddddddddddddddd ddddaka etesro
40. RRA IPRNTE IPRNTS MAXINT MAXSA 1 MSGLVA NVARS REAL MCHEPS RATIOL RATIOS TIME REAL S 8500 REAL ERROR FOUR ONE TWO ZERO c ISI TITIS ISIT TISI III ISIT ITYYYITYTYYTTITYITITTITTITITTITYITTITE C an 7 COMMON SPAUSR MSGLVA IERRA MAXSA NVARS COMMON SPKSYS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME ADDED DADA DADA DAA ADA AAA AAA AAA AAA AAA AAA AAA AA Zn e e e wens e 4943 e Z toua 4 mu e tr a N ta cl gt a Lo a CALL SPRSPK MAX SA 2500 we s e a e MA we G e e m au a w w 46 e e 0 e e e 6 e m a e e e 4 e s a mp e e e e e 6 we V DB e 6 e e s 06 e e e e e 0999 QQQAQ CALL ORDRAS S IF IERRA EQ 0 GO TO 100 CALL SAVEA 3 S CALL STATSA STOP 100 CONTINUE QQQ a z u G H Q KN E MES 2 tal a e e e e s e we e am e e e s 0 de 2 0E0 4 0E0 DO 200 I 1 200 IF 1 GT 1 CALL INAIJ5 I l 1 ONE CALL INAIJS I I FOUR S CALL INBI I TWO S 200 CONTINUE CALL INBI 1 ONE CALL INBI 200 ONE S m e e e e e S a au e e e e e e e we e e ss CALL SOLVES S YW QQQQ 3 Bi y ra m te fer Q Y N p 3 a ta Q te G bj ou lt November 1984 User s Guide 52 4 SPK SSPK 6SPK 7SPK 8SPK
41. SPARSPAK Waterloo Sparse Matrix Package User s Guide for SPARSPAK A Eleanor Chu Alan George Joseph Liu Esmond Ng Department of Computer Science University of Waterloo Waterloo Ontario CANADA Research Report CS 84 36 November 1984 Table of Contents 1 Introduction and basic structure of SPARSPAK A a S 2 Modules of SPARSPAK A and how to use them ee 5 2 1 User mainline program and an example sssesennzzennznnnznszzsnenzonnezsanzaz A 2 2 Modules for input of the matrix structure py D 7 a Input of a nonzero location conte a a 0 b Input of the structure of a row or part of a TOW ooooonconncnicnonnnicnnonconnconacnss ad c Input of a submatrix structure skoku i i tn kaa ot e ene 8 d Input of a full submatrix structure ENEE 8 2 3 Modules for ordering and storage allocation AR ke TAT AT Pat 10 2 4 Modules for inputting numerical values ssssssserenzennenennenezzzonosonenznzarasnan 12 a Input of a single nonzero component REENEN send 12 b Input of a row of nonzeroes Ee E e Input of a sub matrix iii kie osa dk BE 13 2 5 Modules for numerical solution SA ee 14 2 6 Modules for estimating the relative error in the computed solution 14 3 Some guidelines on selecting a method A Ta 16 4 Save and restart facilities 0rrrrrrrreeererereeaereree a Eeer 18
42. SPK END 73SPK Output 1 14 UNIVERSITY OF WATERLOO eee ee ee et SPARSE MATRIX PACKAGE Pee a ae kkt SPARSPAK EZE E RELEASE 8 EXA AAA C JANUARY 1984 eee ee eeeee ANSI FORTRAN ARALAR AR SINGLE PRECISION hkhk kAS LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES 6 OUTPUT UNIT FOR STATISTICS 6 RSTRTA RESTART SVSTEM A INBI INPUT OF RIGHT HAND SIDE SOLVE1 ENVELOPE SOLVE FACTORIZATION ALREADY DONE STATSA SYSTEM A STATISTICS November 1984 ES 45 SPARSPAK A SIZE OF STORAGE ARRAY MAXSA NUMBER OF EQUATIONS NUMBER OF OFF DIAGONAL NONZEROS TIME FOR ORDERING STORAGE FOR ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION STORAGE FOR SOLUTION TIME FOR FACTORIZATION TIME FOR SOLUTION OPERATIONS IN FACTORIZATION OPERATIONS IN SOLUTION TOTAL TIME REQUIRED MAXIMUM STORAGE REQUIRED MAXIMUM ERROR 0 017 0 017 9 050 1 19209e 07 User s Guide November 1984 46 SPARSPAK A User s Guide Example 5 This example consists of four runs of essentiallv the same program illustrating how the SAVEA and RSTRTA modules can be used to avoid repeating successfully completed computations when the execution cannot proceed further because of lack of working storage In the first run MAXSA was too small to accommodate the structure and a message was printed indicating that MAXSA must be at least 999 in order to input the structure A second run with MAXSA 2999
43. SPK MAXSA 250 C GC aaa e mm e mm ss s C INPUT STRUCTURE C ismi EN c 08 CALL IJBEGN I 2 10 DO 100 November 1984 S e 32 1SPK 2SPK SSPK 4SPK 5SPK 6SPR gt 7SPK 8SPK 9SPK 10SPK 11SPK 12SPK 13SPK 14SPK 15SPK 16SPK 17SPK 18SPK 19SPK 20SPK 21SPK 22SPK 23SPK 24SPK 25SPK 26SPK 27SPK 28SPK 29SPK 30SPK 91SPK 32SPK 33SPK 34SPK User s Guide 1PE15 5 SPARSPAK A CALL INIJ 1 1 1 S 100 CONTINUE CALL IJEND S C Ej nn n C DETERMINE SYMMETRIC ORDERING C jj ese eee eee wren eee eee ee mm e een sae CALL ORDRA1 S Cc Um mm mm me mm mm ooo C INPUT NUMERICAL VALUES Uewe eee e e sm ams a e m ZERO 0 0E0 ONE 1 0E0 TWO 2 0E0 FOUR 4 0E0 DO 200 Im 1 10 IF I GT 1 CALL INAIJ1 1 1 1 ONE S CALL INAIJ1 I FOUR CALL INBI I TWO S 200 CONTINUE Si CALL INBI 1 ONE S CALL INBI 10 ONE C l G i A p a ta ee a es a ee ed e C PERFORM NUMERICAL FACTORIZATION AND SOLUTION A SN E CALL SOLVEL c G n EE ee ee C COMPUTE AN ESTIMATE OF THE RELATIVE ERROR C gt IN THE COMPUTED SOLUTION C eee e er tae een E a a a ea a aj oe CALL EREST1 RELERR C CO vb ya Cc OBTAIN STATISTICS o sms rr o e e a CALL STATSA C VVT C COMPUTE THE ACTUAL RELATIVE ERROR IN THE COMPUTED SOLUTION ei SINCE THE TRUE SOLUTION IS KNOWN C TBT da r le i ee ERROR m ZERO DO 300 Im 1 10 ERROR AMAX
44. STATISTICS 61SPK fo ANA nosne 62SPK CALL RSTRTA 3 63SPK CALL ORDRAS S 64SPK CALL STATSA 65SPK C 66SPK STOP 67SPK END 68SPK Output eeeseseees UNIVERSITY OF WATERLOO eeseeeseses SPARSE MATRIX PACKAGE eesesesoss SPARS PAK ttst st ssSSs RELEASE 3 torrsnrr C JANUARY 1984 sestassS82 ANSI FORTRAN sesesssess SINGLE PRECISION d ii ake ES LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES OUTPUT UNIT FOR STATISTICS EN JJBEGN BEGIN STRUCTURE INPUT INIJ INPUT OF ADJACENCY PAIRS IJEND END OF STRUCTURE INPUT SAVEA SAVE STORAGE VECTOR ORDRAL RCM ORDERING STATSA SYSTEM A STATISTICS SIZE OF STORAGE ARRAY MAXSA 9500 NUMBER OF EQUATIONS E 800 NUMBER OF OFF DIAGONAL NONZEROS 1194 TIME FOR ORDERING 0 217 STORAGE FOR ORDERING 2396 November 1984 t 4 e w 58 SPARSPAK A TIME FOR ALLOGATION STORAGE FOR ALLOCATION STORAGE FOR SOLUTION TOTAL TIME REQUIRED MAXIMUM STORAGE REQUIRED RSTRTA RESTART SYSTEM A ORDRAS ONE WAY DISSECTION ORDERING STATSA SYSTEM A STATISTICS SIZE OF STORAGE ARRAY MAXSA NUMBER OF EQUATIONS NUMBER OF OFF DIAGONAL NONZEROS TIME FOR ORDERING STORAGE FOR ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION STORAGE FOR SOLUTION TOTAL TIME REQUIRED MAXIMUM STORAGE REQUIRED RSTRTA RESTART SYSTEM A ORDRAS NESTED DISSECTION ORDERING STATSA SYSTEM A STATISTICS SIZE OF STORAGE ARRAY M
45. TORAGE FOR ORDERING 60 TIME FOR ALLOCATION 0 STORAGE FOR ALLOCATION 760 STORAGE FOR SOLUTION 79 TIME FOR FACTORIZATION 0 017 TIME FOR SOLUTION 0 OPERATIONS IN FACTORIZATION 18 OPERATIONS IN SOLUTION W 28 TIME FOR ESTIMATING RELATIVE ERROR 0 088 OPERATIONS IN ESTIMATING REL ERROR 160 STORAGE FOR ESTIMATING REL ERROR 99 ESTIMATE OF RELATIVE ERROR 1 197e 07 TOTAL TIME REGUIRED 0 067 MAXIMUM STORAGE REGUIRED 99 MAXIMUM ERROR 1 19209e 07 November 1984 37 SPARSPAK A EI User s Guide Example 3 This is similar to Example 1 except that method 3 is used with the A ordering option and two problems differing only in their right hand sides are solved After solving the problem whose solution vector contains all ones a new right hand side is input which corresponds to a different problem whose solution vector contains all twos When the module SOLVES is called a second time the interface detects that the factorization has already been done and only the triangular solution i is performed Program C SPARSPAK A ANSI FORTRAN RELEASE III NAME EX3 C UNIVERSITY OF WATERLOO JANUARY 1984 ABAD AAA AAA AAA AAA AAA AAA AAA AAA ERE RE REA AREA AAA AAA ATT NY Ctrrrrnsrranrannnnnananannanosssnarararar onoossso Corr rro n MAINLINE PROGRAM PEE EEEEET ABBA AAA AAA AAA AAA AAA AAA AAA AAA AAA AAA AA AAA AE RAEE RAE AE AAA AAA ATEN O A EE LA EEN C INTEGER I IERRA IPRNTE IPRNTS MAXINT MAXSA LO MSGLVA NVARS
46. Thus SPRSPK is an installation dependent subroutine November 1984 ikali 63 SPARSPAKA User s Guide 11 References 1 2 Bi 4 5 6 7 E C H Chu and J A George An algorithm to estimate the error in Gaussian elimination without pivoting Research report CS 84 21 Department of Computer Science Universitv of Waterloo 1984 J A George An automatic one way dissection algorithm for irregular finite element problems SIAM J Numer Anal 17 1980 pp 740 751 J A George and J W H Liu Algorithms for matrix partitioning and the numerical solution of finite element systems SIAM J Numer Anal 15 1978 pp 297 327 J A George and J W H Liu An automatic nested dissection algorithm for irregular finite element problems SIAM J Numer Anal 15 1978 pp 1053 1069 J A George and J W H Liu The design of a user interface for a sparse matrix package ACM Trans on Math Software 5 1979 pp 134 162 J W H Liu On multiple elimination i in the minimum degree algorithm Technical Report No 83 03 Department of Computer Science York University Downsview Ontario 1983 JW H Liu and A H Sherman Comparative analysis of the CuthillMcKee and Reverse l Cuthill McKee ordering algorithms for sparse matrices SIAM J Numer Anal 13 1976 pp 198 213 November 1984 _ EN o 64
47. UIRED MAXIMUM ERROR 1400 0 950 2324 2509 0 0 2100 050 959 1168 0 367 2509 2 38419e 07 User s Guide November 1984 56 SPARSPAK A User s Guide Example 6 This is a program to illustrate how one might use SPARSPAK A to choose a method The matrix is 300X300 it has nonzeroes on the diagonal the first column and the last row The structure of the matrix is input using IJBEGN INIJ and IJEND and then saved on FORTRAN unit 3 The modules ORDRA1 ORDRAS and ORDRAS are then executed each one followed by a call to STATSA to obtain the storage information Note that RSTRTA is called after execution of ORDRA1 and ORDRAS to restore the package to the state that existed immediately after the structure inputting routines were executed Note also that SAVEA could have been used after each ordering module with different output unit numbers After one of the methods was chosen RSTRTA with the appropriate unit number could be used to initiate the computation avoiding re executing the ordering module corresponding to the method chosen SE Program C SPARSPAK A ANSI FORTRAN RELEASE III NAME EX6 C C UNIVERSITY OF WATERLOO JANUARY 1984 CARA AAA dk kkt kk kk kk kk bk kk kk dk k kd kd d dk dk dd bkd d AAA CRF tat k tk kd b bb kk k kkt dd db kd b k kk bb bod k db kkt AAA ARA dk d ARA AAA AAA Cert teeters MA I N L I N E P R O G R A M tb dd db bb bt CRAKEAR bo
48. al error in the Kaes approximate solution is determined since the true solution is known Note that the size of the working storage provided was 250 while the maximum amount used by any of the modules was 90 which was the storage requirement for the ERESTI module Thus if the user was going to solve this problem again he could adjust his storage down to 90 Program C SPARSPAK A ANSI FORTRAN RELEASE III NAME EX1 C C UNIVERSITY OF WATERLOO JANUARY 1984 C SRAORAORAOSANAONORARDSDADORANDDAONANAORAODOUDBANAN ANANA AHADHD HNAH AD ke CIAHAORAEARAADROCPDAONFNANAADAAODDOANREAEDRODNRANOODARODRBARAPDBSADBODBAANARANN Crr ssss MAINLINE P R OG RAM t ttstateaKs C OIPARONISRANDAODAODDADANAORADABDOCDDERODODDONONAEOHASBARODAOBHABDOBAHAN l d d d d d d d d dd d d dd ddd ddddd da A d d d C p S IERRA INTEGER I l IPRNTE IPRNTS MAXINT MAXSA 1 MSGLVA NVARS REAL MCHEPS RATIOL RATIOS TIME REAL 5 250 REAL ERROR FOUR ONE RELERR TWO 1 ZERO Ceeee ee tee ee eee eee e e e eee eee e tee ees ee Cc COMMON SPAUSR MSGLVA IERRA MAXSA NVARS COMMON SPKSVS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 c MCHEPS TIME a l ikk d d dd d ddd d d dddddd dddd d ddddddddd d d ddd dadad Z C i a s C INITIALIZE SPARSPAK A MM E ruzi scuscocssen CALL SPR
49. an be selected H e Department of Computer Science University of Waterloo Waterloo Ontario Canada ff Department of Computer Science York University Downsview Ontario Canada November 1984 SPARSPAK A Man get s Guide IMPORTANT NOTE The error estimate provided by the subroutine ERESTi see Section 2 6 is based on estimates of the condition number of the coefficient matrix and the error in the triangular factorization Experience has shown that the error estimate is very good in most cases lt may occasionally overestimate the actual relative error for some ill conditioned problems or badly scaled problems Also ERESTi may not be able to return an error estimate when the estimate of the condition number of the coefficient matrix or the estimate of the error in the triangular factors is large even though the solution may be computed accurately We would appreciate receiving any comments and feedback the user may have in using the error estimators for practical problems Such comments and feedback are important and useful for improving the error estimators Please send comments and feedback to Dr Alan George Department of Computer Science University of Waterloo Waterloo Ontario CANADA N2L 3G1 Telephone 519 885 1211 ext 3473 Dr Alan George November 1984 dl l o 2 SPARSPAK A 1 OO Users Guide 1 Introduction and basic structure of SPARSPAK A SPARSPAK A offers a collection of metho
50. angle to be input Attempt to input an M e element of matrix A where i gt n gon 1 lt 1 or jal NIN to input a numerical value for the 2 7 th element of matrix A into the data structure but the data structure has no space for it Probable cause of error the user has not called INIJ INROW INIJIJ or INCLQ with all the pairs j for which the i j th elements of A are nonzero dda sisien thinks A is sparser than it really is INBI INBIBI INRHS Incorrect execution seguence Probable cause of error routine ORDRzx1i was not executed successfully Index or subseript out of range Probable cause of error attempt to input a numerical value for the 1 th element of b where gt t n or i lt 1 SOLVEi Incorrect execution sequence Probable cause of error the numerical input routines were not executed successfully Incompatible ordering and solution routines Probable cause of error Routines ORDRxi and SOLVEJj were called where 1543 Response execute SAVEA and restart the computation using SOLVE where 2 is the value of METHOD specified in the error message 27 SPARSPAK A IERRA 153 154 User s Guide s SOLVEx Zero pivot or negative sguare root was detected i in the symmetric factorization routine Possible cause of error the matrix may require pivoting in order to preserve numerical stability In this case the use of SPARSPAK A to solve the problem is inappropr
51. ds for solving sparse systems of linear equations Ar b where A is an nXn nonsingular matrix and z and b are vectors of length n We assume the user is aware of the basic issues involved in solving sparse matrix equations and the basic facts about solving systems of linear equations using Gaussian elimination For a discussion on the initial design of this package see 5 For all the methods provided in SPARSPAK A the user and the package interact to solve the matrix problem through the following basic steps Step 1 The user supplies the nonzero structure of A to the package using a set of subroutines described i in Section 2 2 Step 2 The package reorders the gid probl finds a permutation P and allocates storage for the triangular factorization of PAP LU as described i in Section 2 3 Step 3 The user applies the numerical values for the matrix A to the package as described in Section 2 4 Step 4 The package computes the triangular factors L and U of f PAP as described i in Section 2 5 Step 5 The user supplies numerical values for b as described in Section 2 4 This step ka come before Step 4 and may be intermixed with Step 3 Step 6 The package computes the solution x using L U P and b as described in Section 2 5 Step7 The user optionally calls a subroutine which supplies an estimate of the relative error in z The subroutine is described in Section 2 6 The different methods provided in
52. e IJEND has already been called to indicate the end of structure input Insufficient storage was provided in the working storage array The 2 j pairs input to INIJ INROW INIJIJ and INCLQ will be counted and discarded Duplicates which are detected will not be counted but some duplicates mav be missed Input index or subscript is negative or zero IJEND Incorrect execution sequence Probable cause of error routine IJEND was called before new matrix structure has been input Call oy to start a new problem Insufficient storage to transform matrix structure The minimum value of MAXSA required is printed in the error message Number of variables NEQNS is zero 25 SPARSPAK A 7 3 3 Ordering and storage allocation routines TERRA 121 122 123 124 125 126 127 7 3 4 Input of the numerical values IERRA 131 132 November 1984 User s Guide ORDRri Incorrect execution sequence Probable cause of error routine IJEND was not executed successfully Incorrect execution sequence Probable cause of error routine ORDRzi was called after having already been executed successfully iet d Incompatible ordering method Probable cause of error part of the ordering routine ORDRzi was executed and then SAVEA Was executed because of insufficient storage The execution was then restarted using RSTRTA but ORDRrj was called with ij Insufficient storage in working sto
53. e the user can use the routine INRHS which transmits the whole vector to SPARSPAK A It has the form CALL INRHS RHS S where RHS is the vector containing the numerical values November 1984 ZEN Bo SPARSPAK A l o d F User s Guide In all three routines the numbers provided are added to those currentiv held b the package and the use of the routines can be intermixed See example in a above The storage used for the right hand side by SPARSPAK A i is initialized to zero the first time anv of them is executed Important Notes a When the matrix A is symmetric so that method 1 with odd is used SPARSPAK A requires that the elements of the lower triangle be provided Thus for example the to fonn statement will cause an error CALL INAIJS 3 5 1 8 S b The examples which we have steh assume that a single pren at version of SPARSPAK A is being used If the version is in double precision the numerical values and numerical variables should be declared as double precision For example CALL INAIJS 5 8 1 8DO S 2 5 Modules for numerical solution The numerical computation of the solution vector is initiated by the FORTRAN statement CALL SOLVEi S where S is the working storage array for SPARSPAK A Again the last mee 2 is used to distinguish between solvers for different storage methods Internally the routine SOLVE consists of both the factorization and forward backward solution steps If
54. e action which makes the error messages irrelevant Thus all printing by SPARSPAK A can be prevented by setting MSGLVA to zero To summarize we have November 1984 o 22 SPARSPAK A w User s Guide MSGLVA amount of output a 0 no information is provided 1 only warnings and errors are printed 2 warnings errors and summary are printed 3 warnings errors summary and some statistics are printed 4 detailed information for debugging purposes Warning It should be noted that by setting MSGLVA to four a high volume of output may be generated since the input data would also be echoed 7 2 Statistics gathering STATSA SPARSPAK A gathers a number of statistics which the user will find useful if he is comparing various methods or is going to solve numerous similar problems and wants to adjust his working storage to the minimum necessary The package has a common block called SPADTA containing variables whose values can be printed by executing the following statement CALL STATSA The information printed includes the number of equations the number of off diagonal nonzeroes in the matrix the size of the working Reeg the time used to find the ordering the time used for data structure set up the time used for the factorization step the time used for the triangular solution step the time used for the relative error estimation step number of operations required bv the factorization step numb
55. ee Se ee A ee ee C RESTORE STATE OF COMPUTATION ei See A e e e dere Sie N eent CALL RSTRTA 8 November 1984 44 1SPK 2 SPK 9SPK 4SPK 5SPK 6SPK 7SPK 8 SPK 9SPK 10SPK 11SPK 12SPK 13SPK 14SPK 15SPK 16SPK 1 7SPK 18SPK 1 9SPK 2 0SPK 21 SPK 2 2 SPK 23SPK 24 SPK 25SPK 26SPK 2 7SPK 23 SPK 29 SPK 90SPK 31SPK 32SPK 9 9SPK 94 SPK 9 5SPK 96SPK 7SPK SPARSPAK A User s Guide C l 98SPK CO mulas bola as id 39SPK C INPUT RIGHT HAND SIDE VECTOR 40SPK e SE l kou a PPP 41SPK ZERO 0 0E0 4 2SPK ONE 1 0E0 yo Un 49SPK TWO P OEO 44SPK DO 100 I 1 10 45SPK CALL INBI I TWO i 46SPK 100 CONTINUE i 47SPK CALL INBI 1 ONE S la 48SPK CALL INBI 10 ONE S 49SPK C CS 50SPK To O V a 518PK C PERFORM NUMERICAL SOLUTION v 59SPK C 2 DEER 59SPK GALL SOLVE1 S 54SPK C 55SPK Corr 56SPK C OBTAIN STATISTICS 57SPK V D 58 SPK CALL STATSA pout 59SPK C j 60SPK f A ee 61SPK C _ COMPUTE THE ACTUAL RELATIVE ERROR IN THE COMPUTED SOLUTION 62SPK C SINCE THE TRUE SOLUTION IS KNOWN 6SSPK Go RE EA ES EE 64SPK ERROR m ZERO 65SPK DO 200 I 1 10 66SPK ERROR AMAX 1 ERROR ABS S 1 ONE 6 7 SPK 200 CONTINUE 2 a 68 SPK WRITE IPRNTS 11 ERROR ZO a 69SPK 11 FORMAT 10X S5HMAXIMUM ERROR 1PE15 5 TOSPK C NE 71SPK STOP 72
56. ember 1984 y MN pa 6 SPARSPAK A User s Guide _ If these common block names cause conflicts in your program or at your computer installation it is possible to have the package distributed with these common blocks having speci fically requested labels These names should be specified when the package is acquired 2 2 Modules for input of the matrix structure SPARSPAK A has to know the matrix structure before it can determine an appropriate ordering for the system We now describe the group of routines which provide a variety of ways through which the user can inform the package where the nonzero entries are that is those indices j for which the 4 j th element of A is nonzero Before any of these input routines is called the user must execute an initialization routine called IJBEGN which tells the package that the structure of a new matrix problem is about to be input CALL IJBEGN a Input of a nonzero location To tell the package that the 7 7 th element of A is nonzero the user simply executes the statement i CALL INIJ I J S where J and J are respectively the row and column indices of the nonzero and S is the working storage array declared by the user for use by the package In this example I 4 J CALL INIJ 1 J S the package will record a logical nonzero in positions 4 3 and 3 4 of the matrix b Input of the structure of a row or part of a row When the struct
57. emmnnnnnnnnenzznnennzzonznen Bel w 60 10 3 Save and restart implementatioti EE 61 10 4 Method Checking iii eae iii a 62 10 5 Stage sequence checking En ee a Ee ta 62 10 6 Storage allocation of integer and floating point arrays o 63 10 7 Statistics gathering o re MA TI References E SE 64 November 1984 o 2 SPARSPAK Waterloo Sparse Matrix Package User s Guide for SPARSPAK A A collection of e for solving sparse systems of linear equations Eleanor Chu flasi George Jayi Liutt Esmond Ng Research Repor CS 84 36 O November 1984 ABSTRACT This document describes the structure and use of SPARSPAK A a sparse linear equations package which is designed to efficiently solve large sparse systems of linear eguations Computer programs for solving sparse systems of linear eguations typically involve fairly complicated data structures and storage management In many cases the user of such programs simply wants to solve his problem and should not have to understand how the storage management is done or how the matrix components are actually stored One of the attractive features of this package is that it effectively insulates the user from these considerations while still allowing the package to be used in a variety of ways Another important feature of the package is the provision of a variety of methods for solving sparse systems along with convenient means by which the best method for a given problem c
58. er of operations required by the triangular solution step number of operations required bv the relative error estimation step the storage used bv the ordering subroutine k the storage used by the data structure set up subroutine the storage used by the SOLVE module the storage used by the ERESTi module an estimate of the reciprocal of the condition number of the input matrix and November 1984 7 23 SPARSPAK A f E User s Guide an estimate of the relative error in the triangular factorization an estimate of the relative error in the computed solution Since the module STATSA can be called at any time some of the above information may not be available and will not be printed The word operations here means multiplicative operations multiplications and divisions Since most of the arithmetic performed in sparse matrix computation occurs in multiply add pairs the number of Operations as defined here i is a useful measure of the amount of arithmetic performed The reader is referred to the examples in Section 9 for more discussion about the output from STATSA 7 3 Error messages IERRA When a fatal error is detected so that the computation cannot proceed a positive code is assigned to IERRA The user can simply check the value of IERRA to see if the execution of module has been successful This error flag can be used in conjunction with the save restart feature described in Section 4 to retain the re
59. erical Values Call SOLVE for b Input Numerical Values Call SOLVEi for 6 Call SOLVE Flowchart 1 KON Flowchart 2 November 1984 AM 921 SPARSPAKA l User s Guide 7 Output from SPARSPAK A _ As noted earlier in Section 2 the user supplies a one dimensional floating point array S from which all array storage is allocated In particular the interface allocates the first NEQNS storage locations in for the solution vector of the linear system NEQNS is the last variable in the common block SPAUSR After all the interface modules for a particular method have been successfully executed the user can retrieve the solution from these NEQNS locations In addition to the solution vector SPARSPAK A may provide other information about the computation depending upon the value of MSGLVA whether or not errors occur and whether or not the module STATSA 1 is called This section discusses these features of SPARSPAK A Note SPARSPAK A writes output to two FORTRAN logical output units whose numbers are given by IPRNTS and IPRNTE The values for these variables are set in the module SPRSPK when the package is installed Standard output requested by the user is printed on unit IPRNTS while any error messages raised by SPARSPAK A are printed on unit IPRNTE In an interactive environment IPRNTE is usually the user s terminal while IPRNTS is some other output device on which the output of the
60. iate See restrictions in Section 1 Zero pivot was detected in the unsymmetric factorization routine Possible cause of error the matrix may require pivoting in order to preserve numerical stability In this case the use of SPARSPAK A to solve the propiem is inappropriate See restrictions in Section 1 7 3 8 Relative error estimation _IERRA 161 162 163 164 November 1984 Incorrect execution sequence Probable cause of error routine SOLVE was not executed successfully Incompatible condition number estimation routine Probable cause of error Routines ERESTi and SOLVE were called where 77 Response execute SAVEA and restart the computation using ERESTi where i is the value of METHOD specified in the error message Insufficient storage in working storage array to compute an estimate of the relative error in the computed solution Response execute SAVEA and restart the computation using ERESTi with MAXSA at least as large as that indicated in the error message The estimate of the relative error has a value of 1 0 which means that the computed solution may not have any correct significant digits 28 SPARSPAK A M 7 User s Guide 8 Summary listing of interface routines Initialization of SPARSPAK A SPRSPK Structure input sid input l l JJBEGN o IND 1 7 5 INROW I NIR IR S INIJIJ NIJ IL JJ S INCLQ NCLQ CLO S IJEND S ii Se Ordering see next table
61. k kk bb AAA COARA RA hu C C FILE REQUIREMENT C UNIT 3 FOR SAVEA AND RSTRTA Ctra eee ee ee ee ee Te Se dk k k dk RARA c INTEGER I REAL IERRA MAXSA MSGLVA NVARS S 7500 C e l Led bk kkt EE EE C COMMON SPAUSR MSGLVA IERRA MAXSA NVARS MO ra pi C DI DEER C INITIALIZE SPARSPAK A C ere A A eae a ee CALL SPRSPK MAXSA 7500 C e ZE C INPUT STRUCTURE C MAA Bl CALL IJBEGN DO 100 Im 1 300 CALL INIJ I 1 S CALL INIJ 300 I 100 CONTINUE CALL IJEND S C Co ab A NR C SAVE STRUCTURE INFORMATION VE O OA A ee ee Ss November 1984 E ME SH 57 20SPK 21SPK 22SPK 238SPK 24SPK 25SPK 26SPK 275PK 28SPK 29SPK SOSPK 1SPK 92SPK 8 3SPK 34SPK 35SPK 96SPK 3 7SPK 8 SPK 9SPK 40SPK 41SPK SPARSPAK A User s Guide CALL SAVEA 8 42 SPK C 43SPK Le EE EE rr a rr rr a 44SPK C DETERMINE REVERSE CUTHILL MCKEE ORDER ING AND 45SPK C OBTAIN STATISTICS 46SPK c a a a a a a a a e EE rr n a a rr a i a 47SPK CALL ORDRA1 S 48SPK CALL STATSA 49SPK C 50SPK fo A ni i iii rr ei 51SPK C RESTORE STRUCTURE INFORMATION DETERMINE ONE WAY 52SPK C DISSECTION ORDERING AND OBTAIN STATISTICS 59SPK C rr rr rr rr anan nna nannan 54SPK CALL RSTRTA 9 S 55SPK CALL ORDRAS 56SPK CALL STATSA 57SPK C 58SPK Oca rr rr eee ee 59SPK O RESTORE STRUCTURE INFORMATION DETERMINE NESTED 60SPK C DISSECTION ORDERING AND OBTAIN
62. lems say n lt 200 they are probably the best overall methods to use The ir of these methods is to reduce storage requirements but the factorization time will usually be substantially higher than any of the other methods Their storage requirements will usually be substantially less than methods 1 2 unless n is very large The same remark is true about the relative solution times Thus these methods are often useful when storage is restricted and or when many problems which differ only in the right hand side must be solved see Section 6 There are two ordering options provided ORDRAS and ORDRB8 and similarly for the unsymmetric case The A option is specifically tailored for finite element problems typical of those arising in structural analysis and the numerical solution of partial differential equations 2 The B option is effective for less specific problems and uses a refined quotient tree ordermg described mo B3 These methods attempt to find orderings which minimize fill in and they exploit all zeroes Their ordering times are almost always greater than those above but for moderate to large problems the reduced factorization times usually are more than PER 16 SPARSPAK A MN User s Guide 4 Save and restart facilities SPARSPAK A provides two subroutines called SAVEA and RSTRTA which allow the user to stop the calculation at some point save the results on an external sequential file
63. ng sequence is similar to that of INROW described on Section 2 2 CALL INROW 1 NIR IR VALUES S Here the additional parameter VALUES is a floating point array containing the numerical values of the row Again the numerical values are added to the current values in storage c Input of a submatrix The routine that allows the input of a submatrix is INMATI Its parameter list corresponds to that of INIJIJ with the additional parameter VALUES that stores the numerical quantities 8 CALL INMATi EN NIJ II JI VALUES S A Again the numerical values in VALUES are added to ihoss currently held by the package Mixed use of the routines INAIJI INROWi and SE is SES Thus the user is free to use whatever routine 1s most convenient The same convenience is provided in the input of numerical values for the right hand side vector b SPARSPAK A includes the routine INBI which inputs an id of the right hand side vector CALL INBI 1 VALUE S Here I is the index and VALUE is the numerical value Alternatively the routine INBIBI can be used to input a subvector and its calling seguence is CALL INBIBI NI II VALUES SU aj where NI is the number of input numerical values and II and VALUES are vectors containing the indices and numerical values respectively In both routines incremental caleulation of the numerical values is performed In some situations where the entire right hand E vector is availabl
64. r s Guide gt 30 SPARSPAK A 2 User s Guide 9 Examples In this section we provide several programs which illustrate how SPARSPAK A can be used These programs are derived from the one given in Section 2 1 These examples were run using a standard single precision version of SPARSPAK A under the Berkeley f77 compiler on a DEC VAX 11 780 computer All times reported are in seconds It should be noted that the results will be different if a different version of SPARSPAK A or a different computer is used MI November 1984 l 31 SPARSPAK A User s Guide Example 1 This is an example of the simplest use of SPARSPAK A with each of the modules of method 1 used in sequence The problem solved is a 10X10 symmetric tridiagonal system Axr b where the diagonal elements of A are all 4 and the superdiagonal and subdiagonal elements are all 1 The right hand side vector b is chosen so that the entries of the solution vector r are all ones In the program the nonzero structure of A is input using IJBEGN INIJ and IJEND After ORDRAI is executed the interface modules INAIJ1 and INBI are used to transmit the numerical values of A and b to the package respectively The module SOLVE is called to do the numerical solution and ERESTI is called to compute an estimate of the relative error in the computed solution Then STATSA is called to print out the statistics gathered by the interface during execution Finally the actu
65. rage arrav to execute the ordering routine Response execute SAVEA and restart the computation using ORDRz1 with MAX SA at least as large as that indicated in the error message Insufficient storage in working storage arrav to execute the storage allocation routine The ordering routine was successfullv executed Response same as for error 124 Insufficient storage in working storage arrav to hold the data structure pointers The ordering and storage allocation routines were successfullv executed Response same as for error 124 Insufficient storage in working storage array to hold the numerical values The ordering and storage allocation routines were successfully executed and there was enough storage to hold the data structure pointers Response same as for error 124 INAIJi INROWi INMATi Incorrect execution sequence Probable cause of error routine ORDRz1 was not executed successfully Incompatible input routine Probable cause of error attempt to use input routine INAIJi INROWi or INMAT after using ORDRaj where 1 o 26 SPARSPAK A IERRA 133 134 135 IERRA 141 142 7 3 5 Factorization and EE IERRA 151 152 November 1984 User s Guide INAIJi INROWi INMATi Attempt to input the i j th element of matrix A for i lt j This error occurs only for symmetric matrix methods Le when method _ is odd Methods for symmetric matrices expect elements of the o lower tri
66. rea In such a circumstance the array Si in the user s program would have to be declared in blank common 10 4 Method checking As we discussed in the introduction using a particular method means calling the appropriate interface routines ORDRzx1 INAIJi INROWi INMATI SOLVEi and ERESTI where the last character is a numerical digit denoting the method These ordering input solve and relative error estimation modules cannot be mixed since they in general involve different data structures In order to ensure that these modules are not inadvertently mixed by the user ORDRzi sets the variable METHOD in the control block SPACON equal to 10Xi k where k is an integer that distinguishes orderings A and B This variable is checked bv subsequentiv executed input and solve modules 10 5 Stage sequence checking Another control variable that deserves comment is STAGE As its name implies it is used to keep track of the current step or stage of the execution This variable is particularly important in connection with SAVEA and RSTRTA modules In restarting the system using the RSTRTA routine the variable STAGE in the control block SPACON is restored and it indicates the last successfully completed stage or phase before the routine SAVEA was called In this way the execution can be restarted without repeating already successfully completed steps Another function of this variable is to enforce the correct execution seguence of the various
67. sing binary format b If the subroutines SAVEA and RSTRTA are used then before the user executes his program he must define a file for the FORTRAN logical unit K using the appropriate system control statement or command this depends on the environment in which the program is being executed Furthermore this file must be preserved by the user for later access by the RSTRTA subroutine Thus the user must not write to this file l November 1984 _ 7 a Sii 18 SPARSPAK A W era Guide 5 Solving many problems having the same structure M certain applications many problems which have the same sparsity structure but different numerical values must be solved In this case the structure input ordering and data structure set up needs only to be done once This situation can be accommodated perfectly well by SPARSPAK A The control sequence is depicted by the following flowchart SPRSPK Input Structures Call ORDRzi Input Numerical Values of A and b gt Call SOLVEi When the numerical input routines INAIJi INBI se etc are first called after SOLVEi has been called this is detected by SPARSPAK A and fu computer storage used for A and 6 is initialized to zero Note that if such problems must be solved over an extended Ge period i e in different runs the user can execute SAVEA after executing ORDRzi and thus avoiding he input of the structure of A and the execution of ORDRzi i in subsequent eq
68. store the triangular factors L and U of the permuted matrix PAP In general the ordering and allocation subroutines require different amounts of storage Furthermore their storage requirements are often unpredictable because the number of data structure pointers and the number of nonzeroes in the factors L and U are not known until the subroutines have been executed Thus the interface module ORDRzi may terminate in several distinctly different ways a There was not enough storage to execute the ordering subroutine b The ordering was successfully obtained but there was insufficient storage to initiate execution of the data structure set up storage allocation subroutine c The data structure set up subroutine was executed and the amount of storage required for the data structure pointers etc was determined but there was insufficient storage for these pointers d The data structure was successfullv generated but there is insufficient storage for the actual numerical values so the next step input of the numerical values cannot be executed e ORDRzi was successfully executed and there is EEN storage to face to the next step If any of the above conditions occurs the user may execute SAVEA and re initiate the computation after adjusting his storage declarations either up or down and executing RSTRTA If a or b occurs information is supplied indicating the minimum value of MAXSA needed so that c
69. sults of successfullv completed parts of the computation as shown by the program fragment below CALL ORDRA1 S IF IERRA EQ 0 GO TO 100 CALL SAVEA 8 STOP 100 CONTINUE The variable TERRA is set to the value 10Xk l where 0 lt 1 lt 9 distinguishes the error and k is determined by the type of module that sets IERRA positive ok interface modules 10 save and restart modules SAVEA RSTRTA 11 matrix structure input modules INIJ INIJIJ etc 12 matrix ordering and allocation modules ORDRzi 13 matrix numerical input modules INALJi etc 14 right hand side numerical input modules INBI etc 15 factorization and solution modules SOLVE 16 relative error estimation modules ERESTi November 1984 _ E di 24 SPARSPAK A User s Guide 7 3 1 Save and restart routines IERRA 101 102 103 7 3 2 Input of the matrix structure IERRA 111 112 113 114 IERRA 115 116 117 November 1984 SAVEA RSTRTA Output unit given to SAVEA is not positive Input unit given to RSTRTA is not positive Insufficient storage to restart the computational process The minimum value of MAX SA reguired is printed in the error message INIJ INROW INIJIJ INCLQ Incorrect execution sequence Probable cause of error routine IJBEGN was not executed successfully before 1 3 pairs input began Incorrect execution seguence Probable cause of error routin
70. ted by the racire input modules see Section 7 3 1 2 3 Modules for ordering and storage allocation With an internal representation of the nonzero structure of the matrix A available SPARSPAK A is ready to reorder the matrix problem This is initiated by calling an ordering routine whose name always has the form ORDRzi Here i is a numerical digit between 1 and 6 that signifies the storage method The character x can take values A or B which denotes one of two ordering strategies tailored for storage method 2 Executing the statement CALL ORDRAI1 S will imply the use of storage method 1 and the first ordering algorithm for this method See Section 3 for a discussion of the various methods provided and some guidance on which one to November 1984 S 10 SPARSPAK A User s Guide use Section 8 contains a list of ordering strategies provided by the package The routine ORDRzi not only determines an appropriate ordering for the storage method it sets up the data structure for the reordered matrix problem The package is now ready for numerical input Common Errors Just as in the structure input phase the most common cause of premature termination of the ORDRzi module is insufficient working storage As mentioned above this module performs two functions ordering and storage allocation The ordering step determines the permutation P and the allocation step sets up the appropriate data structures to
71. ure of a row or part of a row is available it is more efficient to use the routine INROW The statement to use is CALL INROW I NIR IR S where I denotes the index of the row under consideration IR is an array containing the column indices of some or all of the nonzeroes in the th row NIR is the number of subscripts in IR and S is the user declared working storage November 1984 be gf M G SPARSP AK A 2 MET User s Guide For example in IR 1 2 IR 2 7 IR S 5 CALL INROW I 3 IR the package is informed of nonzeroes in locations 2 5 5 2 5 5 5 7 and 7 5 of the matrix Note that the column indices in the array IR can be in arbitrary order and the rows can be input in any order c Input of a submatrix structure To provide greater flexibility the package allows the user to input the structure of a submatrix The calling statement is CALL INIJIJ NIJ II JJ S where NIJ is the number of input index pairs and Il and JJ are the arrays containing the row and column indices The following example IT 1 JI 1 II 2 JI 8 II 8 JI 8 CALL a JIJ 8 II JI S informs the package that there are nonzeroes in locations 1 1 1 3 3 1 2 3 and 3 2 d Input of a full submatrix structure The structure of an entire matrix is completely specified if all the full submatrices are given In applications where they are readily available the routine
72. usson solutions November 1984 l M d o 19 SPARSPAK A SESCH n User s Guide 8 Solving many problem which differ only in their right hand side In some applications numerous problems which differ only in their right hand sides must be solved In this case we only want to factor A into LU or LLT once and use the factors repeatedly in the calculation of x for each different b Again SPARSPAK A can handle this situation in a straightforward manner as illustrated by the flowcharts below When SPARSPAK A is used as indicated by flowchart 1 the package detects that no right hand side has been provided during the first execution of SOLVEi and only the factorization is performed In subseguent calls to SOLVEi SPARSPAK A detects that the factorization has already been performed and that part of the SOLVE module is bypassed In flowchart 2 both factorization and solution are performed during the first call to SOLVEi with only the solve part performed in subseguent executions of SOLVEr See Example 3 in Section 9 Note that SAVEA can be used after SOLVE has been executed if the user wants to save the factorization for use in some future calculation November 1984 A a eae IM 20 SPARSPAK A a User s Guide SPRSPK SPRSPK i Input Structure of A Input Structure of A Cal ORDRzi Call ORDRxi Input Numerical Values for A Input Numerical Values for A ei Input Num
73. xample when A is symmetric and positive definite or diagonally dominant In general one can detect when the assumption above is invalid by computing an estimate of the relative error in the computed solution using the subroutine ERESTi see Section 2 6 November 1984 IMMA 4 SPARSPAK A i l User s Guide 2 Modules of SPARSPAK A and how to use them 2 1 User mainline program and an example SPARSPAK A allocates all of its storage from a single one dimensional floating point array which for purposes of discussion we will denote by S In addition the user must provide its size MAXSA which is transmitted to the package via a common block SPAUSR SPARSPAK A USER which has four variables COMMON SPAUSR MSGLVA TERRA MAXSA NEQNS Here MSGLVA is the message level indicator which is used to control the amount of information printed by the package The second variable JERRA is an error code which the user can examine in his mainline program for possible errors detected by the package Detailed discussion of the roles of MSGLVA and IERRA is provided in Section 7 The variable NEQNS is the number of equations The following program illustrates how one might use SPARSPAK A The various subroutines referenced are described in the subseguent parts of this section The problem solved is a 10X10 symmetric tridiagonal system Ar b where the diagonal elements of A are all 4 the superdiagonal and subdiagonal elements are all

SPARSPAK: Waterloo Sparse Matrix Package User's Guide for

Contents

Download Pdf Manuals

Related Search

Related Contents