SPARSPAK: Waterloo Splrse Matrix Package User's Guide for
Contents
1. C INTEGER SUBS 100 INTEGER FILE I ICASE IERRB IPRNTE IPRNTS 1 ISEED J K KGRID KMI MAXINT 1 MAXSB MSGLVB NCOLS NDCONS NDEQNS NSCONS 1 NSEQNS NSUBS OUTPUT ROWNUM TYPE TYPTOL REAL MCHEPS RATIOL RATIOS TIME REAL T 7500 VALUES 100 REAL RESCON RESEON RHS TOL WEIGHT C CO RR OE OE OE OE OE ak ak ak OE OE EO ak ak ak ok k ak ak OE Ok k Ok OE kk OR OE GE OE OE OE k ak k kok k k kk k OE OE OR OR ok kk C COMMON SPKSYS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME COMMON SPBUSR MSGLVB IERRB MAXSB NCOLS NSEONS l NDEQNS NSCONS NDCONS C CK N dk kak kk kk kk kok k ok ke ok ek ok ok ok ok 0k k k k ok k k k k k ok ok ak k k koak ERK EK kk KEK EH EE 9 kk November 1984 31 1 SPK 2SPK 3SPK 4SPK 5SPK 6SPK 7SPK 8SPK 9SPK 10SPK 11SPK 12SPK 13SPK 14SPK 15SPK 16SPK 17SPK 18SPK 19SPK 20SPK 21SPK 22SPK 23SPK 24SPK 25SPK 26SPK 27 SPK 28SPK 29SPK 30SPK 31SPK 32SPK 33SPK SPARSPAK B oo jo je a GOGO ana aaa aaa 1 200 300 400 300 u ME ooo INITIALIZATION CALL SPRSPK FILE I OUTPUT IPRNTS TOL MCHEPS TOL 100 0E0 TOL TYPTOL 1 ISEED 1234567 MSGLVB 2 MAXSB 7500 CALL FILEB FILE NSUBS 4 WEIGHT 1 0E0 TYPE I ROWNUM 0 KGRID 7 KM1 KGRID 1 DO 400 E ol DO 300 J SUBS 1 I 1 KGRID J SUBS 2 I 1 KGRID J I SUB2. eese Ge Ge ee Ge ged ge 6 3 2 Problem initialization subroutine id EE N 6 3 3 Row input subroutine esse eise RA GR EG 0000000 Ge ana 6 3 4 Column ordering subroutine sssssesssnnnnnnznnznnnnnzzannananonantnanzzznnansonnznnezenzantnnzzz 6 3 5 Row ordering subroutine oooooooocononncconononanncnanonononcnnanaconnonoccnonoranaconanancccnos 6 3 6 Solution subroutine si oder od EG A ed ee es 6 3 7 Residual calculation subroutine eed ee ee Gee ee de ee ee is 7 Summary listing of SPARSPAK B interface subroutines ss sssenennennnnzenzaneennnnankznntznzni 8 Examples ee NET OR OPP sn e OVO OBOR NC OM E e O Example Mae Example RE IE N OE det a Example eme p Example M M 9 Reteren et is A hea aie November 1984 SPARSPAK Waterloo Sparse Matrix Package User s Guide for SPARSPAK B A collection of modules to be used with SPARSPAK A for solving sparse constrained linear least squares problems Alan George Esmond Ng Research Report CS 84 37 November 1984 SPARSPAK B User s Guide This document describes the use of a system called SPARSPAK B which is an enhancement to the SPARSPAK A package to allow for the solution of large sparse constrained and unconstrained linear least squares problems SPARSPAK B consists of
3. row ordering is larger than that reguired by the other modules This example also illustrates the use of the save and restart facilities After ORCOLB is called SAVEB is executed to save the current state of the computation After the problem is solved with no row ordering the state after the execution of ORCOLB is restored by executing RSTRTB so that the problem can now be solved without invoking FILEB INXYWB and ORCOLB again Program C SPARSPAK B ANSI FORTRAN RELEASE 111 NAME EX3 C C UNIVERSITY OF WATERLOO JANUARY 1984 CM OE RE ak ROER RE oko ok OE RE RARA ako so kk ER EE RA e kk kk kk k CO RR OE RR ak ak EK OE ak OE RR OE OR RE ER RE EE RR RE ER RE OE OR RR OE OR RE OE OE RR RARA M A I N L I NE PROGRAM EXA KAKE x CEE RR OE OE OE OE OE ak ak ok kak OE OE ak ak ok Ok ak OE OE OE OE REE Ok ORR OE GE OE RE k k RE ER JK EEEE EEE E CPE OE OE RR RR EE ak ak ak ak RR ROK OE OE dk RE OE OE EE k RR EE RR RR RR EG EG RR NE RE EE ER EHD OP AA AA C C EXAMPLE 3 SEE USER S GUIDE C C REOUIREMENTS C AN EXTERNAL FILE FOR UNIT 1 C AN EXTERNAL FILE FOR UNIT 2 C A RANDOM NUMBER GENERATOR RANDOM M Xe INTEGER SUBS 100 INTEGER FILE I ICASE IERRB IPRNTE IPRNTS 1 ISEED J K KGRID KMI MAXINT 1 MAXSB MSGLVB NCOLS NDCONS NDEONS NSCONS 1 NSEQNS NSUBS OPTION ROWNUM SAVE TYPE 1 TYPTOL REAL MCHEPS RATIOL RATIOS TIME REAL T 7000 VALUES 100 R
4. A It is well known that if MT M is sparse the sparsity of its Cholesky factor depends crucially on the symmetric ordering of the rows and columns of MTM 7 Note that a symmetric ordering of MTM is the same as a column ordering of M Thus in order to reduce storage and execution time one should find a good column ordering for M before any numerical computation begins so that the Cholesky factor of the reordered matrix is hopefully sparse This can be achieved after the problem has been supplied to November 1984 9 SPARSPAK B User s Guide 4 the package using INXYWB The column ordering process is invoked by executing the following FORTRAN statement CALL ORCOLB T The algorithm used in SPARSPAK B is an implementation of the minimum degree algorithm due to Liu 14 The module ORCOLB is also responsible for setting up the appropriate data structures for the matrices involved in subsequent numerical computations Common Errors The most common cause of premature termination of the ORCOLB module is insufficient working storage As mentioned above this module performs two functions column ordering and storage allocation The ordering step determines the column permutation and the allocation step sets up the appropriate data structures In general the ordering and allocation subroutines require different amounts of storage Furthermore their storage requirements are often unpredictable because the number of data
5. 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 JA George and J W H Liu Computer solution of large sparse positive definite systems Prentice Hall Inc Englewood Cliffs N J 1981 JA George J W H Liu and E G Y Ng Row ordering schemes for sparse Givens transformations I Bipartite graph model to appear in Linear Algebra and its Appl 1984 J A George J W H Liu and E G Y Ng Row ordering schemes for sparse Givens transformations II Implicit graph model to appear in Linear Algebra and its Appl 1984 J A George and E G Y Ng On row and column orderings for sparse least sguares problems SIAM J Numer Anal 20 1983 pp 326 344 J A George and E G Y Ng On the design and implementation of SPARSPAK B Waterloo sparse constrained linear least sguares package in preparation 1984 M T Heath Some extensions of an algorithm for sparse linear least sguares problems SIAM J Sci Stat Comput 3 1982 pp 223 237 C L Lawson and R J Hanson Solving least squares problems Prentice Hall Inc Englewood Cliffs N J 1974 J W H Liu On multiple elimination in the minimum degree algorithm Technical Report No 83 03 Department of Computer Science York University Downsview Ontario 1983 November 1984 44
6. B SPARSPAK B may provide other information about the computation depending upon the value of MSGLVB the first variable in the common block SPBUSR whether or not errors occur and whether or not the module STATSB is called This section discusses these features of SPARSPAK B Notes Two logical output units ZPRNTS and IPRNTE are required by SPARSPAK B Any information and or statistics which the user has requested are recorded on the output unit IPRNTS while any error messages that might be raised by SPARSPAK B during the execution of the program are recorded on the output unit IPRNTE These two logical output units are set in the initialization module SPRSPK and it is the responsibility of the user and or the computer installation to ensure that the files associated with these two logical units are defined before attempting to execute the program The default output units are both 6 These two output units are also used by SPARSPAK A 6 1 Message level indicator MSGLVB The first variable MSGLVB in the common block SPBUSR 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 MSGLVB is set to one by the user only fatal error messages are printed this option could be useful if SPARSPAK B is being used in the inner loop of a large comp
7. row at a time as shown in the example below CALL INXYWB ROWNUM TYPE NSUBS SUBS VALUES RHS WEIGHT T The parameters of the subroutine are as follows ROWNUM An integer variable associated with each row Usually I lt ROWNUM lt m p and the rows of X and Y are labelled from 1 to m p The rows of X and Y can be intermixed in this labelling It should be emphasized that ROWNUM is used only as a label for a row In some contexts it may be reasonable for several rows to have the same value for ROWNUM The use of ROWNUM along with the use of ORROWB with OPTION set to 5 see Section 2 4 allows the user to impose a specific row ordering TYPE An integer variable having value 1 2 3 or 4 indicating the following November 1984 8 SPARSPAK B User s Guide l a row of A sparse equation in X with corresponding element of y and Wy 2 a row of B dense eguation in X with corresponding element of y and Wx See below for an explanation 3 a row of E sparse constraint in Y with corresponding element of z and Wy 4 a row of F dense constraint in Y with corresponding element of z and Wy See below for an explanation NSUBS An integer variable containing the number of nonzeros in the input row SUBS An integer array containing the column indices subscripts of the nonzeros in the input row VALUES A floating point array containing the numerical values of the nonzeros in the input row in positio
8. structure pointers and the number of nonzeros to be stored are not known until the subroutines have been executed Thus the interface module ORCOLB may terminate in several distinctly different ways a There was not enough storage to execute the column 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 successfully generated but there is insufficient storage for the actual numerical values in the upper trapezoidal matrix so the next step numerical computation cannot be executed e ORCOLB was successfully executed and there is sufficient storage to proceed to the next step If any of the above conditions occurs the user may execute SAVEB and re initiate the computation after adjusting the storage declarations either up or down and executing RSTRTB If a or b occurs information is supplied indicating the minimum value of MAXSB needed so that c or d will occur upon re execution If c occurs the minimum value of MAXSB needed for d is provided When c or d occurs after executing SAVEB adjusting the storage declaration then executing RSTRTB one must again ca
9. the solution process the Euclidean norm of the residual vectors Since the module STATSB can be called at any time some of the above information may not be available and will not be printed Furthermore the amount of information printed also depends on the message level MSGLVB 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 is a useful measure of the amount of arithmetic performed The reader is referred to the examples in Section 8 for more discussion about the output from STATSB 6 3 Error messages ERRB When a fatal error is detected so that the computation cannot proceed a positive code is assigned to IERRB The user can simply check the value of JERRB to see if the execution of the module has been successful This error flag can be used in conjunction with the save and restart feature described in Section 4 to retain the results of successfully completed parts of the computation as shown by the program fragment below CALL ORCOLB T IF IERRB EQ 0 GO TO 100 CALL SAVEB 3 T STOP 100 CONTINUE The variable ZERRB is set to the value 10Xk where 1 lt lt 9 distinguishes the error and k is determined by the type of module that sets JERRB positive k interface modules 20 save and restart modules SAVEB and RSTRTB 2
10. DB and STATSB are called to compute the norm of the residual vectors and to print out the statistics gathered by the package Note that in this example the size of the working storage array T was 7500 while the statistic indicates that the maximum amount of storage required by any modules was only 514 which was the storage requirement in ORCOLB Thus if problems with the same structure are to be solved the user can change the size of T to 514 Program C SPARSPAK B ANSI FORTRAN RELEASE III NAME EX C UNIVERSITY OF WATERLOO JANUARY 1984 Q koX ook ok oko ok ke kokok k ak hak ok ak k ak OO OK Oe kk ae sko k k k k k k ok kokok ok sko se k k oje ok k k kk Xe kk kk ME RR k QC Ok Ok ok kk k kokok Ok ok RE kak ok ok kokok ok k k k k k k kokok k k k k kok kok ok k k k ok ok k ke 0k sko OE k k kok k k 9 Xx MAINLINE PROGRAM ck oko kk kok kk C ok ok ae e ke o kk ok ok kok kx kok k ak k k kok k ok kok Oe k k ak k k ok k k E k k Kk kok kk n Gk ok k sko so k k Xx sko Xx x sko sk C kk k K eG HH x x CK ER RE Ak Sk oe ok ke kk kak ok ok kk ok ok k ok ade ok kok k kok k k ak ck ak k ok OK k ok ok kk k k ok k de k ok kok ok ok k k X GR KK k koko C C EXAMPLE 1 SEE USER S GUIDE C C REQUIREMENTS C AN EXTERNAL FILE FOR UNIT 1 C A RANDOM NUMBER GENERATOR RANDOM C C x oko ok oko ok kk ak akak ok k ok ok ok kk k k ok k k k kok k ok Ok k k ko OE OE EE k k ok kok Oe k k ak k kok kok k OR kk k k Xx oko sk INTEGER
11. DO 200 ICASE 1 4 DO 100 K 1 NSUBS VALUES K RANDOM ISEED 100 CONTINUE ROWNUM ROWNUM I RHS RANDOM ISEED CALL INXYWB ROWNUM TYPE NSUBS SUBS 1 VALUES RHS WEIGHT T 200 CONTINUE 300 CONTINUE 400 CONTINUE m oo CALL ORCOLB T m m MEME LM om mm vm mm Sm EE EE ME EE Ma EE EE MA mm EE lS CALL SAVEB SAVE T OPTION 0 300 CONTINUE CALL ORROWB OPTION T SS mm EE am ER ME am w m m EE m m a m a m m m a OM Pew we m EE m w m a a e w o o w am w am am am a m am am ou ew m mm e am am am m en a w m m ME et EE EE am w en m wm m m CALL LSQSLV TOL TYPTOL T CALL RESIDB RESEON RESCON T November 1984 User s Guide 36 44SPK 45SPK 46SPK 47SPK 48SPK 49SPK 50SPK 51SPK 52SPK 53SPK 54SPK 55SPK 56SPK 57SPK 58SPK 59SPK 60SPK 61SPK 62SPK 63SPK 64SPK 65SPK 66SPK 67SPK 68SPK 69SPK 70SPK 71SPK 72SPK 73SPK 74SPK 75SPK 76SPK 77SPK 78SPK 79SPK 80SPK 81SPK 82SPK 83SPK 84SPK 85SPK 86SPK 87SPK 88SPK 89SPK 90SPK 91SPK 92SPK 93SPK 94 SPK 95SPK 96SPK 97SPK 98SPK 99SPK 100SPK 101SPK 102SPK 103SPK 104SPK 105SPK 106SPK 107SPK 108SPK 109SPK SPARSPAK B C PRINT STATISTICS GATHERED BY SPARSPAK B Ce 1 om RR REED a Sa TEC qu ai k se ad AO RE a k AQ Ee v Set CALL STATSB IF OPTION EQ 4 STOP C C EET ETE C RESTORE STATE OF COMPUTATION CU hee EE EE DEAD a CALL RS
12. E ARRAY T C T LTD N et STNA et vee dt a i n VO S EA AF EA WRITE OUTPUT 11 T K K 1 NCOLS 11 FORMAT 10H SOLUTION 1PS5EI35 5 C at ere RR fm me ai AE SEER ES J Me he da DE ie Tam pups DA LAE Ba IDEAT C PRINT STATISTICS GATHERED BY SPARSPAK B C menu wu uu AA BRIL AG NS ER LN OU la de d ESE a OE Hem D CALL STATSB C STOP END Output CO UNIVERSITY OF WATERLOO EKER KERKE SPARSE MATRIX PACKAGE k ohh kk kok k k S P A R S P A K hak sko sko sko sk KK K RELEASE 3 d XOU EEA C JANUARY 1984 x KERKE ER ANS FORTRAN wok SINGLE PRECISION KERR KKRE LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES 6 OUTPUT UNIT FOR STATISTICS 6 FILEB FILE INITIALIZATION INXYWB INPUT ROWS ORCOLB FIND COLUMN ORDERING LSOSLV LEAST SOUARES SOLVE RESIDB COMPUTE RESIDUAL SOLUTION 4 27621e 01 3 45408e 02 1 26385e 01 3 57654e 01 3 1 31271e 01 6 67797e 01 3 41631e 02 4 47427e 01 4 3 17158e 03 1 02797e 01 3 66493e 01 7 21199e 01 3 3 04172e 01 9 93229e 02 5 65309e 01 4 02507e 01 5 6 21579e 02 2 51969e 01 5 16541e 01 1 59730e 01 6 3 53588e 02 4 96227e 01 3 99947e 01 3 46664e 01 3 November 1984 33 100SPK 1I01SPK 102SPK 103SPK 104SPK 105SPK 106SPK 107SPK 108SPK 109SPK 110SPK 111SPK 112SPK 113SPK 114SPK 115SPK 116SPK 117SPK 118SPK 119SPK I20SPK 121SPK 122SPK 123SPK 124SPK 125SPK 126SPK 13655e 01 08074e 01 65426e 01 26985e 02 74474e 02 04207e 01 SPARSPAK B ta SA N NUMBER NU
13. E k EE Xx Xx RE x x C xk k kk kk kk kk kok C C EXAMPLE 4 SEE USER S GUIDE C C REQUIREMENTS C AN EXTERNAL FILE FOR UNIT 1 C AN EXTERNAL FILE FOR UNIT 2 C A RANDOM NUMBER GENERATOR RANDOM C CE RE kk OR Rd dk kk kk ak ok kak ok ok ok ok ok ak kk ER RR kk hak ak kk k k k kok oko sko ER EE Ok OK kk ah kk OE k k sk je X k ok k KEK C INTEGER FILE IERRB IPRNTE IPRNTS MAXINT MAXSB 1 MSGLVB NCOLS NDCONS NDEONS NSCONS NSEONS 1 OPTION SAVE TYPTOL REAL MCHEPS RATIOL RATIOS TIME REAL T 10000 REAL RESCON RESEON TOL C C A RK NE RR kk kk k k kk OR k k kok kk kk kk k k OE kk kk kk REE X C COMMON SPKSYS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME COMMON SPBUSR MSGLVB IERRB MAXSB NCOLS NSEONS 1 NDEQNS NSCONS NDCONS C CA k k ok ok ok OE ok kk ok Gk Ok k ok Gk Gk ok kk Ok ak k Ok Ok RO k Ok Ok k k ok k k Ok kk OR k kk kk kk X Xe OE RR sko k k C EC ees hoe d C INITIALIZATION C oss EE EU WS gre CALL SPRSPK C FILE I SAVE 2 OPTION 4 TOL MCHEPS TOL 100 0E0 TOL TYPTOL I C MSGLVB 2 MAXSB 10000 C CALL RSTRTB SAVE T C CE o Jie Da C ORDER ROWS CU v SANA CALL ORROWB OPTION T C C Oe A EE RE EO A AE C COMPUTE LEAST SOUARES SOLUTION AND RESIDUALS o AE IE OE o AE EA a n E OE O GE AE DONE IE a EI Oe CALL LSOSLV TOL TYPTOL T CALL RESIDB RESEON RESCON T C November 1984 User s Guide 42 1
14. EAL RESCON RESEON RHS TOL WEIGHT C Q o5 ok ook o o Gk kk kkk k kk kk aa alada all Gk GR ko akok x x Xx 0 X 6 G3 0 k 0 kk kkk k ok sok Xx 9x ak le 0 00 C COMMON SPKSYS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME COMMON SPBUSR MSGLVB IERRB MAXSB NCOLS NSEONS 1 E NDEONS NSCONS NDCONS CC 5 5 8 oo ok OE OE OE ok O OE EE obo OE ok ob ok ok kk ME kk ak ok X x 03 ak ak k G0 ak ak ak k k oak ok k k k ok k k RE RR kok ER k k kk kok k C G OE EE DK C INITIALIZATION QU NA E E IN CALL SPRSPK C FILE SAVE 2 November 1984 35 1 SPK 2SPK 3SPK 4 SPK 5SPK 6 SPK 7 SPK 8 SPK 9SPK 10SPK 11SPK 12SPK 13SPK 14SPK 15SPK 16SPK 17SPK 18SPK 19SPK 20SPK 21SPK 22SPK 23SPK 24SPK 25SPK 26SPK 27SPK 28SPK 29SPK 30SPK 31SPK 32SPK 33SPK 34SPK 35SPK 36SPK 37SPK 38SPK 39SPK 40SPK 41SPK 42SPK 43SPK SPARSPAK B MOKO N GO AAADA aA COLO 00 6 66 AAA C TOL MCHEPS TOL 100 0E0 TOL TYPTOL l ISEED 1234567 MSGLVB 2 MAXSB 7000 CALL FILEB FILE m en aw n ak AW a am mm EE MR san Mm EE ME a en aw an sa ewo EE EE N Rm NE NSUBS 4 WEIGHT 1 0E0 TYPE I ROWNUM 0 KGRID 10 KMI KGRID 1 DO 400 oe od DO 300 J I KMI L mm EE va c mo e m m EE mm M SUBS 1 I 1 KGRID J SUBS 2 I 1 KGRID J 1 SUBS 3 I KGRID J SUBS 4 I KGRID J I GENERATE NUMERICAL VALUES USING A RANDOM NUMBER GENERATOR
15. EX2A C C UNIVERSITY OF WATERLOO JANUARY 1984 CK k EE RRR ERROR RR OK Ok OK Ok Ok Ok OR Ok ak ak ak ak ORO OE ak ak ak k k Ok k OE k k k E k RE RR EE ak sko ak ok k k k ak ok k ak k x 0 k C ER ak k ok ak ak ak ak ak Ok OE Ok OE Ok OE ok Ok Ok k OE ak k GO k Ok EE ak OO k OE Ok OE ok k k Gk Ok k kok kok OE EE ER kok k OE OE E e sko ab CO RR ka MAINLINE PROGRAM TEA EH KH CER ERR kok ak Ok Ok OK ak ak ak ak k OK Gk ER RR OE OO Gk OE k ok k k k EE OE OE k kok Ob RR RR k k kok kokok ER EE KE QC OE OE OR OE EEE k k k ak k kok ak EK Gk Gk ok OE k OE Ok EE G OE Gk Gk OE OE Ok Gk G G Ob k ER dk k kok E kk Oe Ok Ok OO A OR M k C C EXAMPLE 2 SEE USER S GUIDE C C REQUIREMENTS C AN EXTERNAL FILE FOR UNIT 1 C A RANDOM NUMBER GENERATOR RANDOM C Q Aaa o kak xXk ok sk ok ok ok ok ok ok ok ak ak akk k k k k Ok OR OK ok ok ok kk ak h ok ok ok kk k k kk kk k k kk k Gk EK HE KH x INTEGER SUBS 100 INTEGER FILE 1 ICASE IERRB IPRNTE IPRNTS 1 ISEED J K KGRID KMI MAXINT 1 MAXSB MSGLVB NCOLS NDCONS NDEQNS NSCONS 1 NSEQNS NSUBS OUTPUT ROWNUM TYPE TYPTOL REAL MCHEPS RATIOL RATIOS TIME REAL T 7500 VALUES 100 REAL RESCON RESEQN RHS TOL WEIGHT C COR Eo Ok ok ak ok OE ak ak OK k OK ok k ak k ak k ak ak k kok k k ak ak ok ak ak kok kok k k k k k k k OE OE OE k k kk kk kk OE k k C COMMON SPKSYS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME COMMON SPBUSR MSGLVB IERRB MAXSB NCOLS NSEO
16. JANUARY 1984 FILEB OUTPUT UNIT FOR ERROR MESSAGES OUTPUT UNIT FOR STATISTICS FILE INITIALIZATION INXYWB INPUT ROWS ORCOLB FIND COLUMN ORDERING ORROWB FIND ROW ORDERING EMSGB SYSTEM B ERROR ORROWB ERROR NUMBER INSUFF SPACE FOR ROW ORDERING MAXSB MUST AT LEAST BE SAVEB SAVE STORAGE VECTOR STATSB SYSTEM B STATISTICS November 1984 SIZE OF STORAGE ARRAY MAXSB NUMBER OF COLUMNS UNKNOWNS NUMBER OF SPARSE EOUATIONS NUMBER OF DENSE EOUATIONS NUMBER OF SPARSE CONSTRAINTS NUMBER OF DENSE CONSTRAINTS TIME FOR COLUMN ORDERING STORAGE FOR COLUMN ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION TOTAL TIME REOUIRED MAXIMUM STORAGE REOUIRED il l User s Guide 244 5850 3000 100 324 0 217 2269 0 050 1966 0 267 2269 41 115SPK 116SPK 117SPK 118SPK 119SPK SPARSPAK B Program 2 C SPARSPAK B ANSI FORTRAN RELEASE III NAME EX4B C C UNIVERSITY OF WATERLOO JANUARY 1984 Ci k rka ki kok Gk ok GE ok ok kk kok ok Rok hr ko k k k koak ak k RR Ok k akok k k AE kk a a kok REE KARR kk kk k k OR x k X Ct k kk kk kk kk kkk ok kokk kok k kk kk akok ak kok ak kokok hak kak ok k k k kk kk k k k k koko kk kk kk sk MA INLINE PROGRAM CPE tk EE EE RE RR kk ko kok k k ok kk dok xxx k hak ak k kk RE k Ok kk OR RR kk ana kk k k k k sko oko kk KK Xx ak Q Ek RR dk ER oko ok RE ok ok ok k ok Ak kok Sko ok ok ok k k Gk OK kk EE RE kk kokok k kok k k k k k kk kk ER O
17. LEAST SOUARES SOLUTION AND RESIDUALS 117SPK CA o 118SPK CALL LSQSLV TOL TYPTOL T 119SPK CALL RESIDB RESEON RESCON T 120SPK C 121SPK Co AAA ER AE EE AA OE Se alee 122SPK C PRINT THE SOLUTION FOUND IN THE FIRST NCOLS 123SPK C LOCATIONS IN THE WORKING STORAGE ARRAY T 124SPK C krie EE a B 8 EE a EN Pre EO N 125SPK WRITE OUTPUT 11 T K K lt l NCOLS 126SPK 11 FORMAT 10H SOLUTION 1PSE15 5 127SPK C MEI ETER o on e NE m 128SPK C PRINT STATISTICS GATHERED BY SPARSPAK B 129SPK C ee mmm mm umu m 130SPK CALL STATSB 131SPK C 132SPK STOP 133SPK END 134SPK Output x E F UNIVERSI TY OF WATERLOO KERE RR SPARSE MATRIX PACKAGE EK he kok Ok kk k k S PAR S PAK RELEASE 3 C JANUARY 1984 REK RKK RE ANSI FORTRAN EER KREET SINGLE PRECISION x CE LAST UPDATE JANUARY 1984 dk Sko kk k sk kkk k k k KKK EK KK k OUTPUT UNIT FOR ERROR MESSAGES 6 OUTPUT UNIT FOR STATISTICS 6 FILEB FILE INITIALIZATION November 1984 26 SPARSPAK B INXYWB INPUT ROWS ORCOLB FIND COLUMN ORDERING LSQSLV LEAST SQUARES SOLVE RESIDB COMPUTE RESIDUAL SOLUTION 2 00000e 00 5 00000e 00 2 64674e 00 3 00018e 00 1 18907e 00 2 28644e 01 3 82073e 01 1 95622e 00 3 16640e 02 5 77518e 01 3 68422e 01 9 56243e 02 6 70207e 01 5 10281e 01 3 93671e 01 STATSB SYSTEM B STATISTICS SIZE OF STORAGE ARRAY MAXSB NUMBER NUMBER NUMBER NUMBER NUMBER TIME FOR ALLOCATION STORAGE FOR ALLOCATION TIME FOR ROW ORDERING ST
18. MBER NUMBER NUMBER NUMBER TIME FOR COLUMN STORAGE FOR COLUMN ORDERING OF OF OF OF OF TIME FOR ALLOCATION STORAGE FOR ALLOCATION TIME FOR ROW ORDERING COLUMNS UNKNOWNS SPARSE EOUATIONS 86252e 01 3 20503e 01 5 32275e 01 90001e 01 4 31955e 01 1 54122e 01 82577e 01 74 59609e 01 1 64609e 00 53660e 01 1 02709e 03 7 23423e 01 STATSB SYSTEM B STATISTICS SIZE OF STORAGE ARRAY MAXSB DENSE EOUATIONS SPARSE CONSTRAINTS DENSE CONSTRAINTS ORDERING STORAGE FOR ROW ORDERING TIME FOR SOLUTION OPERATION COUNT FOR SOLUTION STORAGE FOR SOLUTION TIME FOR COMPUTING RESIDUAL OPN COUNT FOR COMPUTING RESIDUAL STORE FOR COMPUTING RESIDUAL TOTAL TIME REOUIRED MAXIMUM STORAGE REOUIRED RESIDUAL IN EOUATIONS RESIDUAL IN CONSTRAINTS l l li mm N VO t User s Guide 22484e 01 86543e 02 08693e 01 13268e 01 7500 144 0 1066 0 902 0 0 400 128442 1132 0 915 985 850 1132 l 1 117 317 000 3 663e 00 0 000e 01 4 58985e 01 3 88643e 03 6 20310e 02 November 1984 34 SPARSPAK B User s Guide Example 3 The effect of row ordering on execution time is illustrated The problem being considered is the unconstrained problem defined on a k Xk grid Note that when row ordering option OPTION is set to four there is a reduction in execution time for LSOSLV However one should also note that the amount of storage reguired to perform
19. MNS UNKNOWNS 49 NUMBER OF SPARSE EQUATIONS 145 NUMBER OF DENSE EQUATIONS 0 NUMBER OF SPARSE CONSTRAINTS 0 NUMBER OF DENSE CONSTRAINTS 0 TIME FOR COLUMN ORDERING 0 050 STORAGE FOR COLUMN ORDERING 5146 TIME FOR ALLOCATION 0 200 STORAGE FOR ALLOCATION 2795 TIME FOR ROW ORDERING 0 STORAGE FOR ROW ORDERING 0 TIME FOR SOLUTION 2 550 OPERATION COUNT FOR SOLUTION 248739 STORAGE FOR SOLUTION 1814 TIME FOR COMPUTING RESIDUAL 0 367 OPN COUNT FOR COMPUTING RESIDUAL 915 000 STORE FOR COMPUTING RESIDUAL 1667 TOTAL TIME REQUIRED 3 167 MAXIMUM STORAGE REQUIRED 5146 RESIDUAL IN EQUATIONS 3 663e 00 RESIDUAL IN CONSTRAINTS 0 000e 01 Program 2 C SPARSPAK B ANSI FORTRAN RELEASE III NAME EX2B C UNIVERSITY OF WATERLOO JANUARY 1984 Qoo ok ob RR EE EE A Ge kk EE ke dk Gk Gk ok oko k okok Gk e GR e Gb EE OE ae XE RE E GO sk X oak k ak k RE kk OE OR Xx X 0 kk CHHHH AH kk hak OE Ok ok oO OE ok ok kk ER kk OC Ok OE OE X ok Ok Ok kok kk X OE MR AE OE oko ok ok oko kok ok kk EE KEK 0 s CK oko ok x x x x n MAINLINE P RO G RAM kkk sko it sk k kkk CE o5 oko KEER RE EER x G0 n G0 GB G9 95 G0 Xx dod GR GbR kokk kok GR OR A e GR G0 EG kok RR A G0 kok 0 0X RE 0 0x 0X kok k x ah CF koh ok ok ok o Ok kokk dk oe kk hak aha k kok kk k skos sko kok A 3 XA ERE EEE ES C C EXAMPLE 2 SEE USER S GUIDE C C REQUIREMENTS C AN EXTERNAL FILE FOR UNIT 1 C A RANDOM NUMBER GENERATOR RANDOM
20. NS 1 NDEONS NSCONS NDCONS C ook ok kk dk o x x kk hak kk k k kk dk oe kk ok k ok ok ok ok Ok ok k kok ok k k k hi oe k ok ok k k k k ok k k ok kk sko sko koko k OE RK KK OOGOGO aqq MM me ban A M N AN my moun O z mom November 1984 28 ISPK 2SPK 3SPK 4SPK 5SPK 6SPK 7SPK 8SPK 9SPK 10SPK 1ISPK 12SPK 13SPK 14SPK I 5SPK I6SPK 17SPK 18SPK 19SPK 20SPK 21SPK 22SPK 23SPK 24 SPK 25SPK 26 SPK 27SPK 28SPK 29SPK 30SPK 31SPK 32SPK 33SPK 34SPK 35SPK 36SPK 37SPK SPARSPAK B O QUO O 0 BED DO LO OOGOG 100 l 200 300 400 CONTINUE 500 CALL SPRSPK FILE OUTPUT I PRNTS TOL MCHEPS TOL 100 0E0 TOL TYPTOL ISEED MSGLVB MAXSB 1 1234567 2 7500 CALL FILEB FILE KGRID w om m MDO omo wm EE m et mo ME am w an EE am en am a n am a sm EE SE ME w EE aw nn am am m am m a n n 7 KMI KGRID I bo 400 DO SUBS 1 I1 1 KGRID J SUBS 2 I 1 KGRID J 1 SUBS 3 I KGRID J SUBS 4 I KGRID J 1 GENERATE NUMERICAL VALUES USING A RANDOM NUMBER GENERATOR DO 200 ICASE I 4 DO 100 K 1 NSUBS User s Guide VALUES K RANDOM ISEED CONTINUE ROWNUM ROWNUM 1 RHS RANDOM ISEED CALL INXYWB ROWNUM TYPE NSUBS SUBS VALUES RHS CONTINUE CONTINUE m mo om mm a ae m Mm n OM m mo ROWNUM ROW
21. NUM I TYPE NSUBS KGRID KGRID DO 500 I 1 NSUBS SUBS I I VALUES I 1 0E0 CONTINUE RHS 7 0E0 WEIGHT 1 0E0 CALL INXYWB ROWNUM TYPE NSUBS SUBS November 1984 RAS WEIGHT T m WEIGHT T VALUES 29 38SPK 39SPK 40SPK 41SPK 42SPK 43SPK 44SPK 45SPK 46SPK 47SPK 48SPK 49SPK SOSPK SISPK 52SPK 53SPK 54SPK 55SPK 56SPK 37SPK 58SPK 59SPK 60SPK 6 ISPK 62SPK 63SPK 64SPK 65SPK 66SPK 67SPK 68 SPK 69 SPK 70SPK 71SPK 72SPK 73SPK 74SPK 75SPK 76SPK 77SPK 78SPK 79SPK 80SPK 81SPK 82SPK 83SPK 84 SPK 85SPK 86 SPK 8 7 SPK 8 8SPK 8 9SPK 9OSPK 91SPK 92SPK 93SPK 94SPK 95SPK 96SPK 97SPK 98SPK 99SPK 100SPK 101SPK 102SPK 103SPK SPARSPAK B User s Guide C ORDER COLUMNS I n dejes ee eee CALL ORCOLB T C C emcee Haas tis as A RP EL ML P M PU MET PENES C COMPUTE LEAST SQUARES SOLUTION AND RESIDUALS C Ir IDCM RUM S Be es CALL LSQSLV TOL TYPTOL T CALL RESIDB RESEQN RESCON T C C ME NE uel od aa oy ot fas DAL US EJ RE E NA Mads EN Hee g ASA Va bo nan RN a ete tn tk po Sak C PRINT THE SOLUTION FOUND IN THE FIRST NCOLS C LOCATIONS IN THE WORKING STORAGE ARRAY T C TE ag SENSE TAN N KEN ENE a MA N be ee ie a N EKT EA TER WRITE OUTPUT 11 T K K 1 NCOLS 11 FORMAT 10H SOLUTION 1P5E15 5 C a JV m n Yeas Lahm ee A e DS EE BEES et te echt bee no er Se ED CA a onak ae Jes C PRINT STATISTICS GATHERED BY SPARS
22. NUMERICAL VALUES USING A RANDOM NUMBERGENERATOR mem m m mo oma ab A um ee am de am an em am m wm m aw mem w DO 200 ICASE I 4 DO 100 K 1 NSUBS VALUES K RANDOM ISEED CONTINUE ROWNUM ROWNUM RHS RANDOM ISEED CALL INXYWB ROWNUM TYPE NSUBS VALUES RHS CONTINUE CONTINUE 400 CONTINUE ROWNUM TYPE November 1984 3 ROWNUM 1 SUBS T 25 27S PK 28SPK 29SPK 30SPK 31SPK 32SPK 33SPK 34SPK 35SPK 36SPK 37SPK 38S PK 39S PK 40SPK 41SPK 42SPK 43SPK 44S PK 45S PK 46SPK 47SPK 4 8 S PK 49SPK 50SPK 51SPK 52SPK 53SPK 54SPK 55SPK 56SPK 57SPK 58SPK 59SPK 60SPK 61SPK 62SPK 63SPK 64SPK 65SPK 66SPK 67SPK 68SPK 69SPK 70 SPK 71SPK 72SPK 73SPK 74SPK 75SPK 76SPK 77SPK 78SPK 79S PK 80SPK 81SPK 82SPK 83SPK 84SPK 85SPK 86SPK 87SPK 88SPK 89SPK 90SPK 91SPK 92SPK SPARSPAK B User s Guide NSUBS 1 93SPK SUBS 1 1 94SPK VALUES 1 1 0E0 95SPK RHS 2 0E0 96SPK WEIGHT 1 0E0 97SPK CALL INXYWB ROWNUM TYPE NSUBS SUBS VALUES 98SPK RHS WEIGHT T 99SPK C 100SPK ROWNUM ROWNUM 1 101SPK TYPE 3 102SPK NSUBS 1 103SPK SUBS 1 2 104SPK VALUES 1 1 0E0 105SPK RHS 5 0E0 106SPK WEIGHT 1 0E0 107SPK CALL INXYWB ROWNUM TYPE NSUBS SUBS VALUES 108SPK RHS WEIGHT T 109SPK 110SPK A 111SPK C ORDER COLUMNS 112SPK C 113SPK CALL ORCOLB T 114SPK C 115SPK N EE EE EE EE OE AO TE su 116SPK C COMPUTE
23. NXYWB was not executed successfully 232 Incorrect execution seguence Probable cause of error routine ORCOLB was called after having already been executed successfully November 1984 19 SPARSPAK B IERRB 233 234 235 236 237 238 User s Guide ORCOLB Insufficient storage in working storage array to create adjacency structure The minimum value of MAXSB required is printed in the error message Response execute SAVEB and restart the computation using ORCOLB with MAXSB at least as large as that indicated in the error message Number of variables or columns NCOLS is zero Number of equations and constraints is zero i e m p 0 Insufficient storage in working storage array to execute the column ordering routine The minimum value of MAXSB required is printed in the error message Response execute SAVEB and restart the computation using ORCOLB with MAXSB at least as large as that indicated in the error message Insufficient storage in working storage array to execute the storage allocation routine The column ordering routine was successfully executed Response same as for error 236 Insufficient storage in working storage array to hold the data structure pointers The column ordering and storage allocation routines were successfully executed Response same as for error 236 6 3 5 Row ordering subroutine IERRB 241 242 243 244 ORROWB Incorrect execution sequence Probabl
24. O a a mo a o a m m e m w m w a m m U CALL STATSB Or Oy 00 STOP END Note If the SPARSPAK B package available to you is a double precision version the REAL declarations 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 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 shown below CALL SPRSPK In fact SPRSPK is the initialization routine for the entire SPARSPAK package including both SPARSPAK A and SPARSPAK B Note that the only variable in the common block SPBUSR that must be explicitly assigned a value by the user is MAXSB although he may wish to set others as well such as MSGLVB It is assumed that the subroutines which comprise SPARSPAK A and SPARSPAK B have been compiled into a ibrary 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 and SPARSPAK B November 1984 7 SPARSPAK B User s Guide SPARSPAK B reguires an external seguential file for th
25. ON OPERATION COUNT FOR SOLUTION STORAGE FOR SOLUTION TIME FOR COMPUTING RESIDUAL OPN COUNT FOR COMPUTING RESIDUAL STORE FOR COMPUTING RESIDUAL TOTAL TIME REOUIRED MAXIMUM STORAGE REOUIRED RESIDUAL IN EOUATIONS RESIDUAL IN CONSTRAINTS November 1984 l l i ii 10000 100 324 0 2269 0 1966 0 5850 Ja 324939 2306 0 1944 2006 5 5850 User s Guide 217 700 000 000 4 823e 00 0 000e 01 43 6 2SPK 63SPK 64SPK 65SPK 66SPK 67SPK 68SPK SPARSPAK B User s Guide 9 References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 A Bjorck A general updating algorithm for constrained linear least sguares problems SIAM J Sci Stat Comput 5 1984 pp 394 402 E C H Chu J A George J W H Liu and E G Y Ng User s guide for SPARSPAK A A collection of modules for solving sparse systems of linear equations in preparation 1984 J A George and M T Heath Solution of sparse linear least sguares problems using Givens rotations Linear Algebra and its Appl 34 1980 pp 69 83 J A George M T Heath and E G Y Ng A comparison of some methods for solving sparse linear least sguares problems SIAM J Sci Stat Comput 4 1983 pp 177 187 J A George M T Heath and E G Y Ng Solution of sparse underdetermined systems of linear eguations to appear in SIAM J Sci Stat Comput 1984 J
26. ORAGE FOR ROW ORDERING TIME FOR SOLUTION OPERATION COUNT FOR SOLUTION STORAGE FOR SOLUTION TIME FOR COMPUTING RESIDUAL COLUMNS UNKNOWNS SPARSE EOUATIONS DENSE EOUATIONS SPARSE CONSTRAINTS DENSE CONSTRAINTS TIME FOR COLUMN ORDERING STORAGE FOR COLUMN ORDERING td M t vod odW OPN COUNT FOR COMPUTING RESIDUAL STORE FOR COMPUTING RESIDUAL TOTAL TIME REQUIRED MAXIMUM STORAGE REQUIRED RESIDUAL IN EQUATIONS RESIDUAL IN CONSTRAINTS November 1984 ll bm ba o bo User s Guide 11532e 00 77276e 01 93660e 01 84502e 01 24445e 02 7500 0 017 514 0 033 425 0 0 0 650 21758 7 451 0 133 390 000 376 0 833 514 900e 00 0 000e 01 27 95429e 01 61121e 03 61554e 01 24321e 00 10439e 00 SPARSPAK B Users Guide Example 2 The purpose of this example is to illustrate the effect of dense rows The problem being solved is the unconstrained linear least sguares problem min By Il x a Here A is the matrix defined on a k XK grid and B and yg are given below B 111 11 ys 7 where X is partitioned into two portions In the first run all rows of X are treated as sparse TYPE 1 In the second run the rows of A are regarded as sparse TYPE 1 whereas the row of B is treated as dense TYPE 2 Note the difference in storage requirements and execution times in the output Program 1 C SPARSPAK B ANSI FORTRAN RELEASE III NAME
27. PAK B C TE EN cg A A EO KO EE DOE TM CALL STATSB C STOP END Output UNIVERSITY OF WATERLOO xxx tz tt SPARSE MATRIX PACKAGE kok OR x e ER S P A R S P A K kk k k KK kkk ik RELEASE 3 TERRE k CJ JANUARY 1984 sak ANSI FORTRAN HeEEEKERESX SINGLE PRECISION o LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES 6 OUTPUT UNIT FOR STATISTICS 6 FILEB FILE INITIALIZATION INXYWB INPUT ROWS ORCOLB FIND COLUMN ORDERING LSOSLV LEAST SOUARES SOLVE RESIDB COMPUTE RESIDUAL SOLUTION 4 27624e 01 3 45387e 02 1 26388e 01 3 57653e 01 3 1 31271e 01 6 67795e 01 3 41637e 02 4 47426e 01 4 3 17152e 03 1 02795e 01 3 66493e 01 7 21198e 01 3 3 04172e 01 9 93216e 02 5 65310e 01 4 02508e 01 5 6 21576e 02 2 51969e 01 5 16541e 01 1 59729e 01 6 3 53609e 02 4 96227e 01 3 99947e 01 3 46666e 01 3 2 86251e 01 3 20507e 01 5 32274e 01 2 22485e 01 4 6 90003e 01 4 31956e 01 1 54123e 01 9 86527e 02 3 1 82578e 01 4 59611e 01 1 64609e 00 7 08693e 01 6 5 53662e 01 1 02702e 03 7 23422e 01 1 13267e 01 November 1984 30 104SPK 105SPK 106SPK 107SPK 108SPK 109SPK 110SPK 111SPK 112SPK 113SPK 114SPK 115SPK 116SPK 117SPK 118SPK 119SPK 120SPK 121SPK 122SPK 123SPK 124SPK 125SPK 126SPK 13657e 01 08073e 01 65426e 01 26987e 02 74453e 02 04208e 01 58984e 01 88556e 03 20304e 02 SPARSPAK B User s Guide STATSB SYSTEM B STATISTICS SIZE OF STORAGE ARRAY MAXSB 7500 NUMBER OF COLU
28. RESIDB Incorrect execution sequence Probable cause of error routine LSOSLV was not executed successfully Insufficient storage in working storage array to compute residuals The minimum value of MAXSB reguired is printed in the error message Response execute SAVEB and restart the computation using RESIDB with MAXSB at least as large as that indicated in the error message 21 SPARSPAK B User s Guide 7 Summary listing of SPARSPAK B interface subroutines SPARSPAR SPRSPE initialization Probl initialization INABWB ROWNUM TYPE NSUBS SUBS VALUES RHS WEIGHT T Col ordering Row ORROWB OPTION T ordering Numerical LSOSLV TOL TYPTOL T solution Residual esidual gESIDB RESEON RESCON T calculation statistics Save SAVEB K T Restart RSTRTB K T November 1984 22 SPARSPAK B User s Guide 8 Examples In this section we provide several examples which illustrate how SPARSPAK B can be used The programs were compiled using the Berkeley f77 compiler and run on a DEC VAX 11 780 computer A single precision version of SPARSPAK A and SPARSPAK B was used All times reported are in seconds It should be noted that the results may be different if a different version of SPARSPAK A and SPARSPAK B are used or if the programs are compiled and run on a different computer The sample problems are variants of the following unconstrained linear least sguares problem min XB la g
29. ROWB LSOSLV RES IDB STATSB RESTART SYSTEM B FIND ROW ORDERING LEAST SOUARES SOLVE COMPUTE RESIDUAL SYSTEM B STATISTICS SIZE OF STORAGE ARRAY MAXSB NUMBER OF COLUMNS UNKNOWNS NUMBER OF SPARSE EOVATIONS NUMBER OF DENSE EOUATIONS NUMBER OF SPARSE CONSTRAINTS NUMBER OF DENSE CONSTRAINTS TIME FOR COLUMN ORDERING STORAGE FOR COLUMN ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION TIME FOR ROW ORDERING STORAGE FOR ROW ORDERING TIME FOR SOLUTION OPERATION COUNT FOR SOLUTION STORAGE FOR SOLUTION TIME FOR COMPUTING RESIDUAL OPN COUNT FOR COMPUTING RESIDUAL STORE FOR COMPUTING RESIDUAL TOTAL TIME REQUIRED MAXIMUM STORAGE REQUIRED RESIDUAL IN EQUATIONS RESIDUAL IN CONSTRAINTS l l l l User s Guide 0 717 1944 000 2006 5 667 2306 4 823e 00 0 000e 01 0 217 2269 0 050 1966 0 083 5850 3 983 325856 2306 0 717 1944 000 2006 5 050 5850 4 823e 00 0 000e 01 November 1984 38 SPARSPAK B User s Guide Example 4 This example illustrates the use of the save and restart facilities to handle errors detected by SPARSPAK B The problem being solved is the same as the one in Example 3 and row ordering option is set to four The size of the working storage array is initially set to 3000 which will be insufficient for ORROWB as illustrated by Example 3 The state of the computation is saved by calling SAVEB when the error is detected after ORROWB is called After adjusting the si
30. S 3 I KGRID J SUBS 4 I KGRID J 1 GENERATE NUMERICAL VALUES USING A RANDOM NUMBER GENERATOR DO 200 ICASE 1 4 DO 100 K 1 NSUBS VALUES K RANDOM ISEED CONTINUE ROWNUM ROWNUM 1 RHS RANDOM ISEED CALL CONTINUE CONTINUE CONTINUE ROWNUM ROWNUM 1 TYPE 2 NSUBS KGRID KGRID DO 500 I 1 NSUBS SUBS I I VALUES I 1 0E0 CONTINUE RHS 7 0E0 WEIGHT 1 0E0 November 1984 INXYWB ROWNUM TYPE NSUBS VALUES RHS WEIGHT User s Guide SUBS T 32 34 SPK 35SPK 36SPK 37SPK 38SPK 39SPK 40SPK 41SPK 42SPK 43SPK 44SPK 45SPK 46SPK 47SPK 48SPK 49SPK 50SPK 51SPK 52SPK 53SPK 54SPK 55SPK 56SPK 57SPK 58SPK 59SPK 60SPK 6 1SPK 6 2SPK 63SPK 64SPK 65SPK 66 SPK 67 SPK 6 8 SPK 6 9SPK 70SPK 7 1 SPK 7 2SPK 73SPK 74 SPK 7 5SPK 76 SPK 77SPK 78SPK 79SPK 80SPK 81SPK 82SPK 83SPK 84SPK 85SPK 86SPK 87SPK 88SPK 89SPK 90SPK 91SPK 92SPK 93SPK 94 SPK 95SPK 96SPK 97SPK 98SPK 99SPK SPARSPAK B User s Guide CALL INXYWB ROWNUM TYPE NSUBS SUBS VALUES RHS WEIGHT T C Ci ee a m m m manm saw C ORDER COLUMNS C mc IE EDS a a di CALL ORCOLB T C C me M OE EA TUE nca C COMPUTE LEAST SQUARES SOLUTION AND RESIDUALS C ME EE rta MEC AN DU AR EED AE AS ED dole BEE EDS Pn a At CALL LSOSLV TOL TYPTOL T CALL RESIDB RESEQN RESCON T C C did iii EE iii C PRINT THE SOLUTION FOUND IN THE FIRST NCOLS C LOCATIONS IN THE WORKING STORAG
31. SPARSPAK Waterloo Sparse Matrix Package User s Guide for SPARSPAK B Alan George Esmond Ng Department of Computer Science University of Waterloo Waterloo Ontario CANADA Research Report CS 84 37 November 1984 Table of Contents 1 Introduction and basic structure of SPARSPAK B ee ee ee ee ee ee ee ee ee 2 Modules of SPARSPAK B and how to use them ee ee Be ee ee Re Ee de ee Re ee 2 1 User mainline program and an example sssessensenzznnnsanzenzennannnnsnnnpenazzamenanznnnznezzznzznnzz 2 2 Modules for input of the problem s ss ssenssnnazzensennnnennnnzanosenezzansenozzananenazanannenanzonazza Denise POWS ES es EE RN a i es 2 3 Module for ordering the columns sesse ee Ge Ge RA Ge Ee 2 4 Module for reordering the rows Bee ee ed ee Ee ee Ge ee ee once Mn 2 5 Module for numerical solution sssesssenennzennzznenznznnenzzznnananzanmronannanzznannantoszanenazzonana 3 Summary of the use of the basic interface modules PO EE ER ee EER ER 4 Save and restart facilities M AE EE O a A 5 Residual calculation iii 6 Output from SPARSPAK B nana na neo ttna bestens ona kann sana rc cono renos 6 1 Message level indicator MSGLVB anno nana nan Rd ee 6 2 Statistics gathering STATSB A EE A 6 3 Error messages IERRB e EA een eae 6 3 1 Save and restart subroutines
32. SPK 2SPK 3SPK 4SPK 5SPK 6SPK 7SPK 8SPK 9SPK 10SPK 11SPK 12SPK 13SPK 14SPK 15SPK 16SPK 17SPK 18SPK 19SPK 20SPK 21SPK 22SPK 23SPK 24SPK 25SPK 26SPK 27SPK 28SPK 29SPK 30SPK 31SPK 32SPK 33SPK 34SPK 35SPK 36SPK 37SPK 38SPK 39SPK 40SPK 41SPK 42SPK 43SPK 44SPK 45SPK 46SPK 47 SPK 48SPK 49SPK 50SPK SISPK 52SPK 53SPK 54SPK 55SPK 56SPK 57SPK 58SPK 59SPK 60SPK 61SPK SPARSPAK B c ai fai eres eed a ae ee eee C PRINT STATISTICS GATHERED BY SPARSPAK B UU 0 AA eeu a a Eme uere CALL STATSB STOP C END Output Xo k sko ko k k k k k kee KH RE sko sko sk o dk sko sko sko xxx kok sk sko Mk sko kk HE kok HM Sko sko sko sko KOK KH k kk k k k kk k k k k k sko ko sko kk kk kkk UNIVERSITY OF WATERLOO SPARSE MATRIX PACKAGE S P A R S P A K RELEASE 3 C JANUARY 1984 ANSI FORTRAN SINGLE PRECISION LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES OUTPUT UNIT FOR STATISTICS RSTRTB ORROWB LSOSLV RES I DB STATSB RESTART SYSTEM B FIND ROW ORDERING LEAST SOUARES SOLVE COMPUTE RESIDUAL SYSTEM B STATISTICS SIZE OF STORAGE ARRAY MAXSB NUMBER OF NUMBER OF NUMBER OF NUMBER OF NUMBER OF COLUMNS UNKNOWNS SPARSE EOUATIONS DENSE EQUATIONS SPARSE CONSTRAINTS DENSE CONSTRAINTS TIME FOR COLUMN ORDERING STORAGE FOR COLUMN ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION TIME FOR ROW ORDERING STORAGE FOR ROW ORDERING TIME FOR SOLUTI
33. SUBS 100 INTEGER FILE I ICASE IERRB IPRNTE IPRNTS 1 ISEED J 5 K KGRID KMI MAXINT 1 MAXSB MSGLVB NCOLS NDCONS NDEONS NSCONS 1 NSEONS NSUBS OUTPUT ROWNUM TYPE TYPTOL REAL MCHEPS RATIOL RATIOS TIME REAL T 7500 VALUES 100 REAL RESCON RESEON RHS TOL WEIGHT C Qo oko ok ok ook ok ook ok ok ke k kok Ok ok ok ok ok ok ok ok ok ok ke kk k k k ok ok ak k k OK oko Gk OE OE Gk k ek ok ok k k KOK sko pk November 1984 24 I SPK 2SPK 3SPK 4 SPK 3SPK 6 SPK 7 SPK 8 SPK 9SPK 10SPK 11SPK 12SPK 13SPK 14SPK 15SPK I6SPK I7SPK I8SPK 19SPK 20SPK 21SPK 22SPK 23SPK 24SPK 25SPK 26SPK SPARSPAK B 1 1 COMMON SPKSYS IPRNTE IPRNTS MAXINT RATIOS MCHEPS TIME COMMON SPBUSR MSGLVB IERRB NDEQNS NSCONS NDCONS User s Guide RATIOL NSEQNS Q o ox ok oko ok RR skoke b X RE e kkk kakka kk kk kk kk kk RR ak XE ME G0 0 kok kok G0 0x ER G0 OE aa 0X EE oko sk Xx KEE anaa O AAG a QOO O OOGOO CLC CALL SPRSPK FILE I OUTPUT IPRNTS TOL MCHEPS TOL 100 0E0 TOL TYPTOL I 100 1 200 300 ISEED 1234567 MSGLVB 2 MAXSB 7500 CALL FILEB FILE e w A e N EE MO M mom NSUBS 4 WEIGHT TYPE ROWNUM KGRID l 5 KMI KGRID 1 DO 400 DO SUBS 1 I 1 KGRID J SUBS 2 I 1 KGRID J 1 SUBS 3 I KGRID J SUBS 4 I KGRID J 1 GENERATE
34. TRTB SAVE T OPTION 4 GO TO 500 C END Output seen eee eee UNIVERSITY OF WATERLOO rrrxxx SPARSE MATRIX PACKAGE k k sko sko sko ik 0x kok kk o OE KKK sko ko sk k dk kk sk k x kkk SPARS PAK RELEASE 3 C JANUARY 1984 HEeREKEEKE ANS FORTRAN wee SINGLE PRECISION tx x LAST UPDATE JANUARY 1984 OUTPUT UNIT FOR ERROR MESSAGES OUTPUT UNIT FOR STATISTICS FILEB INXYWB ORCOLB SAVEB ORROWB LSOSLV RES I DB STATSB FILE INITIALIZATION INPUT ROWS FIND COLUMN ORDERING SAVE STORAGE VECTOR FIND ROW ORDERING LEAST SOUARES SOLVE COMPUTE RESIDUAL SYSTEM B STATISTICS November 1984 SIZE OF STORAGE ARRAY MAXSB NUMBER OF COLUMNS UNKNOWNS NUMBER OF SPARSE EQUATIONS NUMBER OF DENSE EQUATIONS NUMBER OF SPARSE CONSTRAINTS NUMBER OF DENSE CONSTRAINTS TIME FOR COLUMN ORDERING STORAGE FOR COLUMN ORDERING TIME FOR ALLOCATION STORAGE FOR ALLOCATION TIME FOR ROW ORDERING STORAGE FOR ROW ORDERING TIME FOR SOLUTION OPERATION COUNT FOR SOLUTION STORAGE FOR SOLUTION il l 7000 100 324 0 2269 0 1966 0 0 4 379405 2306 User s Guide 217 050 683 37 110SPK 111SPK 112SPK 113SPK 114SPK 115SPK 116SPK 117SPK 118SPK 119SPK 120SPK 121SPK 122SPK SPARSPAK B TIME FOR COMPUTING RESIDUAL OPN COUNT FOR COMPUTING RESIDUAL STORE FOR COMPUTING RESIDUAL TOTAL TIME REQUIRED MAXIMUM STORAGE REQUIRED RESIDUAL IN EQUATIONS RESIDUAL IN CONSTRAINTS RSTRTB OR
35. comparing SPARSPAK B with other packages or is going to solve numerous similar problems and wants to adjust the working storage to the minimum necessary The package has a common block called SPBDTA containing variables whose values can be printed by executing the following FORTRAN statement CALL STATSB The information printed includes the size of the working storage array 7 the number of columns in X and Y the number of rows in A B E and F the number of nonzeros in the matrices X and Y the maximum number of nonzeros in the rows of X and Y the number of nonzeros in E T AT A the number of off diagonal nonzeros in the upper triangular matrix obtained after the numerical reduction the time used to find the column ordering the storage used by the column ordering subroutine the time used for data structure set up the storage used by the storage allocation subroutine the time used to find the row ordering the storage used by the row ordering subroutine the time used for computing the solution the number of operations reguired by the solution subroutine November 1984 17 SPARSPAK B User s Guide the storage used by the solution subroutine the time used for computing the residuals the number of operations reguired by the residual calculation subroutine the storage used by the residual calculation subroutine the total time used by the solution process the maximum storage reguired by
36. e cause of error routine ORCOLB was not executed successfully Incorrect execution sequence Probable cause of error routine ORROWB was called after having already been executed successfully Input row ordering option OPTION is invalid Insufficient storage in working storage array to execute the row ordering routine The minimum value of MAXSB required is printed in the error message Response execute SAVEB and restart the computation using ORROWB with MAXSB at least as large as that indicated in the error message 6 3 6 Solution subroutine November 1984 SPARSPAK B IERRB 251 252 253 254 255 256 User s Guide LSOSLV Incorrect execution seguence Probable cause of error routine ORCOLB or routine ORROWB was not executed successfully Incorrect execution seguence Probable cause of error routine LSOSLV was called after having already been executed successfully Input tolerance TOL is negative Input tolerance type TYPTOL is invalid Insufficient storage in working storage array to execute the solution routines The minimum value of MAXSB reguired is printed in the error message Response execute SAVEB and restart the computation using LSOSLV with MAXSB at least as large as that indicated in the error message routine LSOSLV fails to compute a singular value decomposition of intermediate small dense matrices 6 3 7 Residual calculation subroutine IERRB 261 262 November 1984
37. e storage of intermediate results Before beginning to solve a problem the user must call a subroutine FILEB to tell SPARSPAK B which FORTRAN unit it should use The FORTRAN statement to be used is CALL FILEB IONUM where ZONUM is the reguired FORTRAN logical unit Important Notes 1 The module FILEB must be called when a new problem is to be solved 2 The user is responsible for defining the external file for the FORTRAN logical unit JONUM using the appropriate system control statement or command This depends on the environment in which the program is being executed Furthermore the file should be preserved throughout the execution of the program Warning The modules of SPARSPAK B communicate with each other through labelled common blocks whose names are SPKSYS SPBUSR SPBCON SPBMAP and SPBDTA Note that SPKSYS is shared by SPARSPAK A and SPARSPAK B Thus the user must not use labelled common blocks with these names in his program 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 specifically reguested names These names should be specified when the package is acguired 2 2 Modules for input of the problem The subroutine ZVXYWB allows the user to provide the structure and the numerical values of X and Y to the package along with the numerical values of y z Wx and Wy They are provided one
38. ll ORCOLB However the interface will detect that the ordering and or storage allocation have already been performed and will skip that part of the computation 2 See Section 4 for details on how to use SAVEB and RSTRTB and Examples 3 and 4 in Section 8 November 1984 10 SPARSPAK B User s Guide 2 4 Module for reordering the rows The execution time or the numerical stability of the module LSOSLV described in Section 2 5 can be affected significantly by the row ordering 3 13 Accordingly SPARSPAK B provides a row ordering module which may be invoked optionally by the following FORTRAN statement CALL ORROWB OPTION T When the weights W and Wy vary widely in magnitude it is important for numerical accuracy that the rows of 4 and E be arranged in order of increasing weight See 13 for details This can be achieved by setting the integer parameter OPTION to 1 If the rows of WyX and WyY do not vary greatly the user may wish to sort the rows in order to reduce execution time although the user should be aware that this sorting might require more storage than would otherwise be required to solve the problem The following options are provided OPTION Details 0 Rows are processed in the order they were supplied regardless of the values of ROWNUM see Section 2 2 l Rows are sorted in order of increasing weight i e the parameter WEIGHT 2 Rows are sorted in order of increasing number of nonzeros i e the para
39. meter NSUBS 3 Rows are sorted in order of increasing minimum column subscripts see below 4 Rows are sorted in order of increasing maximum column subscripts see below 5 Rows are sorted in order of the parameter ROWNUM see Section 2 2 Here the maximum minimum column subscript of a row is the column subscript index of the last first nonzero in that row Moreover the column indices referred to those of the column permuted matrix obtained from ORCOLB The effectiveness of these strategies is not well understood and varies with the problem In the absence of any prior knowledge setting OPTION to 4 is recommended See 8 9 10 for discussions on the row ordering problem in the solution of sparse linear least sguares problems Note The amount of storage reguired to perform row ordering may be substantially larger than those reguired by the other interface subroutines 2 5 Module for numerical solution The actual numerical computation which produces the solution B is initiated by executing the FORTRAN statement CALL LSOSLV TOL TYPTOL T where T is the working storage array November 1984 11 SPARSPAK B User s Guide The parameter TOL is a user specified tolerance which is used to determine when a diagonal element of the upper triangular matrix produced in the numerical computation should be regarded as numerically zero Suppose R denotes the upper triangular matrix A diagonal element R will be regarded a
40. mj me 3 am mm a mm en a m e wm m m CALL SPRSPK MSGLVB 2 MAXSB 5000 CALL FILEB FILE anan N REWIND INPUT 100 CONTINUE READ INPUT ROWNUM TYPE NSUBS 1 SUBS K VALUES K K 1 NSUBS 1 Declared either REAL or DOUBLE PRECISION depending on the version of SPARSPAK A and SPARSPAK B that is available The examples in this manual assume a single precision version is being used November 1984 6 SPARSPAK B User s Guide 1 RHS WEIGHT IF NSUBS EQ 0 GO TO 200 CALL INXYWB ROWNUM TYPE NSUBS SUBS VALUES 1 RHS WEIGHT T GO TO 100 200 CONTINUE Cs fate o TES C ORDER COLUMNS AND ROWS CU usu E a E CALL ORCOLB T CALL ORROWB OPTION T we wa aw a ou a am aw m an m am aw w w a m vw sa am EE m em m m am m m e a m an m a e sm e a ou m am en m au ww Mm a m a m am w am mw m m e vw w am en a w an a am am en m m au m m CALL LSOSLV TOL TYPTOL T CALL RESIDB RESEON RESCON T OON S Cc e tu tu A ta ta S A A tay ta ta O te G aj head O 2 gt gt A ty ta lan 3 A 4 ta PII PRINT THE SOLUTION FOUND IN THE FIRST NCOLS LOCATIONS IN THE WORKING STORAGE ARRAY T AND PRINT THE RESIDUALS WRITE 6 11 T K K lt l NCOLS 11 FORMAT 10H SOLUTION 5F12 5 WRITE 6 22 RESEON RESCON 22 FORMAT 22H EQUATION RESIDUAL FI2 5 d 22H CONSTRAINT RESIDUAL rA cssc S 0 OO
41. n be used to solve a sparse general square system of linear equations p 0 and m n As we have noted above the package does not impose any restrictions on the dimensions of X and Y That is the package also handles the case in which m 0 but p 0 even though such a problem may not have any physical meaning In this situation SPARSPAK B simply treats the nonempty constraints as a linear least sguares problem and computes the minimal norm least sguares solution The basic computational technigue used in SPARSPAK B is due to George and Heath 3 Y Let M denote the m t p Xn matrix d It is assumed that the Cholesky factor of MT M is R sparse Then M is reduced to upper trapezoidal form R by applying Givens transformations and Gaussian eliminations to the rows of M Here R is an nXn upper triangular matrix Finally the upper triangular matrix R is used to compute the solution It can be shown that the structure of R is contained in the structure of the Cholesky factor of the symmetric positive definite matrix MTM Thus using techniques developed for solving sparse symmetric positive definite systems one can predict the structure of R by analyzing the structure of M M 7 This allows a storage scheme for R to be set up before numerical factorization begins and hence the numerical computation of R can be carried out using a static data structure Experience has shown that this approach is efficient both in terms of storage and executio
42. n time compared to schemes that employ dynamic storage allocation 4 In most cases the matrix MTM is sparse when both the least squares matrix X and the constraint matrix Y are sparse However there are instances in which M M is dense even though X and Y are sparse Fortunately this usually occurs when there are only a few dense rows in X and Y Such problems can be handled in a special way by the package More precisely the least squares problem and the constraints may be regarded as being partitioned as follows November 1984 3 3 SPARSPAK B User s Guide xy Hd mee z 2 where A and E contain respectively the sparse equations and constraints and B and F contain respectively the dense equations and constraints It is now assumed that M denotes the matrix E A Furthermore it is assumed that the Cholesky factor of M M is sparse The package reduces M to upper trapezoidal form using the approach described above Then the upper trapezoidal form together with B and F are used to derive the solution to the original problem The algorithm is due to Bjorck 1 Detailed description of the implementation can be found in 11 In general the algorithms the data structures and the storage management for solving sparse matrix problems are quite complicated Thus in order to insulate the user from these considerations a set of simple user interface subroutines is used in the package These interface
43. new interface subroutines and some additional underlying subroutines but makes extensive use of subroutines in the SPARSPAK A package SPARSPAK B cannot run without SPARSPAK A Although SPARSPAK B depends upon SPARSPAK A it is in an important sense independent a user can simultaneously solve a least squares problem and a problem from the class treated by the basic SPARSPAK A package t Department of Computer Science University of Waterloo Waterloo Ontario CANADA November 1984 SPARSPAK B User s Guide IMPORTANT NOTE The numerical algorithm used in LSQSLV see Section 2 5 is a prototype implementation of an algorithm that is due to Bjorck 1 The behaviour of the numerical algorithm in the presence of roundoff errors in finite precision computer arithmetic is not well understood In the course of solving a constrained linear least squares problem some small dense subproblems may have to be solved Some of these resulting subproblems may be sensitive to roundoff errors and numerical solutions to these subproblems may be inaccurate when a small tolerance is supplied to LSQSLV We would appreciate receiving any comments and feedback the user may have in using the package to solve practical problems Such comments and feedback are important and useful since they may allow us to refine the numerical algorithm in the future When the package fails to produce an accurate solution we would be grateful if the user could send us a copy
44. ng point array which for purposes of discussion we will denote by T In addition the user must provide its size MAXSB which is transmitted to the package via a common block SPBUSR SPARSPAK B USER which has eight variables COMMON SPBUSR MSGLVB IERRB MAXSB NCOLS NSEONS NDEONS NSCONS NDCONS Here MSGLVB is the message level indicator which is used to control the amount of information printed by the package The second variable ZERRB 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 MSGLVB and IERRB is provided in Section 6 The variable NCOLS is the number of columns in X and Y and NSEONS NDEONS NSCONS NDCONS are respectively the number of rows in A B E and F Thus NSEONS NDEONS m and NSCONS NDCONS p The following program illustrates how one might use SPARSPAK B The various subroutines referenced are described in the subseguent parts of this section The problem solved is assumed to be stored on an external file FORTRAN unit 1 in a binary format INTEGER SUBS 10 INTEGER FILE IERRB INPUT K MAXSB MSGLVB 1 NCOLS NDCONS NDEONS NSCONS NSEONS 1 NSUBS OPTION ROWNUM TYPE TYPTOL REAL T 5000 VALUES 10 REAL RESCON RESEON RHS TOL WEIGHT COMMON SPBUSR MSGLVB IERRB MAXSB NCOLS 1 NSEONS NDEONS NSCONS NDCONS a m am a aw a m am mm 00m mu Z m man A e bm N M
45. ns corresponding to the column indices stored in SUBS RHS The right hand side element of y or z corresponding to the input row WEIGHT The diagonal element of the weight matrix Wy or Wy corresponding to the input row It is assumed that WEIGHT is positive T The floating point working storage array from which all storage for the package is allocated See Section 2 1 and the example there Dense rows Some problems may yield a few rows of X and of Y that have relatively many nonzeros Such rows wreak havoc with the sparsity preservation techniques used in SPARSPAK B At the moment the package has no robust scheme for deciding which rows will cause unacceptable damage but it does have a way of circumventing problems caused by such rows Accordingly the user can indicate dense rows of X or of Y by setting the corresponding value of the input parameter TYPE to 2 or 4 when such a row is input See Example 2 in Section 8 for an illustration of the use of this feature Important Note The floating point values transmitted to SPARSPAK B by INXYWB are either single or double precision floating point numbers depending on the version of SPARSPAK B being used The examples in this manual assume that a single precision version of the package is being used 2 3 Module for ordering the columns Recall that the Cholesky factor of the symmetric positive definite matrix M M plays an E important role in the solution process where M
46. of the data if possible so that we may locate where problems occur Please send comments and feedback to Dr Alan George or Dr Esmond Ng Department of Computer Science University of Waterloo Waterloo Ontario CANADA N2L 3Gl Telephone 519 885 1211 ext 3473 Dr Alan George 519 885 1211 ext 2517 Dr Esmond Ng November 1984 ho SPARSPAK B User s Guide 1 Introduction and basic structure of SPARSPAK B SPARSPAK B is designed to solve the problem i X 9 Ben Wx B y IE where O B B minimize Wy YB z aj Here X and Y are respectively m Xn and p Xn sparse matrices y and z are vectors of length m and p respectively and B is a vector of length n The matrices Wy and Wy are respectively m Xm and pXp diagonal weight matrices The package is capable of handling the very general problem in which X and Y may be either of full rank or rank deficient and there is no restriction on the dimensions of X and Y Moreover the constraints need not be consistent When the constrained linear least sguares problem does not have a unigue solution the package computes the solution which has minimal Euclidean norm Even though SPARSPAK B is designed to solve a very general constrained linear least sguares problem it can also handle problems that may be much simpler than this general case For example the package is capable of solving an unconstrained linear least sguares problem that is p 0 In particular SPARSPAK B ca
47. ooo ER kk RR EE ER OR kok EE ER kk dk C C C INITIALIZATION C O CALL SPRSPK C FILE 1 SAVE 2 OPTION 4 TOL MCHEPS TOL 100 0E0 TOL TYPTOL I ISEED 1234567 November 1984 39 1 SPK 2SPK 3SPK 4 SPK 5SPK 6 SPK 7 SPK 8 SPK 9SPK 10SPK 11SPK 12SPK 13SPK 14S PK 15SPK 16S PK 17SPK 18SPK 19S PK 20SPK 21SPK 22S PK 23S PK 24SPK 25S PK 26SPK 27SPK 28S PK 29SPK 30SPK 31SPK 32S PK 33S PK 34S PK 35SPK 36SPK 37SPK 38S PK 39S PK 40SPK 41SPK 42SPK 43SPK 44SPK 45SPK 46SPK 47SPK 48SPK SPARSPAK B User s Guide C MSGLVB 2 MAXSB 3000 C l CALL FILEB FILE C Coon nn nr rn m m m o m m m C GENERATE PROBLEM FROM THE GRID CG riada A SS l EN IN RR NSUBS 4 WEIGHT 1 0E0 TYPE 1 ROWNUM 0 C KGRID 10 KM1 KGRID 1 DO 400 I 1 KMI DO 300 J 1 KMI O ll o C GENERATE STRUCTURE C o EG EE kwi m ma anan SUBS 1 I 1 KGRID J SUBS 2 I 1 KGRID J I SUBS 3 I KGRID J SUBS 4 I KGRID J 1 C DS a ak RI N a NA Ned at tan PIT C GENERATE NUMERICAL VALUES USING C A RANDOM NUMBER GENERATOR AAA RADE E E OE OE A aa DO 200 ICASE 1 4 DO 100 K 1 NSUBS VALUES K RANDOM ISEED 100 CONTINUE ROWNUM ROWNUM 1 RHS RANDOM ISEED CALL INXYWB ROWNUM TYPE NSUBS SUBS 1 VALUES RHS WEIGHT T 200 CONTINUE 300 CONTINUE 400 CONTINUE IF IERRB NE 0 GO TO 500 EYE AE ED en a 4 a om e a a
48. problem initialization module FILEB 22 row input module INXYWB 23 column ordering and data structure set up module ORCOLB 24 row reordering module ORROWB 25 solution module LSQSLV 26 residual calculation module RESIDB November 1984 18 SPARSPAK B User s Guide 6 3 1 Save and restart subroutines IERRB SAVEB and RSTRTB 201 Output unit given to SAVEB is not positive 202 Input unit given to RSTRTB is not positive 203 Insufficient storage in working storage array to restart the computational process The minimum value of MAXSB required is printed in the error message 6 3 2 Problem initialization subroutine IERRB FILEB 211 Input output unit given to FILEB is not positive 6 3 3 Row input subroutine JERRB INXYWB 221 Incorrect execution sequence Probable cause of error routine FILEB was not executed successfully 222 Incorrect execution seguence Probable cause of error routine ORCOLB has already been executed To start a new problem FILEB must be called first 223 Number of nonzeros NSUBS is not positive 224 Input row type TYPE is invalid 225 Input index or subscript is not positive 226 Input weight WEIGHT is not positive 22 Insufficient storage in working storage array to form matrix structure The minimum value of MAXSB reguired is printed in the error message 6 3 4 Column ordering subroutine IERRB ORCOLB 231 Incorrect execution sequence Probable cause of error routine I
49. r handling sparse symmetric positive definite systems of linear equations Indeed SPARSPAK B makes extensive use of the subroutines in SPARSPAK A which is a package designed for solving efficiently sparse symmetric positive definite linear systems 2 SPARSPAK B cannot run November 1984 4 SPARSPAK B User s Guide without SPARSPAK A Although SPARSPAK B depends upon SPARSPAK A it is in an important sense independent A user can simultaneously solve a sparse constrained linear least squares problem and a sparse symmetric positive definite linear system in the same program An early version of SPARSPAK B was designed and implemented by Dr M T Heath at the Oak Ridge National Laboratory This early version includes algorithms for solving the basic sparse unconstrained linear least squares problems and for handling constraints dense rows and rank deficiency 12 Those algorithms are special cases of a more general algorithm due to Bjorck for handling more general sparse constrained linear least sguares problems 1 This general algorithm is used in the current version of SPARSPAK B The design of SPARSPAK B is similar to that of SPARSPAK A 2 The reader is referred to 6 for a discussion of the design and implementation issues November 1984 5 SPARSPAK B User s Guide 2 Modules of SPARSPAK B and how to use them 2 1 User mainline program and an example SPARSPAK B allocates all of its storage from a single one dimensional floati
50. riate 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 RSTRTB subroutine c The external file for the logical unit K must be different from the working file specified in the FILEB statement November 1984 l 14 SPARSPAK B User s Guide 5 Residual calculation SPARSPAK B provides a subroutine for computing the Euclidean norm of the residual vectors The FORTRAN statement to be used is as follows CALL RESIDB RESEON RESCON T where T is the working storage array After the subroutine is called RESEON ll y XB Il and RESCON z YB ll I where B is the computed solution Important Note RESIDB should be called only when the least squares solution has been computed that is after LSQSLV has been executed successfully November 1984 15 SPARSPAK B User s Guide 6 Output from SPARSPAK B As noted earlier in Section 2 1 the user supplies a one dimensional floating point array 7 from which all array storage is allocated In particular the interface allocates the first NCOLS storage locations in T for the solution vector of the constrained linear least sguares problem After all the interface modules for a particular problem have been successfully executed the user can retrieve the solution from these NCOLS locations In addition to the least sguares solution
51. s numerically zero if R lt TOL when TYPTOL 0 and Ri AMA TOL when TYPTOL lt 1 The choice of TOL is usually problem dependent In general TOL should be chosen so that it reflects the accuracy of X Y y and z For example if the input numerical values are known to be accurate to t significant digits then one should use a relative test TYPTOL lt 1 and TOL should be set to about 107 November 1984 12 SPARSPAK B User s Guide 3 Summary of the use of the basic interface modules The following flowchart depicts the seguence of calls to SPARSPAK B interface subroutines which occur when a constrained linear least sguares problem is solved After LSOSLV has been executed the solution to the problem can be found in the first NCOLS locations in the working storage array T where NCOLS is the third variable in the labelled common block SPBUSR SPRSPK A second and subseguent problems can be solved by simply starting at the beginning of the module seguence again as implied by the broken line l November 1984 13 SPARSPAK B User s Guide 4 Save and restart facilities As in SPARSPAK A SPARSPAK B provides two subroutines called SAVEB and RSTRTB which allow the user to stop the calculation at some point save the results on an external seguential file and then restart the calculation at exactly the same point some time later To save the results of the computation done thus far the user executes the FORTRAN sta
52. subroutines will be described in detail in the following sections The user and the package interact to solve the problem through the following basic steps Step 1 The user supplies the rows of X y and Y z to the package in any order along with the corresponding diagonal element of Wy and Wy Step 2 E The user calls a subroutine which initiates a reordering of the columns of 4 in order to preserve sparsity in subseguent calculations See Section 2 3 Step 3 The user calls a subroutine which initiates a reordering of the rows of X and Y according to one of several criteria See Section 2 4 Step 4 The user calls a subroutine which computes the solution B The package has facilities to allow the user to compute conveniently the norm of the residual vectors I XB Il and llz rgll where B is the computed solution Important notes l Note that SPARSPAK B is designed to handle problems in which X and Y may have any dimensions However the package will perform more efficiently in terms of storage and execution time when m p2n that is overdetermined problems When m t p n that is underdetermined problems the package will still be able to solve the problems but both the storage and time requirements may be high See 5 for some examples and algorithms for handling sparse unconstrained underdetermined problems 2 The discussions above indicate that the package employs heavily techniques fo
53. t where X is defined below Consider a k Xk grid see figure below where k 4 There is a variable associated with each grid point and for each sguare in the grid there is a set of four eguations involving the variables at the four grid points in that sguare This gives rise to a sparse overdetermined system of linear eguations and it will be solved in the least sguares sense The number of variables is n k and the number of equations is m 4 k 1 For our purpose the numerical values of the nonzeros in X and y are generated using a uniform random number generator November 1984 23 SPARSPAK B User s Guide Example 1 In this example SPARSPAK B is used to solve a constrained linear least sguares problem X A min XB y Ila where O B B minimizes YB z ll The matrix X is the same as the one defined on a k Xk grid The constraints are as follows 100 0 2 010 0 AES RS That is the first and the second variables are assumed to have values 2 and 5 respectively The program begins by calling SPRSPK to initialize SPARSPAK B followed by a call to FILEB which tells SPARSPAK B the unit number to be used to reference the temporary external file required by the package Then the problem is generated and provided to SPARSPAK B using the subroutine INXYWB The column ordering is determined and data structure is set up by calling ORCOLB Finally the solution is obtained by executing LSOSLV The subroutines RESI
54. tement CALL SAVEB K T 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 FORTRAN statement CALL RSTRTB K T Example 3 provided in Section 8 illustrates the use of SAVEB and RSTRTB Note that executing SAVEB does not destroy any information the computation can proceed just as if SAVEB were not executed When errors occur in a module the routines SAVEB and RSTRTB are useful in saving the results of previous modules executed successfully see Section 6 3 and Example 4 in Section 8 Another potential use of the SAVEB and RSTRTB modules is to make the working storage array T available to the user in the middle of a sparse matrix computation After SAVEB has been executed the working storage array T can be used by some other computation Finally the SAVEB and RSTRTB 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 SAVEB and RSTRTB information is either written on or read from the FORTRAN logical unit K using binary format b If the subroutines SAVEB and RSTRTB are used then before the user executes his program he must define a file for the FORTRAN logical unit K using the approp
55. utation where even summary information would generate excessive output Increasing the value of MSGLVB up to 4 provides increasingly detailed information about the computation Note that the module SPRSPK sets MSGLVB to its default value if the user wishes MSGLVB to be different from two he must reset it after SPRSPK has been called In many circumstances SPARSPAK B will be embedded in still another super package which models phenomena producing sparse constrained linear least sguares problems Messages printed by SPARSPAK B 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 corrective action which makes the error messages irrelevant Thus a printing by SPARSPAK B can be prevented by setting MSGLVB to zero We summarize our discussion in this section in the following table November 1984 16 SPARSPAK B User s Guide MSGLVB amount of output 0 no information is provided l 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 a high volume of output may be generated if MSGLVB is set to four since the input data would also be echoed 6 2 Statistics gathering STATSB SPARSPAK B gathers a number of statistics which the user will find useful if he is
56. w e a am am am m m m am wm an an CALL ORCOLB T IF IERRB NE 0 GO TO 500 CALL ORROWB OPTION T IF IERRB NE 0 GO TO 500 m m m m am am m m m m a m m a a m m m m am a m ew am G BOY om em a am w a SE m am m m m a w m m a w am m e m a m eu a am m m m m MM w m m ME m a n CALL LSOSLV TOL TYPTOL T IF LERRB NE 0 GO TO 500 CALL RESIDB RESEQN RESCON T IF IERRB NE 0 GO TO 500 a mm ma EE ME ME EE EE SE EE EE E SR M 021 1 1 1 1 SE EE ME EE anaoa CALL STATSB STOP a 500 CONTINUE November 1984 40 49SPK 50SPK 51SPK 52SPK 53SPK 54SPK 55SPK 56SPK 57SPK 58SPK 59SPK 60SPK 6ISPK 62SPK 63SPK 64SPK 65SPK 66SPK 67SPK 68SPK 69SPK 70SPK 71SPK 72SPK 73SPK 74SPK 75SPK 76SPK 77SPK 78SPK 79SPK 80SPK 81SPK 82SPK 83SPK 84SPK 85SPK 86SPK 87SPK 88SPK 89SPK 90SPK 91SPK 92SPK 93SPK 94 SPK 95SPK 96 SPK 97SPK 98SPK 99SPK 100SPK 101SPK 102SPK 103SPK 104SPK 105SPK 106SPK 107SPK 108SPK 109SPK 110SPK 111SPK 112SPK 113SPK 114SPK SPARSPAK B CALL SAVEB SAVE T CALL STATSB STOP END Output kok kok KKK he k k k kkk KK ok kok kok k hak koko kok k k ak k sko kkk kk dk k k k kk ok kk k ee k kkk KK k k kk kk kkk ho kok k k k k k k ok k k UNIVERSITY OF WATERLOO SPARSE MATRIX PACKAGE S PAR S PAK RELEASE 3 C JANUARY 1984 ANSI FORTRAN SINGLE PRECISION LAST UPDATE
57. ze of the working storage array we then execute the second program The routine RSTRTB is called to restore the state of the computation before ORROWB is called Program 1 C SPARSPAK B ANSI FORTRAN RELEASE III NAME EX4A C C UNIVERSITY OF WATERLOO JANUARY 1984 CHO KOK OE OE EE OE ak k ok OE OKO EG EG EE EE k k GE OE k k eee X AR CEPR EEE ER RE ER RR OE ER RR RE ER RE ER EER RR EE RR EE RR EE OR EER A EE GGG GIG EG Cee RRR RR EEE MAINLINE PROGRAM okck E E GE EE E QN kok AE ak kk k k ak ak EE kkk k k k k kok k kkk k kokok k k OO RE RE EG GG AGE GE RR RE RE ER Ek Ed N EG k CO OE GE OE OE OE OR ER OR OE OE OE OE OE EE k RR RE ER RR O C C EXAMPLE 4 SEE USER S GUIDE C C REOUIREMENTS C AN EXTERNAL FILE FOR UNIT 1 C AN EXTERNAL FILE FOR UNIT 2 C A RANDOM NUMBER GENERATOR RANDOM EA C INTEGER SUBS 100 INTEGER FILE 9d 4CASE IERRB IPRNTE IPRNTS 1 ISEED J gt K KGRID KMI MAXINT 1 MAXSB MSGLVB NCOLS NDCONS NDEONS NSCONS 1 NSEONS NSUBS OPTION ROWNUM SAVE TYPE 1 TYPTOL REAL MCHEPS RATIOL RATIOS TIME REAL T 3000 VALUES 100 REAL RESCON RESEON RHS TOL WEIGHT C CA Ra ER Gk e RE Ek hn dk dk EE OR o Gk RR ok Ok e Ge SR Go ke ke ae ade dk e e OE Gb skok RR EE EE kkk kk ss C COMMON SPKSYS IPRNTE IPRNTS MAXINT RATIOS RATIOL 1 MCHEPS TIME COMMON SPBUSR MSGLVB IERRB MAXSB NCOLS NSEONS 1 NDEONS NSCONS NDCONS C C RE EE ER ER EE RE EE EE Ek EK RE Ek EE
Download Pdf Manuals
Related Search