SUBROUTINE SGEFA(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(*),INFO REAL A(LDA,*) C C SGEFA FACTORS A REAL MATRIX BY GAUSSIAN ELIMINATION. C C SGEFA IS USUALLY CALLED BY SGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR SGECO) = (1 + 9/N)*(TIME FOR SGEFA) . C C ON ENTRY C C A REAL(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 SGESL OR SGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN SGECO 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 SAXPY,SSCAL,ISAMAX C C INTERNAL VARIABLES C REAL T INTEGER ISAMAX,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 = ISAMAX(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.0E0) 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.0E0/A(K,K) CALL SSCAL(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 SAXPY(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.0E0) INFO = N RETURN END SUBROUTINE SPOFA(A,LDA,N,INFO) INTEGER LDA,N,INFO REAL A(LDA,*) C C SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX. C C SPOFA IS USUALLY CALLED BY SPOCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR SPOCO) = (1 + 18/N)*(TIME FOR SPOFA) . C C ON ENTRY C C A REAL(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 SDOT C FORTRAN SQRT C C INTERNAL VARIABLES C REAL SDOT,T REAL S INTEGER J,JM1,K C BEGIN BLOCK WITH ...EXITS TO 40 C C DO 30 J = 1, N INFO = J S = 0.0E0 JM1 = J - 1 IF (JM1 .LT. 1) GO TO 20 DO 10 K = 1, JM1 T = A(K,J) - SDOT(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.0E0) GO TO 40 A(J,J) = SQRT(S) 30 CONTINUE INFO = 0 40 CONTINUE RETURN END SUBROUTINE SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB) INTEGER LDX,N,P,JOB INTEGER JPVT(*) REAL X(LDX,*),QRAUX(*),WORK(*) C C SQRDC 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 REAL(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 REAL(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 REAL(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 SQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C BLAS SAXPY,SDOT,SSCAL,SSWAP,SNRM2 C FORTRAN ABS,AMAX1,MIN0,SQRT C C INTERNAL VARIABLES C INTEGER J,JJ,JP,L,LP1,LUP,MAXJ,PL,PU REAL MAXNRM,SNRM2,TT REAL SDOT,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 SSWAP(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 SSWAP(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) = SNRM2(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.0E0 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 SSWAP(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.0E0 IF (L .EQ. N) GO TO 190 C C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L. C NRMXL = SNRM2(N-L+1,X(L,L),1) IF (NRMXL .EQ. 0.0E0) GO TO 180 IF (X(L,L) .NE. 0.0E0) NRMXL = SIGN(NRMXL,X(L,L)) CALL SSCAL(N-L+1,1.0E0/NRMXL,X(L,L),1) X(L,L) = 1.0E0 + 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 = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL SAXPY(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.0E0) GO TO 150 TT = 1.0E0 - (ABS(X(L,J))/QRAUX(J))**2 TT = AMAX1(TT,0.0E0) T = TT TT = 1.0E0 + 0.05E0*TT*(QRAUX(J)/WORK(J))**2 IF (TT .EQ. 1.0E0) GO TO 130 QRAUX(J) = QRAUX(J)*SQRT(T) GO TO 140 130 CONTINUE QRAUX(J) = SNRM2(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 SGTSL(N,C,D,E,B,INFO) INTEGER N,INFO REAL C(*),D(*),E(*),B(*) C C SGTSL 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 REAL(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 REAL(N) C IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX. C ON OUTPUT D IS DESTROYED. C C E REAL(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 REAL(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 ABS C C INTERNAL VARIABLES C INTEGER K,KB,KP1,NM1,NM2 REAL 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.0E0 E(N) = 0.0E0 C DO 30 K = 1, NM1 KP1 = K + 1 C C FIND THE LARGEST OF THE TWO ROWS C IF (ABS(C(KP1)) .LT. ABS(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.0E0) 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.0E0 B(KP1) = B(KP1) + T*B(K) 30 CONTINUE 40 CONTINUE IF (C(N) .NE. 0.0E0) 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 SPTSL(N,D,E,B) INTEGER N REAL D(*),E(*),B(*) C C SPTSL 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 REAL(N) C IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX. C ON OUTPUT D IS DESTROYED. C C E REAL(N) C IS THE OFFDIAGONAL OF THE TRIDIAGONAL MATRIX. C E(1) THROUGH E(N-1) SHOULD CONTAIN THE C OFFDIAGONAL. C C B REAL(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 REAL 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