593 lines
17 KiB
FortranFixed
593 lines
17 KiB
FortranFixed
SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO)
|
|
INTEGER LDA,N,IPVT(*),INFO
|
|
DOUBLE PRECISION A(LDA,*)
|
|
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 DPOFA(A,LDA,N,INFO)
|
|
INTEGER LDA,N,INFO
|
|
DOUBLE PRECISION A(LDA,*)
|
|
C
|
|
C DPOFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE
|
|
C MATRIX.
|
|
C
|
|
C DPOFA IS USUALLY CALLED BY DPOCO, BUT IT CAN BE CALLED
|
|
C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
|
|
C (TIME FOR DPOCO) = (1 + 18/N)*(TIME FOR DPOFA) .
|
|
C
|
|
C ON ENTRY
|
|
C
|
|
C A DOUBLE PRECISION(LDA, N)
|
|
C THE SYMMETRIC MATRIX TO BE FACTORED. ONLY THE
|
|
C DIAGONAL AND UPPER TRIANGLE ARE USED.
|
|
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 R SO THAT A = TRANS(R)*R
|
|
C WHERE TRANS(R) IS THE TRANSPOSE.
|
|
C THE STRICT LOWER TRIANGLE IS UNALTERED.
|
|
C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
|
|
C
|
|
C INFO INTEGER
|
|
C = 0 FOR NORMAL RETURN.
|
|
C = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR
|
|
C OF ORDER K IS NOT POSITIVE DEFINITE.
|
|
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 DDOT
|
|
C FORTRAN DSQRT
|
|
C
|
|
C INTERNAL VARIABLES
|
|
C
|
|
DOUBLE PRECISION DDOT,T
|
|
DOUBLE PRECISION S
|
|
INTEGER J,JM1,K
|
|
C BEGIN BLOCK WITH ...EXITS TO 40
|
|
C
|
|
C
|
|
DO 30 J = 1, N
|
|
INFO = J
|
|
S = 0.0D0
|
|
JM1 = J - 1
|
|
IF (JM1 .LT. 1) GO TO 20
|
|
DO 10 K = 1, JM1
|
|
T = A(K,J) - DDOT(K-1,A(1,K),1,A(1,J),1)
|
|
T = T/A(K,K)
|
|
A(K,J) = T
|
|
S = S + T*T
|
|
10 CONTINUE
|
|
20 CONTINUE
|
|
S = A(J,J) - S
|
|
C ......EXIT
|
|
IF (S .LE. 0.0D0) GO TO 40
|
|
A(J,J) = DSQRT(S)
|
|
30 CONTINUE
|
|
INFO = 0
|
|
40 CONTINUE
|
|
RETURN
|
|
END
|
|
SUBROUTINE DQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB)
|
|
INTEGER LDX,N,P,JOB
|
|
INTEGER JPVT(*)
|
|
DOUBLE PRECISION X(LDX,*),QRAUX(*),WORK(*)
|
|
C
|
|
C DQRDC USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR
|
|
C FACTORIZATION OF AN N BY P MATRIX X. COLUMN PIVOTING
|
|
C BASED ON THE 2-NORMS OF THE REDUCED COLUMNS MAY BE
|
|
C PERFORMED AT THE USERS OPTION.
|
|
C
|
|
C ON ENTRY
|
|
C
|
|
C X DOUBLE PRECISION(LDX,P), WHERE LDX .GE. N.
|
|
C X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE
|
|
C COMPUTED.
|
|
C
|
|
C LDX INTEGER.
|
|
C LDX IS THE LEADING DIMENSION OF THE ARRAY X.
|
|
C
|
|
C N INTEGER.
|
|
C N IS THE NUMBER OF ROWS OF THE MATRIX X.
|
|
C
|
|
C P INTEGER.
|
|
C P IS THE NUMBER OF COLUMNS OF THE MATRIX X.
|
|
C
|
|
C JPVT INTEGER(P).
|
|
C JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION
|
|
C OF THE PIVOT COLUMNS. THE K-TH COLUMN X(K) OF X
|
|
C IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE
|
|
C VALUE OF JPVT(K).
|
|
C
|
|
C IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL
|
|
C COLUMN.
|
|
C
|
|
C IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN.
|
|
C
|
|
C IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN.
|
|
C
|
|
C BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS
|
|
C ARE MOVED TO THE BEGINNING OF THE ARRAY X AND FINAL
|
|
C COLUMNS TO THE END. BOTH INITIAL AND FINAL COLUMNS
|
|
C ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY
|
|
C FREE COLUMNS ARE MOVED. AT THE K-TH STAGE OF THE
|
|
C REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN
|
|
C IT IS INTERCHANGED WITH THE FREE COLUMN OF LARGEST
|
|
C REDUCED NORM. JPVT IS NOT REFERENCED IF
|
|
C JOB .EQ. 0.
|
|
C
|
|
C WORK DOUBLE PRECISION(P).
|
|
C WORK IS A WORK ARRAY. WORK IS NOT REFERENCED IF
|
|
C JOB .EQ. 0.
|
|
C
|
|
C JOB INTEGER.
|
|
C JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING.
|
|
C IF JOB .EQ. 0, NO PIVOTING IS DONE.
|
|
C IF JOB .NE. 0, PIVOTING IS DONE.
|
|
C
|
|
C ON RETURN
|
|
C
|
|
C X X CONTAINS IN ITS UPPER TRIANGLE THE UPPER
|
|
C TRIANGULAR MATRIX R OF THE QR FACTORIZATION.
|
|
C BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM
|
|
C WHICH THE ORTHOGONAL PART OF THE DECOMPOSITION
|
|
C CAN BE RECOVERED. NOTE THAT IF PIVOTING HAS
|
|
C BEEN REQUESTED, THE DECOMPOSITION IS NOT THAT
|
|
C OF THE ORIGINAL MATRIX X BUT THAT OF X
|
|
C WITH ITS COLUMNS PERMUTED AS DESCRIBED BY JPVT.
|
|
C
|
|
C QRAUX DOUBLE PRECISION(P).
|
|
C QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER
|
|
C THE ORTHOGONAL PART OF THE DECOMPOSITION.
|
|
C
|
|
C JPVT JPVT(K) CONTAINS THE INDEX OF THE COLUMN OF THE
|
|
C ORIGINAL MATRIX THAT HAS BEEN INTERCHANGED INTO
|
|
C THE K-TH COLUMN, IF PIVOTING WAS REQUESTED.
|
|
C
|
|
C LINPACK. THIS VERSION DATED 08/14/78 .
|
|
C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
|
|
C
|
|
C DQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
|
|
C
|
|
C BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2
|
|
C FORTRAN DABS,DMAX1,MIN0,DSQRT
|
|
C
|
|
C INTERNAL VARIABLES
|
|
C
|
|
INTEGER J,JJ,JP,L,LP1,LUP,MAXJ,PL,PU
|
|
DOUBLE PRECISION MAXNRM,DNRM2,TT
|
|
DOUBLE PRECISION DDOT,NRMXL,T
|
|
LOGICAL NEGJ,SWAPJ
|
|
C
|
|
C
|
|
PL = 1
|
|
PU = 0
|
|
IF (JOB .EQ. 0) GO TO 60
|
|
C
|
|
C PIVOTING HAS BEEN REQUESTED. REARRANGE THE COLUMNS
|
|
C ACCORDING TO JPVT.
|
|
C
|
|
DO 20 J = 1, P
|
|
SWAPJ = JPVT(J) .GT. 0
|
|
NEGJ = JPVT(J) .LT. 0
|
|
JPVT(J) = J
|
|
IF (NEGJ) JPVT(J) = -J
|
|
IF (.NOT.SWAPJ) GO TO 10
|
|
IF (J .NE. PL) CALL DSWAP(N,X(1,PL),1,X(1,J),1)
|
|
JPVT(J) = JPVT(PL)
|
|
JPVT(PL) = J
|
|
PL = PL + 1
|
|
10 CONTINUE
|
|
20 CONTINUE
|
|
PU = P
|
|
DO 50 JJ = 1, P
|
|
J = P - JJ + 1
|
|
IF (JPVT(J) .GE. 0) GO TO 40
|
|
JPVT(J) = -JPVT(J)
|
|
IF (J .EQ. PU) GO TO 30
|
|
CALL DSWAP(N,X(1,PU),1,X(1,J),1)
|
|
JP = JPVT(PU)
|
|
JPVT(PU) = JPVT(J)
|
|
JPVT(J) = JP
|
|
30 CONTINUE
|
|
PU = PU - 1
|
|
40 CONTINUE
|
|
50 CONTINUE
|
|
60 CONTINUE
|
|
C
|
|
C COMPUTE THE NORMS OF THE FREE COLUMNS.
|
|
C
|
|
IF (PU .LT. PL) GO TO 80
|
|
DO 70 J = PL, PU
|
|
QRAUX(J) = DNRM2(N,X(1,J),1)
|
|
WORK(J) = QRAUX(J)
|
|
70 CONTINUE
|
|
80 CONTINUE
|
|
C
|
|
C PERFORM THE HOUSEHOLDER REDUCTION OF X.
|
|
C
|
|
LUP = MIN0(N,P)
|
|
DO 200 L = 1, LUP
|
|
IF (L .LT. PL .OR. L .GE. PU) GO TO 120
|
|
C
|
|
C LOCATE THE COLUMN OF LARGEST NORM AND BRING IT
|
|
C INTO THE PIVOT POSITION.
|
|
C
|
|
MAXNRM = 0.0D0
|
|
MAXJ = L
|
|
DO 100 J = L, PU
|
|
IF (QRAUX(J) .LE. MAXNRM) GO TO 90
|
|
MAXNRM = QRAUX(J)
|
|
MAXJ = J
|
|
90 CONTINUE
|
|
100 CONTINUE
|
|
IF (MAXJ .EQ. L) GO TO 110
|
|
CALL DSWAP(N,X(1,L),1,X(1,MAXJ),1)
|
|
QRAUX(MAXJ) = QRAUX(L)
|
|
WORK(MAXJ) = WORK(L)
|
|
JP = JPVT(MAXJ)
|
|
JPVT(MAXJ) = JPVT(L)
|
|
JPVT(L) = JP
|
|
110 CONTINUE
|
|
120 CONTINUE
|
|
QRAUX(L) = 0.0D0
|
|
IF (L .EQ. N) GO TO 190
|
|
C
|
|
C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L.
|
|
C
|
|
NRMXL = DNRM2(N-L+1,X(L,L),1)
|
|
IF (NRMXL .EQ. 0.0D0) GO TO 180
|
|
IF (X(L,L) .NE. 0.0D0) NRMXL = DSIGN(NRMXL,X(L,L))
|
|
CALL DSCAL(N-L+1,1.0D0/NRMXL,X(L,L),1)
|
|
X(L,L) = 1.0D0 + X(L,L)
|
|
C
|
|
C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS,
|
|
C UPDATING THE NORMS.
|
|
C
|
|
LP1 = L + 1
|
|
IF (P .LT. LP1) GO TO 170
|
|
DO 160 J = LP1, P
|
|
T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
|
|
CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
|
|
IF (J .LT. PL .OR. J .GT. PU) GO TO 150
|
|
IF (QRAUX(J) .EQ. 0.0D0) GO TO 150
|
|
TT = 1.0D0 - (DABS(X(L,J))/QRAUX(J))**2
|
|
TT = DMAX1(TT,0.0D0)
|
|
T = TT
|
|
TT = 1.0D0 + 0.05D0*TT*(QRAUX(J)/WORK(J))**2
|
|
IF (TT .EQ. 1.0D0) GO TO 130
|
|
QRAUX(J) = QRAUX(J)*DSQRT(T)
|
|
GO TO 140
|
|
130 CONTINUE
|
|
QRAUX(J) = DNRM2(N-L,X(L+1,J),1)
|
|
WORK(J) = QRAUX(J)
|
|
140 CONTINUE
|
|
150 CONTINUE
|
|
160 CONTINUE
|
|
170 CONTINUE
|
|
C
|
|
C SAVE THE TRANSFORMATION.
|
|
C
|
|
QRAUX(L) = X(L,L)
|
|
X(L,L) = -NRMXL
|
|
180 CONTINUE
|
|
190 CONTINUE
|
|
200 CONTINUE
|
|
RETURN
|
|
END
|
|
SUBROUTINE DGTSL(N,C,D,E,B,INFO)
|
|
INTEGER N,INFO
|
|
DOUBLE PRECISION C(*),D(*),E(*),B(*)
|
|
C
|
|
C DGTSL GIVEN A GENERAL TRIDIAGONAL MATRIX AND A RIGHT HAND
|
|
C SIDE WILL FIND THE SOLUTION.
|
|
C
|
|
C ON ENTRY
|
|
C
|
|
C N INTEGER
|
|
C IS THE ORDER OF THE TRIDIAGONAL MATRIX.
|
|
C
|
|
C C DOUBLE PRECISION(N)
|
|
C IS THE SUBDIAGONAL OF THE TRIDIAGONAL MATRIX.
|
|
C C(2) THROUGH C(N) SHOULD CONTAIN THE SUBDIAGONAL.
|
|
C ON OUTPUT C IS DESTROYED.
|
|
C
|
|
C D DOUBLE PRECISION(N)
|
|
C IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX.
|
|
C ON OUTPUT D IS DESTROYED.
|
|
C
|
|
C E DOUBLE PRECISION(N)
|
|
C IS THE SUPERDIAGONAL OF THE TRIDIAGONAL MATRIX.
|
|
C E(1) THROUGH E(N-1) SHOULD CONTAIN THE SUPERDIAGONAL.
|
|
C ON OUTPUT E IS DESTROYED.
|
|
C
|
|
C B DOUBLE PRECISION(N)
|
|
C IS THE RIGHT HAND SIDE VECTOR.
|
|
C
|
|
C ON RETURN
|
|
C
|
|
C B IS THE SOLUTION VECTOR.
|
|
C
|
|
C INFO INTEGER
|
|
C = 0 NORMAL VALUE.
|
|
C = K IF THE K-TH ELEMENT OF THE DIAGONAL BECOMES
|
|
C EXACTLY ZERO. THE SUBROUTINE RETURNS WHEN
|
|
C THIS IS DETECTED.
|
|
C
|
|
C LINPACK. THIS VERSION DATED 08/14/78 .
|
|
C JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
|
|
C
|
|
C NO EXTERNALS
|
|
C FORTRAN DABS
|
|
C
|
|
C INTERNAL VARIABLES
|
|
C
|
|
INTEGER K,KB,KP1,NM1,NM2
|
|
DOUBLE PRECISION T
|
|
C BEGIN BLOCK PERMITTING ...EXITS TO 100
|
|
C
|
|
INFO = 0
|
|
C(1) = D(1)
|
|
NM1 = N - 1
|
|
IF (NM1 .LT. 1) GO TO 40
|
|
D(1) = E(1)
|
|
E(1) = 0.0D0
|
|
E(N) = 0.0D0
|
|
C
|
|
DO 30 K = 1, NM1
|
|
KP1 = K + 1
|
|
C
|
|
C FIND THE LARGEST OF THE TWO ROWS
|
|
C
|
|
IF (DABS(C(KP1)) .LT. DABS(C(K))) GO TO 10
|
|
C
|
|
C INTERCHANGE ROW
|
|
C
|
|
T = C(KP1)
|
|
C(KP1) = C(K)
|
|
C(K) = T
|
|
T = D(KP1)
|
|
D(KP1) = D(K)
|
|
D(K) = T
|
|
T = E(KP1)
|
|
E(KP1) = E(K)
|
|
E(K) = T
|
|
T = B(KP1)
|
|
B(KP1) = B(K)
|
|
B(K) = T
|
|
10 CONTINUE
|
|
C
|
|
C ZERO ELEMENTS
|
|
C
|
|
IF (C(K) .NE. 0.0D0) GO TO 20
|
|
INFO = K
|
|
C ............EXIT
|
|
GO TO 100
|
|
20 CONTINUE
|
|
T = -C(KP1)/C(K)
|
|
C(KP1) = D(KP1) + T*D(K)
|
|
D(KP1) = E(KP1) + T*E(K)
|
|
E(KP1) = 0.0D0
|
|
B(KP1) = B(KP1) + T*B(K)
|
|
30 CONTINUE
|
|
40 CONTINUE
|
|
IF (C(N) .NE. 0.0D0) GO TO 50
|
|
INFO = N
|
|
GO TO 90
|
|
50 CONTINUE
|
|
C
|
|
C BACK SOLVE
|
|
C
|
|
NM2 = N - 2
|
|
B(N) = B(N)/C(N)
|
|
IF (N .EQ. 1) GO TO 80
|
|
B(NM1) = (B(NM1) - D(NM1)*B(N))/C(NM1)
|
|
IF (NM2 .LT. 1) GO TO 70
|
|
DO 60 KB = 1, NM2
|
|
K = NM2 - KB + 1
|
|
B(K) = (B(K) - D(K)*B(K+1) - E(K)*B(K+2))/C(K)
|
|
60 CONTINUE
|
|
70 CONTINUE
|
|
80 CONTINUE
|
|
90 CONTINUE
|
|
100 CONTINUE
|
|
C
|
|
RETURN
|
|
END
|
|
SUBROUTINE DPTSL(N,D,E,B)
|
|
INTEGER N
|
|
DOUBLE PRECISION D(*),E(*),B(*)
|
|
C
|
|
C DPTSL GIVEN A POSITIVE DEFINITE TRIDIAGONAL MATRIX AND A RIGHT
|
|
C HAND SIDE WILL FIND THE SOLUTION.
|
|
C
|
|
C ON ENTRY
|
|
C
|
|
C N INTEGER
|
|
C IS THE ORDER OF THE TRIDIAGONAL MATRIX.
|
|
C
|
|
C D DOUBLE PRECISION(N)
|
|
C IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX.
|
|
C ON OUTPUT D IS DESTROYED.
|
|
C
|
|
C E DOUBLE PRECISION(N)
|
|
C IS THE OFFDIAGONAL OF THE TRIDIAGONAL MATRIX.
|
|
C E(1) THROUGH E(N-1) SHOULD CONTAIN THE
|
|
C OFFDIAGONAL.
|
|
C
|
|
C B DOUBLE PRECISION(N)
|
|
C IS THE RIGHT HAND SIDE VECTOR.
|
|
C
|
|
C ON RETURN
|
|
C
|
|
C B CONTAINS THE SOULTION.
|
|
C
|
|
C LINPACK. THIS VERSION DATED 08/14/78 .
|
|
C JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
|
|
C
|
|
C NO EXTERNALS
|
|
C FORTRAN MOD
|
|
C
|
|
C INTERNAL VARIABLES
|
|
C
|
|
INTEGER K,KBM1,KE,KF,KP1,NM1,NM1D2
|
|
DOUBLE PRECISION T1,T2
|
|
C
|
|
C CHECK FOR 1 X 1 CASE
|
|
C
|
|
IF (N .NE. 1) GO TO 10
|
|
B(1) = B(1)/D(1)
|
|
GO TO 70
|
|
10 CONTINUE
|
|
NM1 = N - 1
|
|
NM1D2 = NM1/2
|
|
IF (N .EQ. 2) GO TO 30
|
|
KBM1 = N - 1
|
|
C
|
|
C ZERO TOP HALF OF SUBDIAGONAL AND BOTTOM HALF OF
|
|
C SUPERDIAGONAL
|
|
C
|
|
DO 20 K = 1, NM1D2
|
|
T1 = E(K)/D(K)
|
|
D(K+1) = D(K+1) - T1*E(K)
|
|
B(K+1) = B(K+1) - T1*B(K)
|
|
T2 = E(KBM1)/D(KBM1+1)
|
|
D(KBM1) = D(KBM1) - T2*E(KBM1)
|
|
B(KBM1) = B(KBM1) - T2*B(KBM1+1)
|
|
KBM1 = KBM1 - 1
|
|
20 CONTINUE
|
|
30 CONTINUE
|
|
KP1 = NM1D2 + 1
|
|
C
|
|
C CLEAN UP FOR POSSIBLE 2 X 2 BLOCK AT CENTER
|
|
C
|
|
IF (MOD(N,2) .NE. 0) GO TO 40
|
|
T1 = E(KP1)/D(KP1)
|
|
D(KP1+1) = D(KP1+1) - T1*E(KP1)
|
|
B(KP1+1) = B(KP1+1) - T1*B(KP1)
|
|
KP1 = KP1 + 1
|
|
40 CONTINUE
|
|
C
|
|
C BACK SOLVE STARTING AT THE CENTER, GOING TOWARDS THE TOP
|
|
C AND BOTTOM
|
|
C
|
|
B(KP1) = B(KP1)/D(KP1)
|
|
IF (N .EQ. 2) GO TO 60
|
|
K = KP1 - 1
|
|
KE = KP1 + NM1D2 - 1
|
|
DO 50 KF = KP1, KE
|
|
B(K) = (B(K) - E(K)*B(K+1))/D(K)
|
|
B(KF+1) = (B(KF+1) - E(KF)*B(KF))/D(KF+1)
|
|
K = K - 1
|
|
50 CONTINUE
|
|
60 CONTINUE
|
|
IF (MOD(N,2) .EQ. 0) B(1) = (B(1) - E(1)*B(2))/D(1)
|
|
70 CONTINUE
|
|
RETURN
|
|
END
|