C PCON61SUBD.F 25 February 1991 C INTEGER FUNCTION IDAMAX(N,DX,INCX) C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DMAX INTEGER I,INCX,IX,N C IDAMAX = 0 IF( N .LT. 1 ) RETURN IDAMAX = 1 IF(N.EQ.1)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 DMAX = DABS(DX(1)) IX = IX + INCX DO 10 I = 2,N IF(DABS(DX(IX)).LE.DMAX) GO TO 5 IDAMAX = I DMAX = DABS(DX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DMAX = DABS(DX(1)) DO 30 I = 2,N IF(DABS(DX(I)).LE.DMAX) GO TO 30 IDAMAX = I DMAX = DABS(DX(I)) 30 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0.0D0 DTEMP = 0.0D0 IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT = DTEMP RETURN END SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(1),INFO DOUBLE PRECISION A(LDA,1) C C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION. C C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) . C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE C INDICATION OF SINGULARITY. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS DAXPY,DSCAL,IDAMAX C C INTERNAL VARIABLES C DOUBLE PRECISION T INTEGER IDAMAX,J,K,KP1,L,NM1 C C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C INFO = 0 NM1 = N - 1 IF (NM1 .LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = IDAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (A(L,K) .EQ. 0.0D0) GO TO 40 C C INTERCHANGE IF NECESSARY C IF (L .EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE C C COMPUTE MULTIPLIERS C T = -1.0D0/A(K,K) CALL DSCAL(N-K,T,A(K+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C DO 30 J = KP1, N T = A(L,J) IF (L .EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N) .EQ. 0.0D0) INFO = N RETURN END SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(1),JOB DOUBLE PRECISION A(LDA,1),B(1) C C DGESL SOLVES THE DOUBLE PRECISION SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY DGECO OR DGEFA. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGECO OR DGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM DGECO OR DGEFA. C C B DOUBLE PRECISION(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B WHERE C TRANS(A) IS THE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0 C OR DGEFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL DGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL DGESL(A,LDA,N,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS DAXPY,DDOT C C INTERNAL VARIABLES C DOUBLE PRECISION DDOT,T INTEGER K,KB,L,NM1 C NM1 = N - 1 IF (JOB .NE. 0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (NM1 .LT. 1) GO TO 30 DO 20 K = 1, NM1 L = IPVT(K) T = B(L) IF (L .EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL DAXPY(K-1,T,A(1,K),1,B(1),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1, N T = DDOT(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1) L = IPVT(K) IF (L .EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE DGEDI(A,LDA,N,IPVT,DET,WORK,JOB) INTEGER LDA,N,IPVT(1),JOB DOUBLE PRECISION A(LDA,1),DET(2),WORK(1) C C DGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX C USING THE FACTORS COMPUTED BY DGECO OR DGEFA. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGECO OR DGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM DGECO OR DGEFA. C C WORK DOUBLE PRECISION(N) C WORK VECTOR. CONTENTS DESTROYED. C C JOB INTEGER C = 11 BOTH DETERMINANT AND INVERSE. C = 01 INVERSE ONLY. C = 10 DETERMINANT ONLY. C C ON RETURN C C A INVERSE OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE UNCHANGED. C C DET DOUBLE PRECISION(2) C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE NOT REFERENCED. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DABS(DET(1)) .LT. 10.0 C OR DET(1) .EQ. 0.0 . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED. C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY C AND IF DGECO HAS SET RCOND .GT. 0.0 OR DGEFA HAS SET C INFO .EQ. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS DAXPY,DSCAL,DSWAP C FORTRAN DABS,MOD C C INTERNAL VARIABLES C DOUBLE PRECISION T DOUBLE PRECISION TEN INTEGER I,J,K,KB,KP1,L,NM1 C C C COMPUTE DETERMINANT C IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0D0 DET(2) = 0.0D0 TEN = 10.0D0 DO 50 I = 1, N IF (IPVT(I) .NE. I) DET(1) = -DET(1) DET(1) = A(I,I)*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (DABS(DET(1)) .GE. 1.0D0) GO TO 20 DET(1) = TEN*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (DABS(DET(1)) .LT. TEN) GO TO 40 DET(1) = DET(1)/TEN DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C C COMPUTE INVERSE(U) C IF (MOD(JOB,10) .EQ. 0) GO TO 150 DO 100 K = 1, N A(K,K) = 1.0D0/A(K,K) T = -A(K,K) CALL DSCAL(K-1,T,A(1,K),1) KP1 = K + 1 IF (N .LT. KP1) GO TO 90 DO 80 J = KP1, N T = A(K,J) A(K,J) = 0.0D0 CALL DAXPY(K,T,A(1,K),1,A(1,J),1) 80 CONTINUE 90 CONTINUE 100 CONTINUE C C FORM INVERSE(U)*INVERSE(L) C NM1 = N - 1 IF (NM1 .LT. 1) GO TO 140 DO 130 KB = 1, NM1 K = N - KB KP1 = K + 1 DO 110 I = KP1, N WORK(I) = A(I,K) A(I,K) = 0.0D0 110 CONTINUE DO 120 J = KP1, N T = WORK(J) CALL DAXPY(N,T,A(1,J),1,A(1,K),1) 120 CONTINUE L = IPVT(K) IF (L .NE. K) CALL DSWAP(N,A(1,K),1,A(1,L),1) 130 CONTINUE 140 CONTINUE 150 CONTINUE RETURN END SUBROUTINE DGBFA(ABD,LDA,N,ML,MU,IPVT,INFO) INTEGER LDA,N,ML,MU,IPVT(1),INFO DOUBLE PRECISION ABD(LDA,1) C C DGBFA FACTORS A DOUBLE PRECISION BAND MATRIX BY ELIMINATION. C C DGBFA IS USUALLY CALLED BY DGBCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C C ON ENTRY C C ABD DOUBLE PRECISION(LDA, N) C CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS C OF THE MATRIX ARE STORED IN THE COLUMNS OF ABD AND C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS C ML+1 THROUGH 2*ML+MU+1 OF ABD . C SEE THE COMMENTS BELOW FOR DETAILS. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C LDA MUST BE .GE. 2*ML + MU + 1 . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C 0 .LE. ML .LT. N . C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C 0 .LE. MU .LT. N . C MORE EFFICIENT IF ML .LE. MU . C ON RETURN C C ABD AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT DGBSL WILL DIVIDE BY ZERO IF C CALLED. USE RCOND IN DGBCO FOR A RELIABLE C INDICATION OF SINGULARITY. C C BAND STORAGE C C IF A IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT C WILL SET UP THE INPUT. C C ML = (BAND WIDTH BELOW THE DIAGONAL) C MU = (BAND WIDTH ABOVE THE DIAGONAL) C M = ML + MU + 1 C DO 20 J = 1, N C I1 = MAX0(1, J-MU) C I2 = MIN0(N, J+ML) C DO 10 I = I1, I2 C K = I - J + M C ABD(K,J) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C THIS USES ROWS ML+1 THROUGH 2*ML+MU+1 OF ABD . C IN ADDITION, THE FIRST ML ROWS IN ABD ARE USED FOR C ELEMENTS GENERATED DURING THE TRIANGULARIZATION. C THE TOTAL NUMBER OF ROWS NEEDED IN ABD IS 2*ML+MU+1 . C THE ML+MU BY ML+MU UPPER LEFT TRIANGLE AND THE C ML BY ML LOWER RIGHT TRIANGLE ARE NOT REFERENCED. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS DAXPY,DSCAL,IDAMAX C FORTRAN MAX0,MIN0 C C INTERNAL VARIABLES C DOUBLE PRECISION T INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1 C C M = ML + MU + 1 INFO = 0 C C ZERO INITIAL FILL-IN COLUMNS C J0 = MU + 2 J1 = MIN0(N,M) - 1 IF (J1 .LT. J0) GO TO 30 DO 20 JZ = J0, J1 I0 = M + 1 - JZ DO 10 I = I0, ML ABD(I,JZ) = 0.0D0 10 CONTINUE 20 CONTINUE 30 CONTINUE JZ = J1 JU = 0 C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C NM1 = N - 1 IF (NM1 .LT. 1) GO TO 130 DO 120 K = 1, NM1 KP1 = K + 1 C C ZERO NEXT FILL-IN COLUMN C JZ = JZ + 1 IF (JZ .GT. N) GO TO 50 IF (ML .LT. 1) GO TO 50 DO 40 I = 1, ML ABD(I,JZ) = 0.0D0 40 CONTINUE 50 CONTINUE C C FIND L = PIVOT INDEX C LM = MIN0(ML,N-K) L = IDAMAX(LM+1,ABD(M,K),1) + M - 1 IPVT(K) = L + K - M C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (ABD(L,K) .EQ. 0.0D0) GO TO 100 C C INTERCHANGE IF NECESSARY C IF (L .EQ. M) GO TO 60 T = ABD(L,K) ABD(L,K) = ABD(M,K) ABD(M,K) = T 60 CONTINUE C C COMPUTE MULTIPLIERS C T = -1.0D0/ABD(M,K) CALL DSCAL(LM,T,ABD(M+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C JU = MIN0(MAX0(JU,MU+IPVT(K)),N) MM = M IF (JU .LT. KP1) GO TO 90 DO 80 J = KP1, JU L = L - 1 MM = MM - 1 T = ABD(L,J) IF (L .EQ. MM) GO TO 70 ABD(L,J) = ABD(MM,J) ABD(MM,J) = T 70 CONTINUE CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1) 80 CONTINUE 90 CONTINUE GO TO 110 100 CONTINUE INFO = K 110 CONTINUE 120 CONTINUE 130 CONTINUE IPVT(N) = N IF (ABD(M,N) .EQ. 0.0D0) INFO = N RETURN END SUBROUTINE DGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB) INTEGER LDA,N,ML,MU,IPVT(1),JOB DOUBLE PRECISION ABD(LDA,1),B(1) C C DGBSL SOLVES THE DOUBLE PRECISION BAND SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY DGBCO OR DGBFA. C C ON ENTRY C C ABD DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGBCO OR DGBFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM DGBCO OR DGBFA. C C B DOUBLE PRECISION(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B , WHERE C TRANS(A) IS THE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF DGBCO HAS SET RCOND .GT. 0.0 C OR DGBFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS DAXPY,DDOT C FORTRAN MIN0 C C INTERNAL VARIABLES C DOUBLE PRECISION DDOT,T INTEGER K,KB,L,LA,LB,LM,M,NM1 C M = MU + ML + 1 NM1 = N - 1 IF (JOB .NE. 0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (ML .EQ. 0) GO TO 30 IF (NM1 .LT. 1) GO TO 30 DO 20 K = 1, NM1 LM = MIN0(ML,N-K) L = IPVT(K) T = B(L) IF (L .EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1, N K = N + 1 - KB B(K) = B(K)/ABD(M,K) LM = MIN0(K,M) - 1 LA = M - LM LB = K - LM T = -B(K) CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1, N LM = MIN0(K,M) - 1 LA = M - LM LB = K - LM T = DDOT(LM,ABD(LA,K),1,B(LB),1) B(K) = (B(K) - T)/ABD(M,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (ML .EQ. 0) GO TO 90 IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB LM = MIN0(ML,N-K) B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1) L = IPVT(K) IF (L .EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE DGBDI(ABD,LDA,N,ML,MU,IPVT,DET) INTEGER LDA,N,ML,MU,IPVT(1) DOUBLE PRECISION ABD(LDA,1),DET(2) C C DGBDI COMPUTES THE DETERMINANT OF A BAND MATRIX C USING THE FACTORS COMPUTED BY DGBCO OR DGBFA. C IF THE INVERSE IS NEEDED, USE DGBSL N TIMES. C C ON ENTRY C C ABD DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGBCO OR DGBFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM DGBCO OR DGBFA. C C ON RETURN C C DET DOUBLE PRECISION(2) C DETERMINANT OF ORIGINAL MATRIX. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DABS(DET(1)) .LT. 10.0 C OR DET(1) = 0.0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C FORTRAN DABS C C INTERNAL VARIABLES C DOUBLE PRECISION TEN INTEGER I,M C C M = ML + MU + 1 DET(1) = 1.0D0 DET(2) = 0.0D0 TEN = 10.0D0 DO 50 I = 1, N IF (IPVT(I) .NE. I) DET(1) = -DET(1) DET(1) = ABD(M,I)*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (DABS(DET(1)) .GE. 1.0D0) GO TO 20 DET(1) = TEN*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (DABS(DET(1)) .LT. TEN) GO TO 40 DET(1) = DET(1)/TEN DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE RETURN END SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DA INTEGER I,INCX,INCY,M,MP1,N C IF(N.LE.0)RETURN IF (DA .EQ. 0.0D0) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) INTEGER NEXT DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE DATA ZERO, ONE /0.0D0, 1.0D0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C IF(N .GT. 0) GO TO 10 DNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C C PREPARE FOR PHASE 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 85 HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J =I,NN,INCX IF(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END SUBROUTINE DSCAL(N,DA,DX,INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DA,DX(1) INTEGER I,INCX,M,MP1,N,NINCX C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END SUBROUTINE DSWAP (N,DX,INCX,DY,INCY) C C INTERCHANGES TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL C TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,3) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP 30 CONTINUE IF( N .LT. 3 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,3 DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP DTEMP = DX(I + 1) DX(I + 1) = DY(I + 1) DY(I + 1) = DTEMP DTEMP = DX(I + 2) DX(I + 2) = DY(I + 2) DY(I + 2) = DTEMP 50 CONTINUE RETURN END