page@swan.ulowell.edu (Bob Page) (11/03/88)
Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink)
Posting-number: Volume 2, Issue 45
Archive-name: applications/matlab/src.5
# This is a shell archive.
# Remove everything above and including the cut line.
# Then run the rest of the file through sh.
#----cut here-----cut here-----cut here-----cut here----#
#!/bin/sh
# shar: Shell Archiver
# Run the following text with /bin/sh to create:
# src-5
# This archive created: Wed Nov 2 16:20:35 1988
cat << \SHAR_EOF > src-5
IF (N .LT. KP1) GO TO 100
DO 90 J = KP1, N
TR = AR(K,J)
TI = AI(K,J)
AR(K,J) = 0.0D0
AI(K,J) = 0.0D0
CALL WAXPY(K,TR,TI,AR(1,K),AI(1,K),1,AR(1,J),AI(1,J),1)
90 CONTINUE
100 CONTINUE
110 CONTINUE
C
C FORM INVERSE(U)*INVERSE(L)
C
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 150
DO 140 KB = 1, NM1
K = N - KB
KP1 = K + 1
DO 120 I = KP1, N
WORKR(I) = AR(I,K)
WORKI(I) = AI(I,K)
AR(I,K) = 0.0D0
AI(I,K) = 0.0D0
120 CONTINUE
DO 130 J = KP1, N
TR = WORKR(J)
TI = WORKI(J)
CALL WAXPY(N,TR,TI,AR(1,J),AI(1,J),1,AR(1,K),AI(1,K),1)
130 CONTINUE
L = IPVT(K)
IF (L .NE. K)
* CALL WSWAP(N,AR(1,K),AI(1,K),1,AR(1,L),AI(1,L),1)
140 CONTINUE
150 CONTINUE
160 CONTINUE
RETURN
END
SUBROUTINE WPOFA(AR,AI,LDA,N,INFO)
DOUBLE PRECISION AR(LDA,1),AI(LDA,1)
DOUBLE PRECISION S,TR,TI,WDOTCR,WDOTCI
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
TR = AR(K,J)-WDOTCR(K-1,AR(1,K),AI(1,K),1,AR(1,J),AI(1,J),1)
TI = AI(K,J)-WDOTCI(K-1,AR(1,K),AI(1,K),1,AR(1,J),AI(1,J),1)
CALL WDIV(TR,TI,AR(K,K),AI(K,K),TR,TI)
AR(K,J) = TR
AI(K,J) = TI
S = S + TR*TR + TI*TI
10 CONTINUE
20 CONTINUE
S = AR(J,J) - S
IF (S.LE.0.0D0 .OR. AI(J,J).NE.0.0D0) GO TO 40
AR(J,J) = DSQRT(S)
30 CONTINUE
INFO = 0
40 RETURN
END
SUBROUTINE RREF(AR,AI,LDA,M,N,EPS)
DOUBLE PRECISION AR(LDA,1),AI(LDA,1),EPS,TOL,TR,TI,WASUM
TOL = 0.0D0
DO 10 J = 1, N
TOL = DMAX1(TOL,WASUM(M,AR(1,J),AI(1,J),1))
10 CONTINUE
TOL = EPS*DFLOAT(2*MAX0(M,N))*TOL
K = 1
L = 1
20 IF (K.GT.M .OR. L.GT.N) RETURN
I = IWAMAX(M-K+1,AR(K,L),AI(K,L),1) + K-1
IF (DABS(AR(I,L))+DABS(AI(I,L)) .GT. TOL) GO TO 30
CALL WSET(M-K+1,0.0D0,0.0D0,AR(K,L),AI(K,L),1)
L = L+1
GO TO 20
30 CALL WSWAP(N-L+1,AR(I,L),AI(I,L),LDA,AR(K,L),AI(K,L),LDA)
CALL WDIV(1.0D0,0.0D0,AR(K,L),AI(K,L),TR,TI)
CALL WSCAL(N-L+1,TR,TI,AR(K,L),AI(K,L),LDA)
AR(K,L) = 1.0D0
AI(K,L) = 0.0D0
DO 40 I = 1, M
TR = -AR(I,L)
TI = -AI(I,L)
IF (I .NE. K) CALL WAXPY(N-L+1,TR,TI,
$ AR(K,L),AI(K,L),LDA,AR(I,L),AI(I,L),LDA)
40 CONTINUE
K = K+1
L = L+1
GO TO 20
END
SUBROUTINE HILBER(A,LDA,N)
DOUBLE PRECISION A(LDA,N)
C GENERATE INVERSE HILBERT MATRIX
DOUBLE PRECISION P,R
P = DFLOAT(N)
DO 20 I = 1, N
IF (I.NE.1) P = (DFLOAT(N-I+1)*P*DFLOAT(N+I-1))/DFLOAT(I-1)**2
R = P*P
A(I,I) = R/DFLOAT(2*I-1)
IF (I.EQ.N) GO TO 20
IP1 = I+1
DO 10 J = IP1, N
R = -(DFLOAT(N-J+1)*R*(N+J-1))/DFLOAT(J-1)**2
A(I,J) = R/DFLOAT(I+J-1)
A(J,I) = A(I,J)
10 CONTINUE
20 CONTINUE
RETURN
END
SUBROUTINE HTRIDI(NM,N,AR,AI,D,E,E2,TAU)
C
INTEGER I,J,K,L,N,II,NM,JP1
DOUBLE PRECISION AR(NM,N),AI(NM,N),D(N),E(N),E2(N),TAU(2,N)
DOUBLE PRECISION F,G,H,FI,GI,HH,SI,SCALE
DOUBLE PRECISION FLOP,PYTHAG
C
C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
C THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968)
C BY MARTIN, REINSCH, AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX
C TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING
C UNITARY SIMILARITY TRANSFORMATIONS.
C
C ON INPUT.
C
C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C DIMENSION STATEMENT.
C
C N IS THE ORDER OF THE MATRIX.
C
C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX.
C ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C ON OUTPUT.
C
C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
C FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER
C TRIANGLES. THEIR STRICT UPPER TRIANGLES AND THE
C DIAGONAL OF AR ARE UNALTERED.
C
C D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX.
C
C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO.
C
C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
C
C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
C
C MODIFIED TO GET RID OF ALL COMPLEX ARITHMETIC, C. MOLER, 6/27/79.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
TAU(1,N) = 1.0D0
TAU(2,N) = 0.0D0
C
DO 100 I = 1, N
100 D(I) = AR(I,I)
C .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
DO 300 II = 1, N
I = N + 1 - II
L = I - 1
H = 0.0D0
SCALE = 0.0D0
IF (L .LT. 1) GO TO 130
C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
DO 120 K = 1, L
120 SCALE = FLOP(SCALE + DABS(AR(I,K)) + DABS(AI(I,K)))
C
IF (SCALE .NE. 0.0D0) GO TO 140
TAU(1,L) = 1.0D0
TAU(2,L) = 0.0D0
130 E(I) = 0.0D0
E2(I) = 0.0D0
GO TO 290
C
140 DO 150 K = 1, L
AR(I,K) = FLOP(AR(I,K)/SCALE)
AI(I,K) = FLOP(AI(I,K)/SCALE)
H = FLOP(H + AR(I,K)*AR(I,K) + AI(I,K)*AI(I,K))
150 CONTINUE
C
E2(I) = FLOP(SCALE*SCALE*H)
G = FLOP(DSQRT(H))
E(I) = FLOP(SCALE*G)
F = PYTHAG(AR(I,L),AI(I,L))
C .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T ..........
IF (F .EQ. 0.0D0) GO TO 160
TAU(1,L) = FLOP((AI(I,L)*TAU(2,I) - AR(I,L)*TAU(1,I))/F)
SI = FLOP((AR(I,L)*TAU(2,I) + AI(I,L)*TAU(1,I))/F)
H = FLOP(H + F*G)
G = FLOP(1.0D0 + G/F)
AR(I,L) = FLOP(G*AR(I,L))
AI(I,L) = FLOP(G*AI(I,L))
IF (L .EQ. 1) GO TO 270
GO TO 170
160 TAU(1,L) = -TAU(1,I)
SI = TAU(2,I)
AR(I,L) = G
170 F = 0.0D0
C
DO 240 J = 1, L
G = 0.0D0
GI = 0.0D0
C .......... FORM ELEMENT OF A*U ..........
DO 180 K = 1, J
G = FLOP(G + AR(J,K)*AR(I,K) + AI(J,K)*AI(I,K))
GI = FLOP(GI - AR(J,K)*AI(I,K) + AI(J,K)*AR(I,K))
180 CONTINUE
C
JP1 = J + 1
IF (L .LT. JP1) GO TO 220
C
DO 200 K = JP1, L
G = FLOP(G + AR(K,J)*AR(I,K) - AI(K,J)*AI(I,K))
GI = FLOP(GI - AR(K,J)*AI(I,K) - AI(K,J)*AR(I,K))
200 CONTINUE
C .......... FORM ELEMENT OF P ..........
220 E(J) = FLOP(G/H)
TAU(2,J) = FLOP(GI/H)
F = FLOP(F + E(J)*AR(I,J) - TAU(2,J)*AI(I,J))
240 CONTINUE
C
HH = FLOP(F/(H + H))
C .......... FORM REDUCED A ..........
DO 260 J = 1, L
F = AR(I,J)
G = FLOP(E(J) - HH*F)
E(J) = G
FI = -AI(I,J)
GI = FLOP(TAU(2,J) - HH*FI)
TAU(2,J) = -GI
C
DO 260 K = 1, J
AR(J,K) = FLOP(AR(J,K) - F*E(K) - G*AR(I,K)
X + FI*TAU(2,K) + GI*AI(I,K))
AI(J,K) = FLOP(AI(J,K) - F*TAU(2,K) - G*AI(I,K)
X - FI*E(K) - GI*AR(I,K))
260 CONTINUE
C
270 DO 280 K = 1, L
AR(I,K) = FLOP(SCALE*AR(I,K))
AI(I,K) = FLOP(SCALE*AI(I,K))
280 CONTINUE
C
TAU(2,L) = -SI
290 HH = D(I)
D(I) = AR(I,I)
AR(I,I) = HH
AI(I,I) = FLOP(SCALE*DSQRT(H))
300 CONTINUE
C
RETURN
END
SUBROUTINE HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)
C
INTEGER I,J,K,L,M,N,NM
DOUBLE PRECISION AR(NM,N),AI(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M)
DOUBLE PRECISION H,S,SI,FLOP
C
C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
C THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968)
C BY MARTIN, REINSCH, AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN
C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
C REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY HTRIDI.
C
C ON INPUT.
C
C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C DIMENSION STATEMENT.
C
C N IS THE ORDER OF THE MATRIX.
C
C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
C FORMATIONS USED IN THE REDUCTION BY HTRIDI IN THEIR
C FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR.
C
C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
C
C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
C
C ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
C IN ITS FIRST M COLUMNS.
C
C ON OUTPUT.
C
C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
C IN THEIR FIRST M COLUMNS.
C
C NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR
C IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
IF (M .EQ. 0) GO TO 200
C .......... TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC
C TRIDIAGONAL MATRIX TO THOSE OF THE HERMITIAN
C TRIDIAGONAL MATRIX. ..........
DO 50 K = 1, N
C
DO 50 J = 1, M
ZI(K,J) = FLOP(-ZR(K,J)*TAU(2,K))
ZR(K,J) = FLOP(ZR(K,J)*TAU(1,K))
50 CONTINUE
C
IF (N .EQ. 1) GO TO 200
C .......... RECOVER AND APPLY THE HOUSEHOLDER MATRICES ..........
DO 140 I = 2, N
L = I - 1
H = AI(I,I)
IF (H .EQ. 0.0D0) GO TO 140
C
DO 130 J = 1, M
S = 0.0D0
SI = 0.0D0
C
DO 110 K = 1, L
S = FLOP(S + AR(I,K)*ZR(K,J) - AI(I,K)*ZI(K,J))
SI = FLOP(SI + AR(I,K)*ZI(K,J) + AI(I,K)*ZR(K,J))
110 CONTINUE
C .......... DOUBLE DIVISIONS AVOID POSSIBLE UNDERFLOW ..........
S = FLOP((S/H)/H)
SI = FLOP((SI/H)/H)
C
DO 120 K = 1, L
ZR(K,J) = FLOP(ZR(K,J) - S*AR(I,K) - SI*AI(I,K))
ZI(K,J) = FLOP(ZI(K,J) - SI*AR(I,K) + S*AI(I,K))
120 CONTINUE
C
130 CONTINUE
C
140 CONTINUE
C
200 RETURN
END
SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR,JOB)
C
INTEGER I,J,K,L,M,N,II,NM,MML,IERR
DOUBLE PRECISION D(N),E(N),Z(NM,N)
DOUBLE PRECISION B,C,F,G,P,R,S
DOUBLE PRECISION FLOP
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2,
C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON,
C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
C
C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD.
C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS
C FULL MATRIX TO TRIDIAGONAL FORM.
C
C ON INPUT.
C
C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C DIMENSION STATEMENT.
C
C N IS THE ORDER OF THE MATRIX.
C
C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C
C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY.
C
C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS
C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
C THE IDENTITY MATRIX.
C
C ON OUTPUT.
C
C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
C UNORDERED FOR INDICES 1,2,...,IERR-1.
C
C E HAS BEEN DESTROYED.
C
C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE,
C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
C EIGENVALUES.
C
C IERR IS SET TO
C ZERO FOR NORMAL RETURN,
C J IF THE J-TH EIGENVALUE HAS NOT BEEN
C DETERMINED AFTER 30 ITERATIONS.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
C
C*****
C MODIFIED BY C. MOLER TO ELIMINATE MACHEP 11/22/78
C MODIFIED TO ADD JOB PARAMETER 08/27/79
C*****
IERR = 0
IF (N .EQ. 1) GO TO 1001
C
DO 100 I = 2, N
100 E(I-1) = E(I)
C
E(N) = 0.0D0
C
DO 240 L = 1, N
J = 0
C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........
105 DO 110 M = L, N
IF (M .EQ. N) GO TO 120
C*****
P = FLOP(DABS(D(M)) + DABS(D(M+1)))
S = FLOP(P + DABS(E(M)))
IF (P .EQ. S) GO TO 120
C*****
110 CONTINUE
C
120 P = D(L)
IF (M .EQ. L) GO TO 240
IF (J .EQ. 30) GO TO 1000
J = J + 1
C .......... FORM SHIFT ..........
G = FLOP((D(L+1) - P)/(2.0D0*E(L)))
R = FLOP(DSQRT(G*G+1.0D0))
G = FLOP(D(M) - P + E(L)/(G + DSIGN(R,G)))
S = 1.0D0
C = 1.0D0
P = 0.0D0
MML = M - L
C .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
DO 200 II = 1, MML
I = M - II
F = FLOP(S*E(I))
B = FLOP(C*E(I))
IF (DABS(F) .LT. DABS(G)) GO TO 150
C = FLOP(G/F)
R = FLOP(DSQRT(C*C+1.0D0))
E(I+1) = FLOP(F*R)
S = FLOP(1.0D0/R)
C = FLOP(C*S)
GO TO 160
150 S = FLOP(F/G)
R = FLOP(DSQRT(S*S+1.0D0))
E(I+1) = FLOP(G*R)
C = FLOP(1.0D0/R)
S = FLOP(S*C)
160 G = FLOP(D(I+1) - P)
R = FLOP((D(I) - G)*S + 2.0D0*C*B)
P = FLOP(S*R)
D(I+1) = G + P
G = FLOP(C*R - B)
IF (JOB .EQ. 0) GO TO 185
C .......... FORM VECTOR ..........
DO 180 K = 1, N
F = Z(K,I+1)
Z(K,I+1) = FLOP(S*Z(K,I) + C*F)
Z(K,I) = FLOP(C*Z(K,I) - S*F)
180 CONTINUE
185 CONTINUE
C
200 CONTINUE
C
D(L) = FLOP(D(L) - P)
E(L) = G
E(M) = 0.0D0
GO TO 105
240 CONTINUE
C .......... ORDER EIGENVALUES AND EIGENVECTORS ..........
DO 300 II = 2, N
I = II - 1
K = I
P = D(I)
C
DO 260 J = II, N
IF (D(J) .GE. P) GO TO 260
K = J
P = D(J)
260 CONTINUE
C
IF (K .EQ. I) GO TO 300
D(K) = D(I)
D(I) = P
C
IF (JOB .EQ. 0) GO TO 285
DO 280 J = 1, N
P = Z(J,I)
Z(J,I) = Z(J,K)
Z(J,K) = P
280 CONTINUE
285 CONTINUE
C
300 CONTINUE
C
GO TO 1001
C .......... SET ERROR -- NO CONVERGENCE TO AN
C EIGENVALUE AFTER 30 ITERATIONS ..........
1000 IERR = L
1001 RETURN
END
SUBROUTINE CORTH(NM,N,LOW,IGH,AR,AI,ORTR,ORTI)
C
INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW
DOUBLE PRECISION AR(NM,N),AI(NM,N),ORTR(IGH),ORTI(IGH)
DOUBLE PRECISION F,G,H,FI,FR,SCALE
DOUBLE PRECISION FLOP,PYTHAG
C
C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
C THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968)
C BY MARTIN AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
C
C GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE
C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
C UNITARY SIMILARITY TRANSFORMATIONS.
C
C ON INPUT.
C
C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C DIMENSION STATEMENT.
C
C N IS THE ORDER OF THE MATRIX.
C
C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED,
C SET LOW=1, IGH=N.
C
C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX INPUT MATRIX.
C
C ON OUTPUT.
C
C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE HESSENBERG MATRIX. INFORMATION
C ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION
C IS STORED IN THE REMAINING TRIANGLES UNDER THE
C HESSENBERG MATRIX.
C
C ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE
C TRANSFORMATIONS. ONLY ELEMENTS LOW THROUGH IGH ARE USED.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
LA = IGH - 1
KP1 = LOW + 1
IF (LA .LT. KP1) GO TO 200
C
DO 180 M = KP1, LA
H = 0.0D0
ORTR(M) = 0.0D0
ORTI(M) = 0.0D0
SCALE = 0.0D0
C .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
DO 90 I = M, IGH
90 SCALE = FLOP(SCALE + DABS(AR(I,M-1)) + DABS(AI(I,M-1)))
C
IF (SCALE .EQ. 0.0D0) GO TO 180
MP = M + IGH
C .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
DO 100 II = M, IGH
I = MP - II
ORTR(I) = FLOP(AR(I,M-1)/SCALE)
ORTI(I) = FLOP(AI(I,M-1)/SCALE)
H = FLOP(H + ORTR(I)*ORTR(I) + ORTI(I)*ORTI(I))
100 CONTINUE
C
G = FLOP(DSQRT(H))
F = PYTHAG(ORTR(M),ORTI(M))
IF (F .EQ. 0.0D0) GO TO 103
H = FLOP(H + F*G)
G = FLOP(G/F)
ORTR(M) = FLOP((1.0D0 + G)*ORTR(M))
ORTI(M) = FLOP((1.0D0 + G)*ORTI(M))
GO TO 105
C
103 ORTR(M) = G
AR(M,M-1) = SCALE
C .......... FORM (I-(U*UT)/H)*A ..........
105 DO 130 J = M, N
FR = 0.0D0
FI = 0.0D0
C .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
DO 110 II = M, IGH
I = MP - II
FR = FLOP(FR + ORTR(I)*AR(I,J) + ORTI(I)*AI(I,J))
FI = FLOP(FI + ORTR(I)*AI(I,J) - ORTI(I)*AR(I,J))
110 CONTINUE
C
FR = FLOP(FR/H)
FI = FLOP(FI/H)
C
DO 120 I = M, IGH
AR(I,J) = FLOP(AR(I,J) - FR*ORTR(I) + FI*ORTI(I))
AI(I,J) = FLOP(AI(I,J) - FR*ORTI(I) - FI*ORTR(I))
120 CONTINUE
C
130 CONTINUE
C .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) ..........
DO 160 I = 1, IGH
FR = 0.0D0
FI = 0.0D0
C .......... FOR J=IGH STEP -1 UNTIL M DO -- ..........
DO 140 JJ = M, IGH
J = MP - JJ
FR = FLOP(FR + ORTR(J)*AR(I,J) - ORTI(J)*AI(I,J))
FI = FLOP(FI + ORTR(J)*AI(I,J) + ORTI(J)*AR(I,J))
140 CONTINUE
C
FR = FLOP(FR/H)
FI = FLOP(FI/H)
C
DO 150 J = M, IGH
AR(I,J) = FLOP(AR(I,J) - FR*ORTR(J) - FI*ORTI(J))
AI(I,J) = FLOP(AI(I,J) + FR*ORTI(J) - FI*ORTR(J))
150 CONTINUE
C
160 CONTINUE
C
ORTR(M) = FLOP(SCALE*ORTR(M))
ORTI(M) = FLOP(SCALE*ORTI(M))
AR(M,M-1) = FLOP(-G*AR(M,M-1))
AI(M,M-1) = FLOP(-G*AI(M,M-1))
180 CONTINUE
C
200 RETURN
END
SUBROUTINE COMQR3(NM,N,LOW,IGH,ORTR,ORTI,HR,HI,WR,WI,ZR,ZI,IERR
* ,JOB)
C*****
C MODIFICATION OF EISPACK COMQR2 TO ADD JOB PARAMETER
C JOB = 0 OUTPUT H = SCHUR TRIANGULAR FORM, Z NOT USED
C = 1 OUTPUT H = SCHUR FORM, Z = UNITARY SIMILARITY
C = 2 SAME AS COMQR2
C = 3 OUTPUT H = HESSENBERG FORM, Z = UNITARY SIMILARITY
C ALSO ELIMINATE MACHEP
C C. MOLER, 11/22/78 AND 09/14/80
C OVERFLOW CONTROL IN EIGENVECTOR BACKSUBSTITUTION, 3/16/82
C*****
C
INTEGER I,J,K,L,M,N,EN,II,JJ,LL,NM,NN,IGH,IP1,
X ITN,ITS,LOW,LP1,ENM1,IEND,IERR
DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N),
X ORTR(IGH),ORTI(IGH)
DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM
DOUBLE PRECISION FLOP,PYTHAG
C
C THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
C ALGOL PROCEDURE COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS
C AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
C THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
C (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.
C
C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR
C METHOD. THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX
C CAN ALSO BE FOUND IF CORTH HAS BEEN USED TO REDUCE
C THIS GENERAL MATRIX TO HESSENBERG FORM.
C
C ON INPUT.
C
C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C DIMENSION STATEMENT.
C
C N IS THE ORDER OF THE MATRIX.
C
C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED,
C SET LOW=1, IGH=N.
C
C ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
C FORMATIONS USED IN THE REDUCTION BY CORTH, IF PERFORMED.
C ONLY ELEMENTS LOW THROUGH IGH ARE USED. IF THE EIGENVECTORS
C OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND
C ORTI(J) TO 0.0D0 FOR THESE ELEMENTS.
C
C HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
C THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER
C INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE
C REDUCTION BY CORTH, IF PERFORMED. IF THE EIGENVECTORS OF
C THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE
C ARBITRARY.
C
C ON OUTPUT.
C
C ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI
C HAVE BEEN DESTROYED.
C
C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR
C EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
C FOR INDICES IERR+1,...,N.
C
C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVECTORS. THE EIGENVECTORS
C ARE UNNORMALIZED. IF AN ERROR EXIT IS MADE, NONE OF
C THE EIGENVECTORS HAS BEEN FOUND.
C
C IERR IS SET TO
C ZERO FOR NORMAL RETURN,
C J IF THE J-TH EIGENVALUE HAS NOT BEEN
C DETERMINED AFTER A TOTAL OF 30*N ITERATIONS.
C
C MODIFIED TO GET RID OF ALL COMPLEX ARITHMETIC, C. MOLER, 6/27/79.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
IERR = 0
C*****
IF (JOB .EQ. 0) GO TO 150
C*****
C .......... INITIALIZE EIGENVECTOR MATRIX ..........
DO 100 I = 1, N
C
DO 100 J = 1, N
ZR(I,J) = 0.0D0
ZI(I,J) = 0.0D0
IF (I .EQ. J) ZR(I,J) = 1.0D0
100 CONTINUE
C .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS
C FROM THE INFORMATION LEFT BY CORTH ..........
IEND = IGH - LOW - 1
IF (IEND) 180, 150, 105
C .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
105 DO 140 II = 1, IEND
I = IGH - II
IF (ORTR(I) .EQ. 0.0D0 .AND. ORTI(I) .EQ. 0.0D0) GO TO 140
IF (HR(I,I-1) .EQ. 0.0D0 .AND. HI(I,I-1) .EQ. 0.0D0) GO TO 140
C .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
NORM = FLOP(HR(I,I-1)*ORTR(I) + HI(I,I-1)*ORTI(I))
IP1 = I + 1
C
DO 110 K = IP1, IGH
ORTR(K) = HR(K,I-1)
ORTI(K) = HI(K,I-1)
110 CONTINUE
C
DO 130 J = I, IGH
SR = 0.0D0
SI = 0.0D0
C
DO 115 K = I, IGH
SR = FLOP(SR + ORTR(K)*ZR(K,J) + ORTI(K)*ZI(K,J))
SI = FLOP(SI + ORTR(K)*ZI(K,J) - ORTI(K)*ZR(K,J))
115 CONTINUE
C
SR = FLOP(SR/NORM)
SI = FLOP(SI/NORM)
C
DO 120 K = I, IGH
ZR(K,J) = FLOP(ZR(K,J) + SR*ORTR(K) - SI*ORTI(K))
ZI(K,J) = FLOP(ZI(K,J) + SR*ORTI(K) + SI*ORTR(K))
120 CONTINUE
C
130 CONTINUE
C
140 CONTINUE
C*****
IF (JOB .EQ. 3) GO TO 1001
C*****
C .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
150 L = LOW + 1
C
DO 170 I = L, IGH
LL = MIN0(I+1,IGH)
IF (HI(I,I-1) .EQ. 0.0D0) GO TO 170
NORM = PYTHAG(HR(I,I-1),HI(I,I-1))
YR = FLOP(HR(I,I-1)/NORM)
YI = FLOP(HI(I,I-1)/NORM)
HR(I,I-1) = NORM
HI(I,I-1) = 0.0D0
C
DO 155 J = I, N
SI = FLOP(YR*HI(I,J) - YI*HR(I,J))
HR(I,J) = FLOP(YR*HR(I,J) + YI*HI(I,J))
HI(I,J) = SI
155 CONTINUE
C
DO 160 J = 1, LL
SI = FLOP(YR*HI(J,I) + YI*HR(J,I))
HR(J,I) = FLOP(YR*HR(J,I) - YI*HI(J,I))
HI(J,I) = SI
160 CONTINUE
C*****
IF (JOB .EQ. 0) GO TO 170
C*****
DO 165 J = LOW, IGH
SI = FLOP(YR*ZI(J,I) + YI*ZR(J,I))
ZR(J,I) = FLOP(YR*ZR(J,I) - YI*ZI(J,I))
ZI(J,I) = SI
165 CONTINUE
C
170 CONTINUE
C .......... STORE ROOTS ISOLATED BY CBAL ..........
180 DO 200 I = 1, N
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
WR(I) = HR(I,I)
WI(I) = HI(I,I)
200 CONTINUE
C
EN = IGH
TR = 0.0D0
TI = 0.0D0
ITN = 30*N
C .......... SEARCH FOR NEXT EIGENVALUE ..........
220 IF (EN .LT. LOW) GO TO 680
ITS = 0
ENM1 = EN - 1
C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
C FOR L=EN STEP -1 UNTIL LOW DO -- ..........
240 DO 260 LL = LOW, EN
L = EN + LOW - LL
IF (L .EQ. LOW) GO TO 300
C*****
XR = FLOP(DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))
X + DABS(HR(L,L)) +DABS(HI(L,L)))
YR = FLOP(XR + DABS(HR(L,L-1)))
IF (XR .EQ. YR) GO TO 300
C*****
260 CONTINUE
C .......... FORM SHIFT ..........
300 IF (L .EQ. EN) GO TO 660
IF (ITN .EQ. 0) GO TO 1000
IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
SR = HR(EN,EN)
SI = HI(EN,EN)
XR = FLOP(HR(ENM1,EN)*HR(EN,ENM1))
XI = FLOP(HI(ENM1,EN)*HR(EN,ENM1))
IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 340
YR = FLOP((HR(ENM1,ENM1) - SR)/2.0D0)
YI = FLOP((HI(ENM1,ENM1) - SI)/2.0D0)
CALL WSQRT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)
IF (YR*ZZR + YI*ZZI .GE. 0.0D0) GO TO 310
ZZR = -ZZR
ZZI = -ZZI
310 CALL WDIV(XR,XI,YR+ZZR,YI+ZZI,ZZR,ZZI)
SR = FLOP(SR - ZZR)
SI = FLOP(SI - ZZI)
GO TO 340
C .......... FORM EXCEPTIONAL SHIFT ..........
320 SR = FLOP(DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2)))
SI = 0.0D0
C
340 DO 360 I = LOW, EN
HR(I,I) = FLOP(HR(I,I) - SR)
HI(I,I) = FLOP(HI(I,I) - SI)
360 CONTINUE
C
TR = FLOP(TR + SR)
TI = FLOP(TI + SI)
ITS = ITS + 1
ITN = ITN - 1
C .......... REDUCE TO TRIANGLE (ROWS) ..........
LP1 = L + 1
C
DO 500 I = LP1, EN
SR = HR(I,I-1)
HR(I,I-1) = 0.0D0
NORM = FLOP(DABS(HR(I-1,I-1)) + DABS(HI(I-1,I-1)) + DABS(SR))
NORM = FLOP(NORM*DSQRT((HR(I-1,I-1)/NORM)**2 +
X (HI(I-1,I-1)/NORM)**2 + (SR/NORM)**2))
XR = FLOP(HR(I-1,I-1)/NORM)
WR(I-1) = XR
XI = FLOP(HI(I-1,I-1)/NORM)
WI(I-1) = XI
HR(I-1,I-1) = NORM
HI(I-1,I-1) = 0.0D0
HI(I,I-1) = FLOP(SR/NORM)
C
DO 490 J = I, N
YR = HR(I-1,J)
YI = HI(I-1,J)
ZZR = HR(I,J)
ZZI = HI(I,J)
HR(I-1,J) = FLOP(XR*YR + XI*YI + HI(I,I-1)*ZZR)
HI(I-1,J) = FLOP(XR*YI - XI*YR + HI(I,I-1)*ZZI)
HR(I,J) = FLOP(XR*ZZR - XI*ZZI - HI(I,I-1)*YR)
HI(I,J) = FLOP(XR*ZZI + XI*ZZR - HI(I,I-1)*YI)
490 CONTINUE
C
500 CONTINUE
C
SI = HI(EN,EN)
IF (SI .EQ. 0.0D0) GO TO 540
NORM = PYTHAG(HR(EN,EN),SI)
SR = FLOP(HR(EN,EN)/NORM)
SI = FLOP(SI/NORM)
HR(EN,EN) = NORM
HI(EN,EN) = 0.0D0
IF (EN .EQ. N) GO TO 540
IP1 = EN + 1
C
DO 520 J = IP1, N
YR = HR(EN,J)
YI = HI(EN,J)
HR(EN,J) = FLOP(SR*YR + SI*YI)
HI(EN,J) = FLOP(SR*YI - SI*YR)
520 CONTINUE
C .......... INVERSE OPERATION (COLUMNS) ..........
540 DO 600 J = LP1, EN
XR = WR(J-1)
XI = WI(J-1)
C
DO 580 I = 1, J
YR = HR(I,J-1)
YI = 0.0D0
ZZR = HR(I,J)
ZZI = HI(I,J)
IF (I .EQ. J) GO TO 560
YI = HI(I,J-1)
HI(I,J-1) = FLOP(XR*YI + XI*YR + HI(J,J-1)*ZZI)
560 HR(I,J-1) = FLOP(XR*YR - XI*YI + HI(J,J-1)*ZZR)
HR(I,J) = FLOP(XR*ZZR + XI*ZZI - HI(J,J-1)*YR)
HI(I,J) = FLOP(XR*ZZI - XI*ZZR - HI(J,J-1)*YI)
580 CONTINUE
C*****
IF (JOB .EQ. 0) GO TO 600
C*****
DO 590 I = LOW, IGH
YR = ZR(I,J-1)
YI = ZI(I,J-1)
ZZR = ZR(I,J)
ZZI = ZI(I,J)
ZR(I,J-1) = FLOP(XR*YR - XI*YI + HI(J,J-1)*ZZR)
ZI(I,J-1) = FLOP(XR*YI + XI*YR + HI(J,J-1)*ZZI)
ZR(I,J) = FLOP(XR*ZZR + XI*ZZI - HI(J,J-1)*YR)
ZI(I,J) = FLOP(XR*ZZI - XI*ZZR - HI(J,J-1)*YI)
590 CONTINUE
C
600 CONTINUE
C
IF (SI .EQ. 0.0D0) GO TO 240
C
DO 630 I = 1, EN
YR = HR(I,EN)
YI = HI(I,EN)
HR(I,EN) = FLOP(SR*YR - SI*YI)
HI(I,EN) = FLOP(SR*YI + SI*YR)
630 CONTINUE
C*****
IF (JOB .EQ. 0) GO TO 240
C*****
DO 640 I = LOW, IGH
YR = ZR(I,EN)
YI = ZI(I,EN)
ZR(I,EN) = FLOP(SR*YR - SI*YI)
ZI(I,EN) = FLOP(SR*YI + SI*YR)
640 CONTINUE
C
GO TO 240
C .......... A ROOT FOUND ..........
660 HR(EN,EN) = FLOP(HR(EN,EN) + TR)
WR(EN) = HR(EN,EN)
HI(EN,EN) = FLOP(HI(EN,EN) + TI)
WI(EN) = HI(EN,EN)
EN = ENM1
GO TO 220
C .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND
C VECTORS OF UPPER TRIANGULAR FORM ..........
C
C***** THE FOLLOWING SECTION CHANGED FOR OVERFLOW CONTROL
C C. MOLER, 3/16/82
C
680 IF (JOB .NE. 2) GO TO 1001
C
NORM = 0.0D0
DO 720 I = 1, N
DO 720 J = I, N
TR = FLOP(DABS(HR(I,J))) + FLOP(DABS(HI(I,J)))
IF (TR .GT. NORM) NORM = TR
720 CONTINUE
IF (N .EQ. 1 .OR. NORM .EQ. 0.0D0) GO TO 1001
C .......... FOR EN=N STEP -1 UNTIL 2 DO -- ..........
DO 800 NN = 2, N
EN = N + 2 - NN
XR = WR(EN)
XI = WI(EN)
HR(EN,EN) = 1.0D0
HI(EN,EN) = 0.0D0
ENM1 = EN - 1
C .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
DO 780 II = 1, ENM1
I = EN - II
ZZR = 0.0D0
ZZI = 0.0D0
IP1 = I + 1
DO 740 J = IP1, EN
ZZR = FLOP(ZZR + HR(I,J)*HR(J,EN) - HI(I,J)*HI(J,EN))
ZZI = FLOP(ZZI + HR(I,J)*HI(J,EN) + HI(I,J)*HR(J,EN))
740 CONTINUE
YR = FLOP(XR - WR(I))
YI = FLOP(XI - WI(I))
IF (YR .NE. 0.0D0 .OR. YI .NE. 0.0D0) GO TO 765
YR = NORM
760 YR = FLOP(YR/100.0D0)
YI = FLOP(NORM + YR)
IF (YI .NE. NORM) GO TO 760
YI = 0.0D0
765 CONTINUE
CALL WDIV(ZZR,ZZI,YR,YI,HR(I,EN),HI(I,EN))
TR = FLOP(DABS(HR(I,EN))) + FLOP(DABS(HI(I,EN)))
IF (TR .EQ. 0.0D0) GO TO 780
IF (TR + 1.0D0/TR .GT. TR) GO TO 780
DO 770 J = I, EN
HR(J,EN) = FLOP(HR(J,EN)/TR)
HI(J,EN) = FLOP(HI(J,EN)/TR)
770 CONTINUE
780 CONTINUE
C
800 CONTINUE
C*****
C .......... END BACKSUBSTITUTION ..........
ENM1 = N - 1
C .......... VECTORS OF ISOLATED ROOTS ..........
DO 840 I = 1, ENM1
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
IP1 = I + 1
C
DO 820 J = IP1, N
ZR(I,J) = HR(I,J)
ZI(I,J) = HI(I,J)
820 CONTINUE
C
840 CONTINUE
C .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
C VECTORS OF ORIGINAL FULL MATRIX.
C FOR J=N STEP -1 UNTIL LOW+1 DO -- ..........
DO 880 JJ = LOW, ENM1
J = N + LOW - JJ
M = MIN0(J,IGH)
C
DO 880 I = LOW, IGH
ZZR = 0.0D0
ZZI = 0.0D0
C
DO 860 K = LOW, M
ZZR = FLOP(ZZR + ZR(I,K)*HR(K,J) - ZI(I,K)*HI(K,J))
ZZI = FLOP(ZZI + ZR(I,K)*HI(K,J) + ZI(I,K)*HR(K,J))
860 CONTINUE
C
ZR(I,J) = ZZR
ZI(I,J) = ZZI
880 CONTINUE
C
GO TO 1001
C .......... SET ERROR -- NO CONVERGENCE TO AN
C EIGENVALUE AFTER 30 ITERATIONS ..........
1000 IERR = EN
1001 RETURN
END
SUBROUTINE WSVDC(XR,XI,LDX,N,P,SR,SI,ER,EI,UR,UI,LDU,VR,VI,LDV,
* WORKR,WORKI,JOB,INFO)
INTEGER LDX,N,P,LDU,LDV,JOB,INFO
DOUBLE PRECISION XR(LDX,1),XI(LDX,1),SR(1),SI(1),ER(1),EI(1),
* UR(LDU,1),UI(LDU,1),VR(LDV,1),VI(LDV,1),
* WORKR(1),WORKI(1)
C
C
C WSVDC IS A SUBROUTINE TO REDUCE A DOUBLE-COMPLEX NXP MATRIX X BY
C UNITARY TRANSFORMATIONS U AND V TO DIAGONAL FORM. THE
C DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X. THE
C COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
C AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
C
C ON ENTRY
C
C X DOUBLE-COMPLEX(LDX,P), WHERE LDX.GE.N.
C X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
C DECOMPOSITION IS TO BE COMPUTED. X IS
C DESTROYED BY WSVDC.
C
C LDX INTEGER.
C LDX IS THE LEADING DIMENSION OF THE ARRAY X.
C
C N INTEGER.
C N IS THE NUMBER OF COLUMNS OF THE MATRIX X.
C
C P INTEGER.
C P IS THE NUMBER OF ROWS OF THE MATRIX X.
C
C LDU INTEGER.
C LDU IS THE LEADING DIMENSION OF THE ARRAY U
C (SEE BELOW).
C
C LDV INTEGER.
C LDV IS THE LEADING DIMENSION OF THE ARRAY V
C (SEE BELOW).
C
C WORK DOUBLE-COMPLEX(N).
C WORK IS A SCRATCH ARRAY.
C
C JOB INTEGER.
C JOB CONTROLS THE COMPUTATION OF THE SINGULAR
C VECTORS. IT HAS THE DECIMAL EXPANSION AB
C WITH THE FOLLOWING MEANING
C
C A.EQ.0 DO NOT COMPUTE THE LEFT SINGULAR
C VECTORS.
C A.EQ.1 RETURN THE N LEFT SINGULAR VECTORS
C IN U.
C A.GE.2 RETURNS THE FIRST MIN(N,P)
C LEFT SINGULAR VECTORS IN U.
C B.EQ.0 DO NOT COMPUTE THE RIGHT SINGULAR
C VECTORS.
C B.EQ.1 RETURN THE RIGHT SINGULAR VECTORS
C IN V.
C
C ON RETURN
C
C S DOUBLE-COMPLEX(MM), WHERE MM=MIN(N+1,P).
C THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
C SINGULAR VALUES OF X ARRANGED IN DESCENDING
C ORDER OF MAGNITUDE.
C
C E DOUBLE-COMPLEX(P).
C E ORDINARILY CONTAINS ZEROS. HOWEVER SEE THE
C DISCUSSION OF INFO FOR EXCEPTIONS.
C
C U DOUBLE-COMPLEX(LDU,K), WHERE LDU.GE.N.
C IF JOBA.EQ.1 THEN K.EQ.N,
C IF JOBA.EQ.2 THEN K.EQ.MIN(N,P).
C U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
C U IS NOT REFERENCED IF JOBA.EQ.0. IF N.LE.P
C OR IF JOBA.GT.2, THEN U MAY BE IDENTIFIED WITH X
C IN THE SUBROUTINE CALL.
C
C V DOUBLE-COMPLEX(LDV,P), WHERE LDV.GE.P.
C V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
C V IS NOT REFERENCED IF JOBB.EQ.0. IF P.LE.N,
C THEN V MAY BE IDENTIFIED WHTH X IN THE
C SUBROUTINE CALL.
C
C INFO INTEGER.
C THE SINGULAR VALUES (AND THEIR CORRESPONDING
C SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
C ARE CORRECT (HERE M=MIN(N,P)). THUS IF
C INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
C VECTORS ARE CORRECT. IN ANY EVENT, THE MATRIX
C B = CTRANS(U)*X*V IS THE BIDIAGONAL MATRIX
C WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
C ELEMENTS OF E ON ITS SUPER-DIAGONAL (CTRANS(U)
C IS THE CONJUGATE-TRANSPOSE OF U). THUS THE
C SINGULAR VALUES OF X AND B ARE THE SAME.
C
C LINPACK. THIS VERSION DATED 07/03/79 .
C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C WSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
C
C BLAS WAXPY,PYTHAG,WDOTCR,WDOTCI,WSCAL,WSWAP,WNRM2,RROTG
C FORTRAN DABS,DIMAG,DMAX1
C FORTRAN MAX0,MIN0,MOD,DSQRT
C
C INTERNAL VARIABLES
C
INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
* MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
DOUBLE PRECISION PYTHAG,WDOTCR,WDOTCI,TR,TI,RR,RI
DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,WNRM2,SCALE,SHIFT,SL,SM,SN,
* SMM1,T1,TEST,ZTEST,SMALL,FLOP
LOGICAL WANTU,WANTV
C
DOUBLE PRECISION ZDUMR,ZDUMI
DOUBLE PRECISION CABS1
CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)
C
C SET THE MAXIMUM NUMBER OF ITERATIONS.
C
MAXIT = 75
C
C SMALL NUMBER, ROUGHLY MACHINE EPSILON, USED TO AVOID UNDERFLOW
C
SHAR_EOF
# End of shell archive
exit 0
--
Bob Page, U of Lowell CS Dept. page@swan.ulowell.edu ulowell!page
Have five nice days.