277 lines
8.3 KiB
FortranFixed
277 lines
8.3 KiB
FortranFixed
SUBROUTINE ORTHES(NM,N,LOW,IGH,A,ORT)
|
|
C
|
|
C
|
|
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTHES,
|
|
C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
|
|
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
|
|
C
|
|
C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE
|
|
C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
|
|
C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
|
|
C ORTHOGONAL 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 BALANC. IF BALANC HAS NOT BEEN USED,
|
|
C SET LOW=1, IGH=N.
|
|
C
|
|
C A CONTAINS THE INPUT MATRIX.
|
|
C
|
|
C ON OUTPUT
|
|
C
|
|
C A CONTAINS THE HESSENBERG MATRIX. INFORMATION ABOUT
|
|
C THE ORTHOGONAL TRANSFORMATIONS USED IN THE REDUCTION
|
|
C IS STORED IN THE REMAINING TRIANGLE UNDER THE
|
|
C HESSENBERG MATRIX.
|
|
C
|
|
C ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
|
|
C ONLY ELEMENTS LOW THROUGH IGH ARE USED.
|
|
C
|
|
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
|
|
C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
|
|
C
|
|
C THIS VERSION DATED AUGUST 1983.
|
|
C
|
|
C ------------------------------------------------------------------
|
|
C
|
|
C .. Scalar Arguments ..
|
|
INTEGER IGH,LOW,N,NM
|
|
C ..
|
|
C .. Array Arguments ..
|
|
DOUBLE PRECISION A(NM,N),ORT(IGH)
|
|
C ..
|
|
C .. Local Scalars ..
|
|
DOUBLE PRECISION F,G,H,SCALE
|
|
INTEGER I,II,J,JJ,KP1,LA,M,MP
|
|
C ..
|
|
C .. Intrinsic Functions ..
|
|
INTRINSIC DABS,DSIGN,DSQRT
|
|
C ..
|
|
LA = IGH - 1
|
|
KP1 = LOW + 1
|
|
IF (LA.GE.KP1) THEN
|
|
C
|
|
DO 90 M = KP1,LA
|
|
H = 0.0D0
|
|
ORT(M) = 0.0D0
|
|
SCALE = 0.0D0
|
|
C .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
|
|
DO 10 I = M,IGH
|
|
SCALE = SCALE + DABS(A(I,M-1))
|
|
10 CONTINUE
|
|
C
|
|
IF (SCALE.NE.0.0D0) THEN
|
|
MP = M + IGH
|
|
C .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
|
|
DO 20 II = M,IGH
|
|
I = MP - II
|
|
ORT(I) = A(I,M-1)/SCALE
|
|
H = H + ORT(I)*ORT(I)
|
|
20 CONTINUE
|
|
C
|
|
G = -DSIGN(DSQRT(H),ORT(M))
|
|
H = H - ORT(M)*G
|
|
ORT(M) = ORT(M) - G
|
|
C .......... FORM (I-(U*UT)/H) * A ..........
|
|
DO 50 J = M,N
|
|
F = 0.0D0
|
|
C .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
|
|
DO 30 II = M,IGH
|
|
I = MP - II
|
|
F = F + ORT(I)*A(I,J)
|
|
30 CONTINUE
|
|
C
|
|
F = F/H
|
|
C
|
|
DO 40 I = M,IGH
|
|
A(I,J) = A(I,J) - F*ORT(I)
|
|
40 CONTINUE
|
|
50 CONTINUE
|
|
C
|
|
C .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) ..........
|
|
DO 80 I = 1,IGH
|
|
F = 0.0D0
|
|
C .......... FOR J=IGH STEP -1 UNTIL M DO -- ..........
|
|
DO 60 JJ = M,IGH
|
|
J = MP - JJ
|
|
F = F + ORT(J)*A(I,J)
|
|
60 CONTINUE
|
|
C
|
|
F = F/H
|
|
C
|
|
DO 70 J = M,IGH
|
|
A(I,J) = A(I,J) - F*ORT(J)
|
|
70 CONTINUE
|
|
80 CONTINUE
|
|
C
|
|
C
|
|
ORT(M) = SCALE*ORT(M)
|
|
A(M,M-1) = SCALE*G
|
|
END IF
|
|
90 CONTINUE
|
|
END IF
|
|
RETURN
|
|
*$st$ Unreachable comments ...
|
|
C
|
|
END
|
|
SUBROUTINE TRED1(NM,N,A,D,E,E2)
|
|
C
|
|
C
|
|
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1,
|
|
C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
|
|
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
|
|
C
|
|
C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX
|
|
C TO A SYMMETRIC TRIDIAGONAL MATRIX USING
|
|
C ORTHOGONAL 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 A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE
|
|
C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
|
|
C
|
|
C ON OUTPUT
|
|
C
|
|
C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
|
|
C FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER
|
|
C TRIANGLE. THE FULL UPPER TRIANGLE OF A IS UNALTERED.
|
|
C
|
|
C D CONTAINS THE DIAGONAL ELEMENTS OF 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 QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
|
|
C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
|
|
C
|
|
C THIS VERSION DATED AUGUST 1983.
|
|
C
|
|
C ------------------------------------------------------------------
|
|
C
|
|
C .. Scalar Arguments ..
|
|
INTEGER N,NM
|
|
C ..
|
|
C .. Array Arguments ..
|
|
DOUBLE PRECISION A(NM,N),D(N),E(N),E2(N)
|
|
C ..
|
|
C .. Local Scalars ..
|
|
DOUBLE PRECISION F,G,H,SCALE
|
|
INTEGER I,II,J,JP1,K,L
|
|
C ..
|
|
C .. Intrinsic Functions ..
|
|
INTRINSIC DABS,DSIGN,DSQRT
|
|
C ..
|
|
DO 10 I = 1,N
|
|
D(I) = A(N,I)
|
|
A(N,I) = A(I,I)
|
|
10 CONTINUE
|
|
C .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
|
|
DO 130 II = 1,N
|
|
I = N + 1 - II
|
|
L = I - 1
|
|
H = 0.0D0
|
|
SCALE = 0.0D0
|
|
IF (L.GE.1) THEN
|
|
C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
|
|
DO 20 K = 1,L
|
|
SCALE = SCALE + DABS(D(K))
|
|
20 CONTINUE
|
|
C
|
|
IF (SCALE.NE.0.0D0) THEN
|
|
C
|
|
DO 30 K = 1,L
|
|
D(K) = D(K)/SCALE
|
|
H = H + D(K)*D(K)
|
|
30 CONTINUE
|
|
C
|
|
E2(I) = SCALE*SCALE*H
|
|
F = D(L)
|
|
G = -DSIGN(DSQRT(H),F)
|
|
E(I) = SCALE*G
|
|
H = H - F*G
|
|
D(L) = F - G
|
|
IF (L.NE.1) THEN
|
|
C .......... FORM A*U ..........
|
|
DO 40 J = 1,L
|
|
E(J) = 0.0D0
|
|
40 CONTINUE
|
|
C
|
|
DO 60 J = 1,L
|
|
F = D(J)
|
|
G = E(J) + A(J,J)*F
|
|
JP1 = J + 1
|
|
IF (L.GE.JP1) THEN
|
|
C
|
|
DO 50 K = JP1,L
|
|
G = G + A(K,J)*D(K)
|
|
E(K) = E(K) + A(K,J)*F
|
|
50 CONTINUE
|
|
END IF
|
|
C
|
|
E(J) = G
|
|
60 CONTINUE
|
|
C .......... FORM P ..........
|
|
F = 0.0D0
|
|
C
|
|
DO 70 J = 1,L
|
|
E(J) = E(J)/H
|
|
F = F + E(J)*D(J)
|
|
70 CONTINUE
|
|
C
|
|
H = F/ (H+H)
|
|
C .......... FORM Q ..........
|
|
DO 80 J = 1,L
|
|
E(J) = E(J) - H*D(J)
|
|
80 CONTINUE
|
|
C .......... FORM REDUCED A ..........
|
|
DO 100 J = 1,L
|
|
F = D(J)
|
|
G = E(J)
|
|
C
|
|
DO 90 K = J,L
|
|
A(K,J) = A(K,J) - F*E(K) - G*D(K)
|
|
90 CONTINUE
|
|
100 CONTINUE
|
|
C
|
|
END IF
|
|
C
|
|
DO 110 J = 1,L
|
|
F = D(J)
|
|
D(J) = A(L,J)
|
|
A(L,J) = A(I,J)
|
|
A(I,J) = F*SCALE
|
|
110 CONTINUE
|
|
GO TO 130
|
|
ELSE
|
|
C
|
|
DO 120 J = 1,L
|
|
D(J) = A(L,J)
|
|
A(L,J) = A(I,J)
|
|
A(I,J) = 0.0D0
|
|
120 CONTINUE
|
|
END IF
|
|
END IF
|
|
C
|
|
E(I) = 0.0D0
|
|
E2(I) = 0.0D0
|
|
130 CONTINUE
|
|
C
|
|
C
|
|
END
|