C APPENDIX 2 1/15/81 C C SUBROUTINES FOR SOLVING SPARSE NONSYMMETRIC SYSTEMS C OF LINEAR EQUATIONS (TRACK NONZEROES DYNAMICALLY) C C C*** SUBROUTINE TDRV C*** DRIVER FOR SUBROUTINE FOR SOLVING SPARSE NONSYMMETRIC SYSTEMS OF C LINEAR EQUATIONS (TRACK NONZEROES DYNAMICALLY) C SUBROUTINE TDRV * (N, R,IC, IA,JA,A, B, Z, NSP,ISP,RSP,ESP, FLAG) C C PARAMETERS C CLASS ABBREVIATIONS ARE -- C N - INTEGER VARIABLE C F - REAL VARIABLE C V - SUPPLIES A VALUE TO THE DRIVER C R - RETURNS A RESULT FROM THE DRIVER C I - USED INTERNALLY BY THE DRIVER C A - ARRAY C C CLASS \ PARAMETER C ------+---------- C \ C THE NONZERO ENTRIES OF THE COEFFICIENT MATRIX M ARE STORED C ROW-BY-ROW IN THE ARRAY A. TO IDENTIFY THE INDIVIDUAL NONZERO C ENTRIES IN EACH ROW, WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY C LIES. THE COLUMN INDICES WHICH CORRESPOND TO THE NONZERO ENTRIES C OF M ARE STORED IN THE ARRAY JA; I.E., IF A(K) = M(I,J), THEN C JA(K) = J. IN ADDITION, WE NEED TO KNOW WHERE EACH ROW STARTS AND C HOW LONG IT IS. THE INDEX POSITIONS IN JA AND A WHERE THE ROWS OF C M BEGIN ARE STORED IN THE ARRAY IA; I.E., IF M(I,J) IS THE FIRST C NONZERO ENTRY (STORED) IN THE I-TH ROW AND A(K) = M(I,J), THEN C IA(I) = K. MOREOVER, THE INDEX IN JA AND A OF THE FIRST LOCATION C FOLLOWING THE LAST ELEMENT IN THE LAST ROW IS STORED IN IA(N+1). C THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS GIVEN BY C IA(I+1) - IA(I), THE NONZERO ENTRIES OF THE I-TH ROW ARE STORED C CONSECUTIVELY IN C A(IA(I)), A(IA(I)+1), ..., A(IA(I+1)-1), C AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN C JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1). C FOR EXAMPLE, THE 5 BY 5 MATRIX C ( 1. 0. 2. 0. 0.) C ( 0. 3. 0. 0. 0.) C M = ( 0. 4. 5. 6. 0.) C ( 0. 0. 0. 7. 0.) C ( 0. 0. 0. 8. 9.) C WOULD BE STORED AS C \ 1 2 3 4 5 6 7 8 9 C ---+-------------------------- C IA \ 1 3 4 7 8 10 C JA \ 1 3 2 2 3 4 4 4 5 C A \ 1. 2. 3. 4. 5. 6. 7. 8. 9. . C C NV \ N - NUMBER OF VARIABLES/EQUATIONS. C FVA \ A - NONZERO ENTRIES OF THE COEFFICIENT MATRIX M, STORED C \ BY ROWS. C \ SIZE = NUMBER OF NONZERO ENTRIES IN M. C NVA \ IA - POINTERS TO DELIMIT THE ROWS IN A. C \ SIZE = N+1. C NVA \ JA - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF A. C \ SIZE = SIZE OF A. C FVA \ B - RIGHT-HAND SIDE B; B AND Z CAN THE SAME ARRAY. C \ SIZE = N. C FRA \ Z - SOLUTION X; B AND Z CAN BE THE SAME ARRAY. C \ SIZE = N. C C THE ROWS AND COLUMNS OF THE ORIGINAL MATRIX M CAN BE C REORDERED (E.G., TO REDUCE FILLIN OR ENSURE NUMERICAL STABILITY) C BEFORE CALLING THE DRIVER. IF NO REORDERING IS DONE, THEN SET C R(I) = C(I) = IC(I) = I FOR I=1,...,N. THE SOLUTION Z IS RETURNED C IN THE ORIGINAL ORDER. C C NVA \ R - ORDERING OF THE ROWS OF M. C \ SIZE = N. C NVA \ IC - INVERSE OF THE ORDERING OF THE COLUMNS OF M; I.E., C \ IC(C(I)) = I FOR I=1,...,N, WHERE C IS THE C \ ORDERING OF THE COLUMNS OF M. C \ SIZE = N. C C VARIOUS ERRORS ARE DETECTED BY THE DRIVER AND THE INDIVIDUAL C SUBROUTINES. C C NR \ FLAG - ERROR FLAG; VALUES AND THEIR MEANINGS ARE -- C \ 0 NO ERRORS DETECTED C \ N+K NULL ROW IN A -- ROW = K C \ 2N+K DUPLICATE ENTRY IN A -- ROW = K C \ 5N+K NULL PIVOT -- ROW = K C \ 8N+K ZERO PIVOT -- ROW = K C \ 10N+1 INSUFFICIENT STORAGE IN TDRV C \ 12N+K INSUFFICIENT STORAGE IN TRK C C WORKING STORAGE IS NEEDED FOR THE FACTORED FORM OF THE MATRIX C M PLUS VARIOUS TEMPORARY VECTORS. THE ARRAYS ISP AND RSP SHOULD BE C EQUIVALENCED; INTEGER STORAGE IS ALLOCATED FROM THE BEGINNING OF C ISP AND REAL STORAGE FROM THE END OF RSP. C C NV \ NSP - DECLARED DIMENSION OF RSP; NSP GENERALLY MUST BE C \ LARGER THAN 6N+2 + 2*K (WHERE K = (NUMBER OF NONZERO C \ ENTRIES IN THE UPPER TRIANGLE OF M)). C NVIRA \ ISP - INTEGER WORKING STORAGE DIVIDED UP INTO VARIOUS ARRAYS C \ NEEDED BY THE SUBROUTINES; ISP AND RSP SHOULD BE C \ EQUIVALENCED. C \ SIZE = LRATIO*NSP C FVIRA \ RSP - REAL WORKING STORAGE DIVIDED UP INTO VARIOUS ARRAYS C \ NEEDED BY THE SUBROUTINES; ISP AND RSP SHOULD BE C \ EQUIVALENCED. C \ SIZE = NSP. C NR \ ESP - IF NSP IS SUFFICIENTLY LARGE TO ALLOCATE SPACE, THEN C \ ESP IS SET TO THE AMOUNT OF EXCESS STORAGE PROVIDED. C C C CONVERSION TO DOUBLE PRECISION C C TO CONVERT THESE ROUTINES FOR DOUBLE PRECISION ARRAYS, SIMPLY USE C THE DOUBLE PRECISION DECLARATIONS IN PLACE OF THE REAL DECLARATIONS C IN EACH SUBPROGRAM; IN ADDITION, THE DATA VALUE OF THE INTEGER C VARIABLE LRATIO MUST BE SET AS INDICATED IN SUBROUTINE TDRV C