C ALGORITHM 664, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 14, NO. 3, P.257. C*********************************************************************** C C TESTDRIVER FOR THE SUBROUTINE GBSOL JANUARY 15, 1988 C C THIS PROGRAM IS A 1977 AMERICAN NATIONAL STANDARD FORTRAN IMPLE- C MENTATION OF A TESTDRIVER FOR THE SUBROUTINE GBSOL. C C THERE ARE SIX TESTS BASED ON THE FOLLOWING BAND MATRIX: C C C 20. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. 0. 0. 0. ... C C -1. 0. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. 0. 0. ... C C 2. -1. 20. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. 0. ... C C -3. 2. -1. 0. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. ... C C 4. -3. 2. -1. 20. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. ... C C -5. 4. -3. 2. -1. 0. 1. -2. 3. -4. 5. -6. 7. -8. 9. ... C C 6. -5. 4. -3. 2. -1. 20. 1. -2. 3. -4. 5. -6. 7. -8. ... C C -7. 6. -5. 4. -3. 2. -1. 0. 1. -2. 3. -4. 5. -6. 7. ... C C 8. -7. 6. -5. 4. -3. 2. -1. 20. 1. -2. 3. -4. 5. -6. ... C C -9. 8. -7. 6. -5. 4. -3. 2. -1. 0. 1. -2. 3. -4. 5. ... C C 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 20. 1. -2. 3. -4. ... C C 0. 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 0. 1. -2. 3. ... C C 0. 0. 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 20. 1. -2. ... C C 0. 0. 0. 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 0. 1. ... C C 0. 0. 0. 0. 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 20. ... C . . . . . . . . . . . . . . . C . . . . . . . . . . . . . . . C . . . . . . . . . . . . . . . C C C THE ROWS OF THIS MATRIX ARE CALCULATED BY THE SUBROUTINE ROW WHICH C IS SUPPLIED TOGETHER WITH THE TESTDRIVER. C C IN ALL TESTS, THE RIGHT HAND SIDE IS CHOSEN SO THAT ALL COMPONENTS C OF THE SOLUTION VECTOR HAVE THE VALUE ONE. C C IN THE FIRST FOUR TESTS A SYSTEM WITH ONE THOUSAND EQUATIONS IS C SOLVED. TO BEGIN WITH, THE SECOND DIMENSION NA OF A IS CHOSEN AS C 1000 SO THAT NO I/O OPERATIONS ARE PERFORMED. THEN, NA IS SET TO C THE VALUE OF 511 SO THAT THE MATRIX IS DIVIDED INTO TWO BLOCKS AND C THE FORTRAN 77 VERSIONS OF THE SUBROUTINES GOPEN, GCLOSE, AND C GRWRAN ARE TESTED (THE USER SHOULD REWRITE THESE SUBROUTINES BY C APPLYING ASSEMBLER UTILITY PROGRAMS OF HIS INSTALLATION. FORTRAN C SUPPORTED I/O IS TOO INEFFICIENT AND SHOULD BE USED ONLY FOR C TESTING). IN THE THIRD TEST, THE VALUE OF NA IS TOO SMALL, AND C GBSOL RETURNS TO THE CALLING PROGRAM WITH THE ERROR CODE "JERR=1." C IN THE FOURTH TEST THE PARAMETER EPS IS CHOSEN TO BE SLIGHTLY C LARGER THAN THE SMALLEST PIVOT SO THAT THE PIVOT IS REPLACED BY C EPS. C C FOR THE LAST TWO TESTS, THE FIRST LINE OF THE MATRIX IS CHANGED C IN ORDER TO PRODUCE A NEARLY SINGULAR MATRIX WITH A SMALL PIVOT: C C C 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. DEL 0. 0. 0. 0. ... C C -1. 0. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. 0. 0. ... C C 2. -1. 20. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. 0. ... C C -3. 2. -1. 0. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. ... C . . . . . . . . . . . . . . . C . . . . . . . . . . . . . . . C . . . . . . . . . . . . . . . C C C DEL IS SET TO THE VALUE OF 1.D-12. C C IN THE FIFTH TEST, THE DIMENSION N IS 750, NA IS SET TO 241, AND C THE SYSTEM IS SOLVED. FINALLY, THE PARAMETER EPS IS AGAIN CHOSEN C TO BE SLIGHTLY LARGER THAN THE SMALLEST PIVOT SO THAT THE PIVOT C IS REPLACED BY EPS. C C AUTHOR: GEZA SCHRAUF, CAL TECH C C*********************************************************************** PARAMETER (NP=1000,NAP=1000) C DOUBLE PRECISION A(21,NAP),B(NP,1) DOUBLE PRECISION ALOGDT,DEL,EPS,ERROR,PIVMAX,PIVMIN,SIGNDT INTEGER IFLAG,INFO,IPIV,JERR,L,M,MD,N,NA,NFSUA,NFSUAB,NOEOA, . NOEOB,NREAL8,NRHS,N1 C COMMON /TOROW/ DEL,IFLAG,N C IFLAG = 1 M = 21 MD = (M+1)/2 N = NP NA = NAP NRHS = 1 N1 = N - 1 C----------------------------------------------------------------------- C FIRST TEST: SOLVE THE SYSTEM WITHOUT DIRECT ACCESS I/O. C----------------------------------------------------------------------- WRITE(6,9001) NOEOB = M*N WRITE (6,9011) M,NA,N,NOEOB C NOEOA = M*NA NFSUA = 2*NOEOA NREAL8 = NOEOA + N NFSUAB = 2*NREAL8 C WRITE (6,9021) NFSUAB,NFSUA C C ASSIGN VALUES TO THE RIGHT HAND SIDE. C DO 10 L = 1,N,2 B(L,1) = 20.D0 B(L+1,1) = 0.D0 10 CONTINUE C B(1,1) = 15.D0 B(2,1) = -6.D0 B(3,1) = 16.D0 B(4,1) = -7.D0 B(5,1) = 17.D0 B(6,1) = -8.D0 B(7,1) = 18.D0 B(8,1) = -9.D0 B(9,1) = 19.D0 B(10,1) =-10.D0 C B(N-9,1) = 30.D0 B(N-8,1) = 1.D0 B(N-7,1) = 29.D0 B(N-6,1) = 2.D0 B(N-5,1) = 28.D0 B(N-4,1) = 3.D0 B(N-3,1) = 27.D0 B(N-2,1) = 4.D0 B(N-1,1) = 26.D0 B(N,1) = 5.D0 C EPS = 1.D-12 INFO = 0 C C SOLVE THE SYSTEM A*X = B. C CALL GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,EPS, . INFO,JERR) C WRITE (6,9041) ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,INFO,JERR C C CALCULATE THE ERROR IN L2 - NORM. ERROR = 0.D0 C DO 20 L = 1,N ERROR = ERROR + (B(L,1)-1.D0)**2 20 CONTINUE C ERROR = DSQRT(ERROR/N) C WRITE (6,9051) ERROR C CALCULATE THE ERROR IN MAXIMUM - NORM. ERROR = 0.D0 C DO 30 L = 1,N TEMP = DABS(B(L,1)-1.D0) IF (TEMP.GT.ERROR) ERROR = TEMP 30 CONTINUE C WRITE (6,9052) ERROR C----------------------------------------------------------------------- C SECOND TEST: SOLVE THE SYSTEM WITH DIRECT ACCESS I/O. C----------------------------------------------------------------------- WRITE(6,9002) C NA = 511 C NOEOB = M*N WRITE (6,9011) M,NA,N,NOEOB C NOEOA = M*NA NFSUA = 2*NOEOA NREAL8 = NOEOA + N NFSUAB = 2*NREAL8 C WRITE (6,9021) NFSUAB,NFSUA C C ASSIGN VALUES TO THE RIGHT HAND SIDE. C DO 40 L = 1,N,2 B(L,1) = 20.D0 B(L+1,1) = 0.D0 40 CONTINUE C B(1,1) = 15.D0 B(2,1) = -6.D0 B(3,1) = 16.D0 B(4,1) = -7.D0 B(5,1) = 17.D0 B(6,1) = -8.D0 B(7,1) = 18.D0 B(8,1) = -9.D0 B(9,1) = 19.D0 B(10,1) =-10.D0 C B(N-9,1) = 30.D0 B(N-8,1) = 1.D0 B(N-7,1) = 29.D0 B(N-6,1) = 2.D0 B(N-5,1) = 28.D0 B(N-4,1) = 3.D0 B(N-3,1) = 27.D0 B(N-2,1) = 4.D0 B(N-1,1) = 26.D0 B(N,1) = 5.D0 C EPS = 1.D-12 INFO = 0 C C SOLVE THE SYSTEM A*X = B. C CALL GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,EPS, . INFO,JERR) C WRITE (6,9041) ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,INFO,JERR C C CALCULATE THE ERROR IN L2 - NORM. ERROR = 0.D0 C DO 50 L = 1,N ERROR = ERROR + (B(L,1)-1.D0)**2 50 CONTINUE C ERROR = DSQRT(ERROR/N) C WRITE (6,9051) ERROR C CALCULATE THE ERROR IN MAXIMUM - NORM. ERROR = 0.D0 C DO 60 L = 1,N TEMP = DABS(B(L,1)-1.D0) IF (TEMP.GT.ERROR) ERROR = TEMP 60 CONTINUE C WRITE (6,9052) ERROR C----------------------------------------------------------------------- C THIRD TEST: VALUE OF NA IS TOO SMALL. C----------------------------------------------------------------------- WRITE(6,9003) C NA = 7 C NOEOB = M*N WRITE (6,9011) M,NA,N,NOEOB C C ASSIGN VALUES TO THE RIGHT HAND SIDE. C DO 70 L = 1,N,2 B(L,1) = 20.D0 B(L+1,1) = 0.D0 70 CONTINUE C B(1,1) = 15.D0 B(2,1) = -6.D0 B(3,1) = 16.D0 B(4,1) = -7.D0 B(5,1) = 17.D0 B(6,1) = -8.D0 B(7,1) = 18.D0 B(8,1) = -9.D0 B(9,1) = 19.D0 B(10,1) =-10.D0 C B(N-9,1) = 30.D0 B(N-8,1) = 1.D0 B(N-7,1) = 29.D0 B(N-6,1) = 2.D0 B(N-5,1) = 28.D0 B(N-4,1) = 3.D0 B(N-3,1) = 27.D0 B(N-2,1) = 4.D0 B(N-1,1) = 26.D0 B(N,1) = 5.D0 C EPS = 1.D-12 INFO = 0 C C SOLVE THE SYSTEM A*X = B. C CALL GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,EPS, . INFO,JERR) C WRITE (6,9043) JERR C C----------------------------------------------------------------------- C FOURTH TEST: SOLVE A SYSTEM WITH REPLACED SMALLEST PIVOT. C----------------------------------------------------------------------- WRITE(6,9004) C NA = 511 C NOEOB = M*N WRITE (6,9011) M,NA,N,NOEOB C NOEOA = M*NA NFSUA = 2*NOEOA NREAL8 = NOEOA + N NFSUAB = 2*NREAL8 C WRITE (6,9021) NFSUAB,NFSUA C C ASSIGN VALUES TO THE RIGHT HAND SIDE. C DO 80 L = 1,N,2 B(L,1) = 20.D0 B(L+1,1) = 0.D0 80 CONTINUE C B(1,1) = 15.D0 B(2,1) = -6.D0 B(3,1) = 16.D0 B(4,1) = -7.D0 B(5,1) = 17.D0 B(6,1) = -8.D0 B(7,1) = 18.D0 B(8,1) = -9.D0 B(9,1) = 19.D0 B(10,1) =-10.D0 C B(N-9,1) = 30.D0 B(N-8,1) = 1.D0 B(N-7,1) = 29.D0 B(N-6,1) = 2.D0 B(N-5,1) = 28.D0 B(N-4,1) = 3.D0 B(N-3,1) = 27.D0 B(N-2,1) = 4.D0 B(N-1,1) = 26.D0 B(N,1) = 5.D0 C EPS = 0.9667910D+1 INFO = 0 C C SOLVE THE SYSTEM A*X = B. C CALL GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,EPS, . INFO,JERR) C WRITE (6,9041) ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,INFO,JERR C C CALCULATE THE ERROR IN L2 - NORM. ERROR = 0.D0 C DO 90 L = 1,N ERROR = ERROR + (B(L,1)-1.D0)**2 90 CONTINUE C ERROR = DSQRT(ERROR/N) C WRITE (6,9051) ERROR C CALCULATE THE ERROR IN MAXIMUM - NORM. ERROR = 0.D0 C DO 100 L = 1,N TEMP = DABS(B(L,1)-1.D0) IF (TEMP.GT.ERROR) ERROR = TEMP 100 CONTINUE C WRITE (6,9052) ERROR C----------------------------------------------------------------------- C FIFTH TEST: SOLVE A SYSTEM WITH A SMALL PIVOT. C----------------------------------------------------------------------- WRITE(6,9005) C IFLAG = 2 C N = 750 NA = 241 C NOEOB = M*N WRITE (6,9011) M,NA,N,NOEOB C NOEOA = M*NA NFSUA = 2*NOEOA NREAL8 = NOEOA + N NFSUAB = 2*NREAL8 C WRITE (6,9021) NFSUAB,NFSUA C DEL = 1.D-12 WRITE (6,9031) DEL C ASSIGN VALUES TO THE RIGHT HAND SIDE. DO 110 L = 1,N,2 B(L,1) = 20.D0 B(L+1,1) = 0.D0 110 CONTINUE C B(1,1) = DEL B(2,1) = -6.D0 B(3,1) = 16.D0 B(4,1) = -7.D0 B(5,1) = 17.D0 B(6,1) = -8.D0 B(7,1) = 18.D0 B(8,1) = -9.D0 B(9,1) = 19.D0 B(10,1) =-10.D0 C B(N-9,1) = 30.D0 B(N-8,1) = 1.D0 B(N-7,1) = 29.D0 B(N-6,1) = 2.D0 B(N-5,1) = 28.D0 B(N-4,1) = 3.D0 B(N-3,1) = 27.D0 B(N-2,1) = 4.D0 B(N-1,1) = 26.D0 B(N,1) = 5.D0 C EPS = 1.D-12 WRITE (6,9032) EPS INFO = 0 C C SOLVE THE SYSTEM A*X = B. C CALL GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,EPS, . INFO,JERR) C WRITE (6,9041) ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,INFO,JERR C C CALCULATE THE ERROR IN L2 - NORM. ERROR = 0.D0 C DO 120 L = 1,N ERROR = ERROR + (B(L,1)-1.D0)**2 120 CONTINUE C ERROR = DSQRT(ERROR/N) C WRITE (6,9051) ERROR C C CALCULATE THE ERROR IN MAXIMUM - NORM. ERROR = 0.D0 C DO 130 L = 1,N TEMP = DABS(B(L,1)-1.D0) IF (TEMP.GT.ERROR) ERROR = TEMP 130 CONTINUE C WRITE (6,9052) ERROR C----------------------------------------------------------------------- C SIXTH TEST: SOLVE A SYSTEM WITH REPLACED SMALLEST PIVOT. C----------------------------------------------------------------------- WRITE(6,9006) C IFLAG = 2 C N = 750 NA = 241 C NOEOB = M*N WRITE (6,9011) M,NA,N,NOEOB C NOEOA = M*NA NFSUA = 2*NOEOA NREAL8 = NOEOA + N NFSUAB = 2*NREAL8 C WRITE (6,9021) NFSUAB,NFSUA C DEL = 1.D-12 WRITE (6,9031) DEL C ASSIGN VALUES TO THE RIGHT HAND SIDE. DO 140 L = 1,N,2 B(L,1) = 20.D0 B(L+1,1) = 0.D0 140 CONTINUE C B(1,1) = DEL B(2,1) = -6.D0 B(3,1) = 16.D0 B(4,1) = -7.D0 B(5,1) = 17.D0 B(6,1) = -8.D0 B(7,1) = 18.D0 B(8,1) = -9.D0 B(9,1) = 19.D0 B(10,1) =-10.D0 C B(N-9,1) = 30.D0 B(N-8,1) = 1.D0 B(N-7,1) = 29.D0 B(N-6,1) = 2.D0 B(N-5,1) = 28.D0 B(N-4,1) = 3.D0 B(N-3,1) = 27.D0 B(N-2,1) = 4.D0 B(N-1,1) = 26.D0 B(N,1) = 5.D0 C EPS = 0.22837D-9 WRITE (6,9032) EPS INFO = 0 C C SOLVE THE SYSTEM A*X = B. C CALL GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,EPS, . INFO,JERR) C WRITE (6,9041) ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,INFO,JERR C C CALCULATE THE ERROR IN L2 - NORM. ERROR = 0.D0 C DO 150 L = 1,N ERROR = ERROR + (B(L,1)-1.D0)**2 150 CONTINUE C ERROR = DSQRT(ERROR/N) C WRITE (6,9051) ERROR C CALCULATE THE ERROR IN MAXIMUM - NORM. ERROR = 0.D0 C DO 160 L = 1,N TEMP = DABS(B(L,1)-1.D0) IF (TEMP.GT.ERROR) ERROR = TEMP 160 CONTINUE C WRITE (6,9052) ERROR C STOP C 9001 FORMAT(//' ',71('-')/' F I R S T T E S T W I T H O U T ', . ' D I R E C T - A C C E S S I / O'/' ',71('-')//) 9002 FORMAT(//' ',68('-')/' S E C O N D T E S T W I T H ', . ' D I R E C T - A C C E S S I / O'/' ',68('-')//) 9003 FORMAT(//' ',71('-')/' T H I R D T E S T W I T H T O O ', . ' S M A L L V A L U E OF NA'/' ',71('-')//) 9004 FORMAT(//' ',79('-')/' F O U R T H T E S T W I T H ', . ' S M A L L E S T P I V O T R E P L A C E D'/ . ' ',79('-')//) 9005 FORMAT(//' ',69('-')/' F I F T H T E S T W I T H ', . ' S M A L L P I V O T E L E M E N T'/' ',69('-')//) 9006 FORMAT(//' ',77('-')/' S I X T H T E S T W I T H ', . ' S M A L L E S T P I V O T R E P L A C E D'/ . ' ',77('-')//) 9011 FORMAT (//' FIRST DIMENSION OF A: ', . ' M =',I10/ . ' SECOND DIMENSION OF A: ', . ' NA =',I10// . ' NUMBER OF EQUATIONS: ', . ' N =',I10/ . ' NUMBER OF ELEMENTS IN THE BAND OF THE MATRIX A: ', . ' M*N =',I10///) 9021 FORMAT (' MEMORY REQUIREMENTS:'// . ' NUMBER OF FORTRAN STORAGE UNITS NEEDED '/ . ' TO STORE THE MATRIX A AND THE RHS B: ', . ' 2*(M*NA+N) =',I10// . ' NUMBER OF FORTRAN STORAGE UNITS NEEDED '/ . ' TO STORE THE MATRIX A: ', . ' 2*(M*NA) =',I10//) 9031 FORMAT (' DEL = ',E14.7/) 9032 FORMAT (' EPS = ',E14.7) 9041 FORMAT (//' ALOGDT = ',E14.7/' SIGNDT = ',F3.0//' PIVMAX = ', . E14.7/' PIVMIN = ',E14.7/' IPIV = ',I10//' INFO = ', . I10//' JERR = ',I10) 9043 FORMAT (//' JERR = ',I10// . ' RETURN TO CALLING PROGRAM WITH ERROR CODE', . ' "JERR = 1,"'/' BECAUSE NA = 7 ', . ' IS LESS THAN (M+1)/2 + 2 = 13.') 9051 FORMAT (//' ERROR IN THE L2-NORM: ',E15.7) 9052 FORMAT (' AND IN THE MAXIMUM-NORM:',E15.7//) END C C C C C SUBROUTINE ROW(AROW,L) DOUBLE PRECISION AROW(*),DEL,ATEMP INTEGER I,IFLAG,ISIGN,ITEMP,L,MD,MDMI,MDPI,N C COMMON /TOROW/ DEL,IFLAG,N C ISIGN = -1 MD = 11 C ITEMP = L/2 ATEMP = L - ITEMP*2 C AROW(MD) = 20.D0 * ATEMP C DO 10 I = 1,10 ISIGN = ISIGN* (-1) MDMI = MD - I MDPI = MD + I C AROW(MDMI) = -ISIGN*I AROW(MDPI) = ISIGN*I 10 CONTINUE C IF (L.GT.10) GO TO 30 C IF (L.LE.10) AROW(1) = 0.D0 IF (L.LE.9) AROW(2) = 0.D0 IF (L.LE.8) AROW(3) = 0.D0 IF (L.LE.7) AROW(4) = 0.D0 IF (L.LE.6) AROW(5) = 0.D0 IF (L.LE.5) AROW(6) = 0.D0 IF (L.LE.4) AROW(7) = 0.D0 IF (L.LE.3) AROW(8) = 0.D0 IF (L.LE.2) AROW(9) = 0.D0 IF (L.GT.1) GO TO 30 AROW(10) = 0.D0 IF (IFLAG.EQ.1) GO TO 30 DO 20 I=11,20 AROW(I) = 0.D0 20 CONTINUE AROW(21) = DEL C 30 IF (L.LT.N-9) GO TO 40 C IF (L.GE.N) AROW(12) = 0.D0 IF (L.GE.N-1) AROW(13) = 0.D0 IF (L.GE.N-2) AROW(14) = 0.D0 IF (L.GE.N-3) AROW(15) = 0.D0 IF (L.GE.N-4) AROW(16) = 0.D0 IF (L.GE.N-5) AROW(17) = 0.D0 IF (L.GE.N-6) AROW(18) = 0.D0 IF (L.GE.N-7) AROW(19) = 0.D0 IF (L.GE.N-8) AROW(20) = 0.D0 IF (L.GE.N-9) AROW(21) = 0.D0 40 CONTINUE C RETURN END C*********************************************************************** C C NAME: GBSOL FORTRAN 77 JANUARY 15, 1988 C C PURPOSE: C THIS ALGORITHM IS A 1977 AMERICAN NATIONAL STANDARD FORTRAN C IMPLEMENTATION OF A GAUSSIAN ELIMINATION PROCEDURE WITH PARTIAL C PIVOTING TO SOLVE SYSTEMS WITH LARGE BANDED MATRICES AND SEVERAL C RIGHT-HAND SIDES. C C THE ALGORITHM KEEPS ONLY THE DATA NECESSARY FOR THE ACTUAL COMPU- C TATION IN THE MEMORY AND SWAPS THE OTHER DATA ON A DISK. IN ORDER C TO MINIMIZE I/O OPERATIONS, THE ELEMENTS OF THE MATRIX ARE CALCU- C LATED WHEN THEY ARE NEEDED FOR THE COMPUTATION OF THE UPPER TRIAN- C GULAR MATRIX. HENCE, THE USER HAS TO PROVIDE A SUBROUTINE WHICH C CALCULATES THE ROWS OF THE MATRIX. C C THE WHOLE PROGRAM CONFORMS TO STANDARD, BUT IT IS RECOMMENDED THAT C THE SUBROUTINES FOR HANDLING RANDOM-ACCESS STORAGE (GOPEN, GCLOSE, C AND GRWRAN) BE REPLACED BY INSTALLATION-DEPENDENT VERSIONS. C C THE ELEMENTS OF THE BANDMATRIX ARE STORED IN ROWS AS ILLUSTRATED C IN THE FOLLOWING EXAMPLE: C C C C * * 31 41 51 0 0 0 0 0 0 C * 22 32 42 52 0 0 0 0 0 C 13 23 33 43 53 0 0 0 0 C 0 14 24 34 44 54 0 0 0 C 0 0 15 25 35 45 55 0 0 C 0 0 0 16 26 36 46 56 0 C 0 0 0 0 17 27 37 47 57 C 0 0 0 0 0 18 28 38 48 * C 0 0 0 0 0 0 19 29 39 * * C C C THE FIRST INDEX SHOWS THE POSITION OF THE ELEMENT IN THE BAND, C AND THE SECOND ONE INDICATES THE ROW. ONLY THE BAND IS STORED, C AND THE USER HAS TO SUPPLY ZEROS FOR THE POSITIONS INDICATED WITH C AN ASTERISK IN THE ABOVE EXAMPLE. THE MATRIX A DECLARED IN GBSOL C CONTAINS THE FOLLOWING ELEMENTS: C C 0 0 31 41 51 C 0 22 32 42 52 C 13 23 33 43 53 C 14 24 34 44 54 C 15 25 35 45 55 C 16 26 36 46 56 C 17 27 37 47 57 C 18 28 38 48 0 C 19 29 39 0 0 C C THE USER HAS TO PROVIDE A SUBROUTINE "ROW" THAT CALCULATES THE C ROWS OF A. THIS SUBROUTINE IS CALLED BY GBSOL IN THE FORM C C CALL ROW( A(1,L),L ) C C WHERE A(1,L) POINTS TO THE ADDRESS OF THE FIRST ELEMENT IN THE C L-TH ROW. C C THE ELEMENTS OF THE ROWS OF A ARE STORED IN CONSECUTIVE ADDRESSES C SO THAT THE PROGRAM CAN EASILY BE MODIFIED FOR MACHINES WITH PIPE- C LINE ARCHITECTURE, SUCH AS THE CYBER 205 OR THE FPS164. THE MODI- C FICATION CONSISTS OF REPLACING THE LOOPS MARKED IN THE CODE BY C VECTOR OPERATIONS. C C CALLING SEQUENCE: C C CALL GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV, C EPS,INFO,JERR) C PARAMETERS: C A: DOUBLE PRECISION. A(M,NA), IS THE ARRAY DECLARED IN C THE SUBROUTINE TO PROVIDE MEMORY FOR THE COMPUTATION. C M: INTEGER, INPUT, NUMBER OF ELEMENTS IN EACH ROW OF THE C BAND. C NA: INTEGER, INPUT, SECOND DIMENSION OF A. C BY SPECIFYING NA, THE USER CHOOSES THE SIZE OF THE C PART OF THE MATRIX TO BE KEPT IN THE MEMORY: C NA>=N: THE WHOLE MATRIX IS KEPT IN THE MEMORY. C THERE ARE NO I/O OPERATIONS. C NA = (M+1)/2 + 2. C B: DOUBLE PRECISION. B(N,NRHS), CONTAINS THE RIGHT-HAND C SIDES AS INPUT AND THE SOLUTION VECTORS AS OUTPUT. C N: INTEGER, INPUT, NUMBER OF EQUATIONS. C NRHS: INTEGER, INPUT, NUMBER OF RIGHT-HAND SIDES. C ALOGDT: DOUBLE PRECISION, OUTPUT, LOGARITHM OF THE ABSOLUTE C VALUE OF THE DETERMINANT. C SIGNDT: DOUBLE PRECISION, OUTPUT, SIGN OF THE DETERMINANT. C THE DETERMINANT ITSELF CAN BE CALCULATED BY C DET( A ) = SIGNDT * 10.0D0**ALOGDT. C PIVMAX: DOUBLE PRECISION, OUTPUT, MAXIMUM OF THE ABSOLUTE C VALUES OF THE PIVOT ELEMENTS. C PIVMIN: DOUBLE PRECISION, OUTPUT, MINIMUM OF THE ABSOLUTE C VALUES OF THE PIVOT ELEMENTS. C IPIV: INTEGER, OUTPUT, POSITION OF THE PIVOT WITH MINIMAL C ABSOLUTE VALUE. C EPS: DOUBLE PRECISION, INPUT. ALL PIVOT ELEMENTS WHOSE C ABSOLUTE VALUES ARE LESS THAN EPS ARE REPLACED BY C SIGN( PIVOT ) * EPS. C INFO: INTEGER,OUTPUT, C = 0, NO PIVOT ELEMENT IS REPLACED. C = L, POSITION OF THE LAST ELEMENT WHICH IS REPLACED. C C FOR NORMAL USE WE RECOMMEND SETTING C EPS = DSQRT( SMALLEST DOUBLE PRECISION CONSTANT ) C IN ORDER TO AVOID OVERFLOW IN CASE OF A SINGULAR C MATRIX. WITH THIS SETTING AND A PROPERLY SCALED MATRIX C (I.E. PIVMAX = O(1)), A RETURNED VALUE OF INFO DIFFER- C ENT FROM ZERO INDICATES THAT THE SOLUTIONS HAVE NO C SIGNIFICANT FIGURES. C C JERR: INTEGER, OUTPUT,ERROR CODE: C JERR=0: NA >= (M+1)/2+2. C JERR=1: NA < (M+1)/2+2. CHOOSE A LARGER VALUE FOR NA. C C USER SUPPLIED SUBROUTINES CALLED: C C ROW CALCULATES A ROW OF THE BAND OF MATRIX AS DESCRIBED C ABOVE. C C GBSOL CONTAINS THE FOLLOWING CALLS IN ORDER TO PERFORM THE C INSTALLATION-DEPENDENT DIRECT-ACCESS IN- AND OUTPUT OPERATIONS: C C C CALL GOPEN(M) TO OPEN A TEMPORARY FILE FOR DIRECT-ACCESS I/O. C C CALL GCLOSE TO CLOSE AND FREE THE FILE OPENED WITH GOPEN. C C CALL GRWRAN(IOP,A,NREAL8,NSTART,M) C C IOP=1: TO READ NREAL8 DOUBLE PRECISION NUMBERS C FROM A DISK TO THE ARRAY A, STARTING WITH C THE NSTART-TH NUMBER ON THE DISK, OR C C IOP=2: TO WRITE THE FIRST NREAL8 DOUBLE PRECISION C NUMBERS OF THE ARRAY A ON A DISK, STARTING C WITH THE POSITION NSTART. C C THE USER SHOULD APPLY UTILITY PROGRAMS WRITTEN IN ASSEMBLER IN C ORDER TO GUARANTEE RAPID I/O OPERATIONS. STANDARD FORTRAN VERSIONS C OF THE ABOVE SUBROUTINES ARE SUPPLIED, BUT FORTRAN SUPPLIED I/O IS C TOO SLOW. C C AUTHOR: GEZA SCHRAUF, CAL TECH C C*********************************************************************** SUBROUTINE GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV, . EPS,INFO,JERR) C DOUBLE PRECISION A(M,NA),B(N,NRHS) DOUBLE PRECISION ABSPIV,AH,ALOGDT,EPS,PIV,PIVMAX,PIVMIN,RA1L,RA1N, . RPIV,S,SIGNDT INTEGER I,IH,INFO,IPIV,ISTART,J,JA,JB,JEND,JERR,JH,JPREL,JSTART, . K,L,LABS, LABS1,LAMLR,LB,LBLK,LBLKM1,LDONE,LEND,LENGTH, . LREL,LRELP1,LREL1,LSTABS,LSTREL,M,MD,MD1,M1,N,NA,NBLKM1, . NBLKS,NEWABS,NEWREL,NREAL8,NRHS,NRLB,NRPB,NRPBPJ,NRPS, . NRPSM1,NSTART,NTEMP,N1 ALOGDT = 0.D0 SIGNDT = 1.D0 PIVMAX = 0.D0 INFO = 0 IPIV = 1 JERR = 0 C MD = (M+1)/2 MD1 = MD - 1 M1 = M - 1 N1 = N - 1 C IF (NA.LT.N) GO TO 10 C NRPB = N NRLB = N NBLKS = 1 NBLKM1 = NBLKS - 1 GO TO 30 C 10 NRPB = NA - MD NBLKS = N/NRPB NTEMP = NBLKS*NRPB IF(N.NE.NTEMP) NBLKS = NBLKS + 1 C NBLKM1 = NBLKS - 1 NRLB = N - NBLKM1*NRPB C IF (NRPB.GE.2) GO TO 20 JERR = 1 RETURN 20 CONTINUE C NREAL8 = M * NRPB C OPEN TEMPORARY FILE FOR DIRECT ACCESS. CALL GOPEN(M) 30 CONTINUE C C CALCULATE THE FIRST MD ROWS OF THE MATRIX A. DO 40 L = 1,MD CALL ROW(A(1,L),L) 40 CONTINUE C SHIFT THE MATRIX ELEMENTS OF THE FIRST MD1 ROWS TO THE LEFT. IH = MD ISTART = IH C DO 70 J = 1,MD1 DO 50 I = 1,IH A(I,J) = A(ISTART+I-1,J) 50 CONTINUE IH = IH + 1 DO 60 I = IH,M A(I,J) = 0.D0 60 CONTINUE ISTART = ISTART - 1 70 CONTINUE C C INITIALIZE PIVMIN. C PIVMIN = 0.D0 C DO 75 J=1,MD AH = DABS(A(1,J)) IF (AH.GT.PIVMIN) PIVMIN = AH 75 CONTINUE C C CALCULATION OF THE UPPER TRIANGULAR MATRIX. C DO 200 LBLK = 1,NBLKS LEND = NRPB IF (LBLK.EQ.NBLKS) LEND = NRLB - 1 C LDONE = (LBLK-1)*NRPB C DO 170 L = 1,LEND LREL = L LABS = LDONE + LREL LREL1 = LREL - 1 LABS1 = LABS - 1 LSTABS = MIN0(LABS+MD1,N) NRPS = LSTABS - LABS1 C PIVOT SEARCH. ABSPIV = 0.D0 JPREL = 1 DO 80 J = 1,NRPS AH = DABS(A(1,J+LREL1)) IF (AH.LE.ABSPIV) GO TO 80 JPREL = J ABSPIV = AH 80 CONTINUE C INTERCHANGE THE ROWS OF A AND THE CORRESPONDING RHS'S. IF (JPREL.EQ.1) GO TO 110 C JA = JPREL + LREL1 JB = JPREL + LABS1 C DO 90 I = 1,M S = A(I,JA) A(I,JA) = A(I,LREL) A(I,LREL) = S 90 CONTINUE DO 100 I = 1,NRHS S = B(JB,I) B(JB,I) = B(LABS,I) B(LABS,I) = S 100 CONTINUE C 110 PIV = A(1,LREL) IF (DABS(PIV).GE.EPS) GO TO 120 IF (PIV.NE.0.D0) PIV = DSIGN(EPS,PIV) IF (PIV.EQ.0.D0) PIV = EPS A(1,LREL) = PIV INFO = LABS 120 CONTINUE C IF (JPREL.NE.1) SIGNDT = -SIGNDT IF (PIV.LT.0.D0) SIGNDT = -SIGNDT ABSPIV = DABS(PIV) ALOGDT = ALOGDT + DLOG10(ABSPIV) IF (ABSPIV.GT.PIVMAX) PIVMAX = ABSPIV IF (ABSPIV.GT.PIVMIN) GO TO 130 PIVMIN = ABSPIV IPIV = LABS 130 CONTINUE C C ADD A SUITABLE MULTIPLE OF THE LREL-TH ROW OF A TO THE FOLLOWING C ROWS OF A, AND SHIFT THE MATRIX ROWS TO THE LEFT, SO THAT THE C NEXT PIVOT IS TO BE SEARCHED AGAIN IN THE FIRST COLUMN OF A. C ADD THE SAME MULTIPLE OF THE LABS-TH ROW OF B TO THE FOLLOWING C ROWS OF B. C NRPSM1 = NRPS - 1 RPIV = 1.D0/PIV C JSTART = LABS + 1 JEND = LABS + NRPSM1 LAMLR = LABS - LREL C DO 150 I = 1,NRHS S = -B(LABS,I)*RPIV C ************ VECTORIZE THE FOLLOWING LOOP, IF POSSIBLE *********** DO 161 J = JSTART,JEND B(J,I) = A(1,J-LAMLR)*S + B(J,I) 161 CONTINUE 150 CONTINUE C DO 160 JH = 1,NRPSM1 JA = LREL + JH JB = LABS + JH S = -A(1,JA)*RPIV C ************ VECTORIZE THE FOLLOWING LOOP, IF POSSIBLE *********** DO 140 I = 1,M1 A(I,JA) = A(I+1,LREL)*S + A(I+1,JA) 140 CONTINUE A(M,JA) = 0.D0 160 CONTINUE C NEWABS = LABS + MD NEWREL = NEWABS - LDONE C IF (NEWABS.LE.N) CALL ROW(A(1,NEWREL),NEWABS) 170 CONTINUE C IF (LBLK.EQ.NBLKS) GO TO 200 C C WRITE THE NREAL8 = M*NRPB NUMBERS JUST CALCULATED ON THE TEMPO- C RARY FILE. C NSTART = LDONE*M + 1 CALL GRWRAN(2,A,NREAL8,NSTART,M) C SHIFT THE REMAINING MD ROWS UP. DO 190 J = 1,MD NRPBPJ = NRPB + J DO 180 I = 1,M A(I,J) = A(I,NRPBPJ) 180 CONTINUE 190 CONTINUE C 200 CONTINUE C THE LAST ELEMENT. LRELP1 = LREL + 1 PIV = A(1,LRELP1) IF (DABS(PIV).GE.EPS) GO TO 210 IF (PIV.NE.0.D0) PIV = DSIGN(EPS,PIV) IF (PIV.EQ.0.D0) PIV = EPS A(1,LRELP1) = PIV INFO = N 210 CONTINUE C IF (PIV.LT.0.D0) SIGNDT = -SIGNDT ABSPIV = DABS(PIV) ALOGDT = ALOGDT + DLOG10(ABSPIV) IF (ABSPIV.GT.PIVMAX) PIVMAX = ABSPIV IF (ABSPIV.GT.PIVMIN) GO TO 220 PIVMIN = ABSPIV IPIV = N 220 CONTINUE C C SOLVE THE SYSTEMS WITH THE UPPER TRIANGULAR MATRIX. C RA1N = 1.D0/A(1,LRELP1) DO 230 I = 1,NRHS B(N,I) = B(N,I)*RA1N 230 CONTINUE C DO 290 LB = 1,NBLKS LBLK = NBLKS + 1 - LB LBLKM1 = LBLK - 1 LEND = NRPB C IF (LBLK.NE.NBLKS) GO TO 240 LEND = NRLB - 1 GO TO 250 * 240 CONTINUE C READ NREAL8 = M*NRPB NUMBERS FROM THE TEMPORARY FILE. NSTART = LBLKM1*NREAL8 + 1 CALL GRWRAN(1,A,NREAL8,NSTART,M) C 250 CONTINUE C BACKSOLVE. DO 280 L = 1,LEND LREL = LEND + 1 - L LABS = LBLKM1*NRPB + LREL LENGTH = M1 IF (LBLK.EQ.NBLKS) LENGTH = MIN0(L,M1) RA1L = 1.D0/A(1,LREL) C DO 270 K = 1,NRHS S = 0.D0 C ************ VECTORIZE THE FOLLOWING LOOP, IF POSSIBLE *********** DO 260 I = 1,LENGTH S = A(I+1,LREL)*B(I+LABS,K) + S 260 CONTINUE C B(LABS,K) = (B(LABS,K)-S)*RA1L 270 CONTINUE 280 CONTINUE C 290 CONTINUE C CLOSE AND FREE TEMPORARY FILE. IF (NA.EQ.N) GO TO 300 CALL GCLOSE 300 RETURN END C C C C C C*********************************************************************** C C NAME: GOPEN FORTRAN 77 JANUARY 15, 1988 C C PURPOSE: C THE SUBROUTINE "GOPEN" OPENS A TEMPORARY FILE FOR DIRECT-ACCESS. C C CALLING SEQUENCE: CALL GOPEN(M) C C PARAMETERS: C M: INTEGER, INPUT, NUMBER OF DOUBLE-PRECISION NUMBERS IN C EACH DIRECT-ACCESS RECORD. USED ONLY FOR THE STANDARD C FORTRAN 77 VERSION. C C AUTHOR: GEZA SCHRAUF, CAL TECH C C*********************************************************************** C SUBROUTINE GOPEN(M) INTEGER IUNIT,IRECL,M IUNIT = 9 C C----------------------------------------------------------------------- C C STANDARD FORTRAN 77 VERSION. THIS VERSION SHOULD BE REPLACED BY C A FASTER MACHINE-DEPENDENT ROUTINE. C C ******** CHOOSE THE RECORD LENGTH "IRECL" SO THAT ONE ******** C ******** RECORD CONTAINS M DOUBLE PRECISION NUMBERS. ******** C IRECL = M * 8 C OPEN(IUNIT,STATUS='SCRATCH',ACCESS='DIRECT',RECL=IRECL) C C----------------------------------------------------------------------- C C INSTALLATION FOR THE IBM 4331 OF THE DEPARTMENT OF APPLIED MATHE- C MATICS AT THE CALIFORNIA INSTITUTE OF TECHNOLOGY. C C IBM FORTRAN UTILITIES FOR VM/370 *** (PROGRAM NUMBER: 5798-DFH) C C UTILITY USED: UOPEN IN FORTRAN 77 LEVEL LANGUAGE (MANUAL P. 34). C C- CALL UOPEN(IUNIT,'TEMP DATA A',5,3,0,0,' ',0,IERR) C C- IF (IERR.EQ.0) GO TO 10 C C- WRITE (6,9001) IERR C9001 FORMAT (///' AN ERROR OCCURRED USING THE IBM FORTRAN UTILITY', C- . ' "UOPEN" WHICH WAS '/' CALLED BY "GOPEN":'//' IBM ERROR CODE:', C- . I4//' 0: OPENED AS REQUESTED.'/ C- . ' 1: FILE ALREADY OPENED WITH ANOTHER UNIT NUMBER.'/ C- . ' 2: SHOULD NOT OCCUR.'/ C- . ' 3: UNIT NUMBER ALREADY ASSOCIATED WITH A FILE.'/ C- . ' 8: DATA INPUT ERROR. FILE NOT OPENED.'///) C C- 10 CONTINUE C C----------------------------------------------------------------------- C RETURN END C C C C C C*********************************************************************** C C NAME: GCLOSE FORTRAN 77 JANUARY 15, 1988 C C PURPOSE: C THE SUBROUTINE "GCLOSE" CLOSES AND FREES A TEMPORARY FILE PREVI- C OUSLY OPENED WITH "GOPEN." C C CALLING SEQUENCE: CALL GCLOSE C C AUTHOR: GEZA SCHRAUF, CAL TECH C C*********************************************************************** C SUBROUTINE GCLOSE INTEGER IUNIT IUNIT = 9 C C----------------------------------------------------------------------- C C STANDARD FORTRAN 77 VERSION. THIS VERSION SHOULD BE REPLACED BY C A MACHINE-DEPENDENT ROUTINE. C C CLOSE(IUNIT) C C----------------------------------------------------------------------- C C INSTALLATION FOR THE IBM 4331 OF THE DEPARTMENT OF APPLIED MATHE- C MATICS AT THE CALIFORNIA INSTITUTE OF TECHNOLOGY. C C IBM FORTRAN UTILITIES FOR VM/370 *** (PROGRAM NUMBER: 5798-DFH) C C UTILITY USED: UCLOSE IN FORTRAN 77 LEVEL LANGUAGE (MANUAL P. 32). C C- CALL UCLOSE(IUNIT,2,IERR) C C- IF (IERR.EQ.0) GO TO 10 C C- WRITE (6,9001) IERR C9001 FORMAT (///' AN ERROR OCCURRED USING THE IBM FORTRAN UTILITY', C- . ' "UCLOSE" WHICH WAS '/' CALLED BY "GCLOSE":'// C- . ' IBM ERROR CODE:',I4//' 0: CLOSED AS REQUESTED.'/ C- . ' 4: UNIT IS NOT OPENED.'/' 8: DATA INPUT ERROR.'///) C C- 10 CONTINUE C C----------------------------------------------------------------------- C RETURN END C C C C C C*********************************************************************** C C NAME: GRWRAN FORTRAN 77 JANUARY 15, 1988 C C PURPOSE: C THE SUBROUTINE "GRWRAN" ALLOWS THE USER TO READ FROM OR WRITE TO A C FILE ON A DISK USING A DIRECT-ACCESS METHOD. C C CALLING SEQUENCE: CALL GRWRAN(IOP,ARRAY,NREAL8,NSTART,M) C C PARAMETERS: C IOP: INTEGER, OPERATION MODE: C IOP = 1: READ FROM THE FILE TO THE ARRAY "ARRAY." C IOP = 2: WRITE FROM THE ARRAY "ARRAY" TO THE FILE. C ARRAY: DOUBLE PRECISION, ARRAY TO OR FROM WHICH DATA IS C TRANSFERRED. C NREAL8: INTEGER, NUMBER OF DOUBLE-PRECISION NUMBERS TO BE C TRANSFERRED. C NSTART: INTEGER, THE FIRST DOUBLE-PRECISION NUMBER ON THE DISK C WHICH IS TO BE TRANSFERRED. C M: INTEGER, INPUT, NUMBER OF DOUBLE PRECISION NUMBERS IN C EACH DIRECT-ACCESS RECORD. ONLY USED FOR THE STANDARD C FORTRAN 77 VERSION. C C AUTHOR: GEZA SCHRAUF, CAL TECH C C*********************************************************************** C SUBROUTINE GRWRAN(IOP,ARRAY,NREAL8,NSTART,M) DOUBLE PRECISION ARRAY(NREAL8) INTEGER IOP,IUNIT,K,KFIRST,KLAST,L,LREC,LSTART, . M,NRPB,NREAL8,NSTART IUNIT = 9 C C----------------------------------------------------------------------- C C STANDARD FORTRAN 77 VERSION. THIS VERSION SHOULD BE REPLACED BY C A FASTER MACHINE-DEPENDENT ROUTINE. C LSTART = (NSTART - 1)/M NRPB = NREAL8/M C IF (IOP.EQ.1) THEN DO 10 L=1,NRPB LREC = LSTART + L KFIRST = (L-1)*M + 1 KLAST = KFIRST + M - 1 READ(IUNIT,REC=LREC) (ARRAY(K),K=KFIRST,KLAST) 10 CONTINUE ELSE IF (IOP.EQ.2) THEN DO 20 L=1,NRPB LREC = LSTART + L KFIRST = (L-1)*M + 1 KLAST = KFIRST + M - 1 WRITE(IUNIT,REC=LREC) (ARRAY(K),K=KFIRST,KLAST) 20 CONTINUE ENDIF C C----------------------------------------------------------------------- C C ACCELERATED STANDARD FORTRAN 77 VERSION WITH LARGER RECORD LENGTH. C NREAL8 HAS TO BE S M A L L E R THAN THE MAXIMAL POSSIBLE RECORD C LENGTH. TO USE THIS VERSION REPLACE THE STATEMENT "CALL GOPEN(M)" C THE BY STATEMENT "CALL GOPEN(NREAL8)" IN SUBROUTINE GBSOL. C C- LREC = (NSTART - 1)/NREAL8 + 1 C C- IF (IOP.EQ.1) THEN C- READ(IUNIT,REC=LREC) ARRAY C- ELSE IF (IOP.EQ.2) THEN C- WRITE(IUNIT,REC=LREC) ARRAY C- ENDIF C C----------------------------------------------------------------------- C C INSTALLATION FOR THE IBM 4331 OF THE DEPARTMENT OF APPLIED MATHE- C MATICS AT THE CALIFORNIA INSTITUTE OF TECHNOLOGY. C C IBM FORTRAN UTILITIES FOR VM/370 *** (PROGRAM NUMBER: 5798-DFH) C C UTILITY USED: RWFIL IN FORTRAN 77 LEVEL LANGUAGE (MANUAL P. 25). C C- IFIRST = 2*NSTART - 1 C- IWORDS = 2*NREAL8 C C- CALL RWFIL(IUNIT,IOP,IFIRST,ARRAY(1),IWORDS,IERR) C C- IF (IERR.EQ.0) GO TO 10 C C- WRITE (6,9001) IERR,IOP,NREAL8,NSTART C9001 FORMAT (///' AN ERROR OCCURRED USING THE IBM FORTRAN UTILITY', C- . ' "RWFIL" WHICH WAS '/' CALLED BY "GRWRAN":'// C- . ' IBM ERROR CODE:',I4//' 0: NO ERROR.'/ C- . ' 5: END OF FILE ENCOUNTERED.'/ C- . ' 6: UNIT HAS NOT BEEN OPENED AS A DIRECT-ACCESS FILE.'/ C- . ' 8: DATA INPUT ERROR.'//' PARAMETER VALUES OF "GRWRAN":'// C- . ' OPERATION CODE IOP =', C- . I4/' NUMBER OF REAL*8 NUMBERS NREAL8 =', C- . I10/' FIRST REAL*8 NUMBER NSTART =',I10///) C C- 10 CONTINUE C C----------------------------------------------------------------------- C RETURN END * * ----------------------------------------------------------------------- F I R S T T E S T W I T H O U T D I R E C T - A C C E S S I / O ----------------------------------------------------------------------- * * * * FIRST DIMENSION OF A: M = 21 SECOND DIMENSION OF A: NA = 1000 * NUMBER OF EQUATIONS: N = 1000 NUMBER OF ELEMENTS IN THE BAND OF THE MATRIX A: M*N = 21000 * * * MEMORY REQUIREMENTS: * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A AND THE RHS B: 2*(M*NA+N) = 44000 * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A: 2*(M*NA) = 42000 * * * * CALOGDT = 0.1211436E+04 SIGNDT = 1. * PIVMAX = 0.3727772E+02 PIVMIN = 0.9667904E+01 IPIV = 994 * INFO = 0 * JERR = 0 * * ERROR IN THE L2-NORM: 0.1497275E-14 AND IN THE MAXIMUM-NORM: 0.1081080E-13 * * * * -------------------------------------------------------------------- S E C O N D T E S T W I T H D I R E C T - A C C E S S I / O -------------------------------------------------------------------- * * * * FIRST DIMENSION OF A: M = 21 SECOND DIMENSION OF A: NA = 511 * NUMBER OF EQUATIONS: N = 1000 NUMBER OF ELEMENTS IN THE BAND OF THE MATRIX A: M*N = 21000 * * * MEMORY REQUIREMENTS: * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A AND THE RHS B: 2*(M*NA+N) = 23462 * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A: 2*(M*NA) = 21462 * * * * ALOGDT = 0.1211436E+04 SIGNDT = 1. * PIVMAX = 0.3727772E+02 PIVMIN = 0.9667904E+01 IPIV = 994 * INFO = 0 * JERR = 0 * * ERROR IN THE L2-NORM: 0.1497275E-14 AND IN THE MAXIMUM-NORM: 0.1081080E-13 * * * * ----------------------------------------------------------------------- T H I R D T E S T W I T H T O O S M A L L V A L U E OF NA ----------------------------------------------------------------------- * * * * FIRST DIMENSION OF A: M = 21 SECOND DIMENSION OF A: NA = 7 * NUMBER OF EQUATIONS: N = 1000 NUMBER OF ELEMENTS IN THE BAND OF THE MATRIX A: M*N = 21000 * * * * * JERR = 1 * RETURN TO CALLING PROGRAM WITH ERROR CODE "JERR = 1," BECAUSE NA = 7 IS LESS THAN (M+1)/2 + 2 = 13. * * ----------------------------------------------------------------------- F O U R T H T E S T W I T H S M A L L E S T P I V O T R E P L ----------------------------------------------------------------------- * * * * FIRST DIMENSION OF A: M = 21 SECOND DIMENSION OF A: NA = 511 * NUMBER OF EQUATIONS: N = 1000 NUMBER OF ELEMENTS IN THE BAND OF THE MATRIX A: M*N = 21000 * * * MEMORY REQUIREMENTS: * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A AND THE RHS B: 2*(M*NA+N) = 23462 * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A: 2*(M*NA) = 21462 * * * * ALOGDT = 0.1211436E+04 SIGNDT = 1. * PIVMAX = 0.3727772E+02 PIVMIN = 0.9667910E+01 IPIV = 994 * INFO = 994 * JERR = 0 * * ERROR IN THE L2-NORM: 0.4620833E-07 AND IN THE MAXIMUM-NORM: 0.6845441E-06 * * * * --------------------------------------------------------------------- F I F T H T E S T W I T H S M A L L P I V O T E L E M E N T --------------------------------------------------------------------- * * * * FIRST DIMENSION OF A: M = 21 SECOND DIMENSION OF A: NA = 241 * NUMBER OF EQUATIONS: N = 750 NUMBER OF ELEMENTS IN THE BAND OF THE MATRIX A: M*N = 15750 * * * MEMORY REQUIREMENTS: * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A AND THE RHS B: 2*(M*NA+N) = 11622 * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A: 2*(M*NA) = 10122 * * DEL = 0.1000000E-11 * EPS = 0.1000000E-11 * * ALOGDT = 0.8937705E+03 SIGNDT = -1. * PIVMAX = 0.3542299E+02 PIVMIN = 0.2283666E-09 IPIV = 750 * INFO = 0 * JERR = 0 * * ERROR IN THE L2-NORM: 0.2320288E-12 AND IN THE MAXIMUM-NORM: 0.3152451E-11 * * * * ----------------------------------------------------------------------- S I X T H T E S T W I T H S M A L L E S T P I V O T R E P L A ----------------------------------------------------------------------- * * * * FIRST DIMENSION OF A: M = 21 SECOND DIMENSION OF A: NA = 241 * NUMBER OF EQUATIONS: N = 750 NUMBER OF ELEMENTS IN THE BAND OF THE MATRIX A: M*N = 15750 * * * MEMORY REQUIREMENTS: * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A AND THE RHS B: 2*(M*NA+N) = 11622 * NUMBER OF FORTRAN STORAGE UNITS NEEDED TO STORE THE MATRIX A: 2*(M*NA) = 10122 * * DEL = 0.1000000E-11 * EPS = 0.2283700E-09 * * ALOGDT = 0.8937705E+03 SIGNDT = -1. * PIVMAX = 0.3542299E+02 PIVMIN = 0.2283700E-09 IPIV = 750 * INFO = 750 * JERR = 0 * * ERROR IN THE L2-NORM: 0.2057232E-02 AND IN THE MAXIMUM-NORM: 0.2803182E-01 * *