real a(500,500),aa(500,500),d(500),e(500),e2(500)
real z(500,500),zz(500,500)
real dd(500),ee(500)
real yy(500),y2(500),d2(500)
nm=500
do 100 n=50,500,50
do 20 i=1,n
do 20 j=1,i
a(i,j)=1./float(i+j-1)
aa(j,i)=a(i,j)
a(j,i)=a(i,j)
aa(i,j)=a(i,j)
20 continue
t1=second(0)
call tred2(nm,n,a,d,e,z)
t2=second(0)
call tred2v(nm,n,aa,dd,ee,zz,yy,d2,y2)
t3=second(0)
t3=t3-t2
t2=t2-t1
write(6,45)n,t2,t3
45 format(' n, old, new times for tred',i5,2e15.5)
c do 55 i=1,n
c write(6,56)d(i),dd(i),e(i),ee(i)
c56 format(' d,dd,ee,ee',4e15.5)
c55 continue
c do 70 i=1,n
c write(6,66)(z(i,j),j=1,n)
c66 format(' old z',5e15.5)
c write(6,67)(zz(i,j),j=1,n)
c67 format(' new z', 5e15.5)
c70 continue
100 continue
stop
end
subroutine symtr6(a,lda,n,x,y)
real a(lda,n),x(n),y(n)
do 30 j=1,n,2
xj=x(j)
if (j.eq.n) go to 27
xjp1=x(j+1)
do 20 i=j+1,n
y(i)=(y(i)+a(i,j)*xj)+a(i,j+1)*xjp1
20 continue
y(j)=y(j)+a(j+1,j)*xjp1+a(j,j)*xj
if (j.eq.1) go to 30
do 25 i=1,j-1
25 y(i)=(y(i)+a(j,i)*xj)+a(j+1,i)*xjp1
go to 30
27 do 28 i=1,n
28 y(i)=y(i)+a(n,i)*xj
30 continue
return
end
SUBROUTINE SXMPY(N1,LDY,Y,N2,LDX,X,LDM,M)
C ** SXMPY16.F -- SXMPY UNROLLED TO A DEPTH OF 16
C
INTEGER N1,LDY,N2,LDX,LDM
REAL Y(LDY,N1),X(LDX,N2),M(LDM,N1)
C
C Cleanup odd vector
C
J = MOD(N2,2)
IF (J .GE. 1) THEN
DO 10 I = 1, N1
Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I)
10 CONTINUE
ENDIF
C
C Cleanup odd group of two vectors
C
J = MOD(N2,4)
IF (J .GE. 2) THEN
DO 20 I = 1, N1
Y(1,I) = ( (Y(1,I))
$ + X(1,J-1)*M(J-1,I)) + X(1,J)*M(J,I)
20 CONTINUE
ENDIF
C
C Cleanup odd group of four vectors
C
J = MOD(N2,8)
IF (J .GE. 4) THEN
DO 30 I = 1, N1
Y(1,I) = ((( (Y(1,I))
$ + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I))
$ + X(1,J-1)*M(J-1,I)) + X(1,J) *M(J,I)
30 CONTINUE
ENDIF
C
C Cleanup odd group of eight vectors
C
J = MOD(N2,16)
IF (J .GE. 8) THEN
DO 40 I = 1, N1
Y(1,I) = ((((((( (Y(1,I))
$ + X(1,J-7)*M(J-7,I)) + X(1,J-6)*M(J-6,I))
$ + X(1,J-5)*M(J-5,I)) + X(1,J-4)*M(J-4,I))
$ + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I))
$ + X(1,J-1)*M(J-1,I)) + X(1,J) *M(J,I)
40 CONTINUE
ENDIF
C
C Main loop - groups of sixteen vectors
C
JMIN = J+16
DO 60 J = JMIN, N2, 16
DO 50 I = 1, N1
Y(1,I) = ((((((((((((((( (Y(1,I))
$ + X(1,J-15)*M(J-15,I)) + X(1,J-14)*M(J-14,I))
$ + X(1,J-13)*M(J-13,I)) + X(1,J-12)*M(J-12,I))
$ + X(1,J-11)*M(J-11,I)) + X(1,J-10)*M(J-10,I))
$ + X(1,J- 9)*M(J- 9,I)) + X(1,J- 8)*M(J- 8,I))
$ + X(1,J- 7)*M(J- 7,I)) + X(1,J- 6)*M(J- 6,I))
$ + X(1,J- 5)*M(J- 5,I)) + X(1,J- 4)*M(J- 4,I))
$ + X(1,J- 3)*M(J- 3,I)) + X(1,J- 2)*M(J- 2,I))
$ + X(1,J- 1)*M(J- 1,I)) + X(1,J) *M(J,I)
50 CONTINUE
60 CONTINUE
RETURN
END
SUBROUTINE TRED2v(NM,N,A,D,E,Z,yy,d2,y2)
C
real d2(n),y2(n)
INTEGER I,J,K,L,N,II,NM,JP1
REAL A(NM,N),D(N),E(N),Z(NM,N) ,yy(n)
REAL F,G,H,HH,SCALE
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
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 TO A
C SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
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 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 Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
C PRODUCED IN THE REDUCTION.
C
C A AND Z MAY COINCIDE. IF DISTINCT, A IS UNALTERED.
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 APRIL 1983.
c this is a modification of the eispack routine
c this routine should perform better than tred2 from
c eispack. To make this routine run even faster implement
c symtr6 and sxmpy in assembler.
c
c Jack Dongarra and Linda Kaufman
C
C ------------------------------------------------------------------
C
DO 100 I = 1, N
C
DO 80 J = I, N
80 Z(J,I) = A(J,I)
C
D(I) = A(N,I)
100 CONTINUE
C
IF (N .EQ. 1) GO TO 510
C .......... FOR I=N STEP -1 UNTIL 2 DO -- ..........
DO 300 II = 2, N
I = N + 2 - II
L = I - 1
H = 0.0E0
SCALE = 0.0E0
IF (L .LT. 2) GO TO 130
C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
DO 120 K = 1, L
120 SCALE = SCALE + ABS(D(K))
C
IF (SCALE .NE. 0.0E0) GO TO 140
130 E(I) = D(L)
C
DO 135 J = 1, L
D(J) = Z(L,J)
Z(I,J) = 0.0E0
Z(J,I) = 0.0E0
135 CONTINUE
C
GO TO 290
C
140 DO 150 K = 1, L
D(K) = D(K) / SCALE
H = H + D(K) * D(K)
150 CONTINUE
C
F = D(L)
G = -SIGN(SQRT(H),F)
E(I) = SCALE * G
H = H - F * G
D(L) = F - G
C .......... FORM A*U ..........
DO 170 J = 1, L
z(j,i)=d(j)
170 E(J) = 0.0E0
C
call symtr6(z,nm,l,d,e)
C .......... FORM P ..........
F = 0.0E0
C
DO 245 J = 1, L
E(J) = E(J) / H
F = F + E(J) * D(J)
245 CONTINUE
C
HH = F / (H + H)
C .......... FORM Q ..........
DO 250 J = 1, L
250 E(J) = E(J) - HH * D(J)
C .......... FORM REDUCED A ..........
DO 280 J = 1, L
F = D(J)
G = E(J)
C
DO 260 K = J, L
260 Z(K,J) = Z(K,J) - F * E(K) - G * D(K)
C
D(J) = Z(L,J)
Z(I,J) = 0.0E0
280 CONTINUE
C
290 D(I) = H
300 CONTINUE
C .......... ACCUMULATION OF TRANSFORMATION MATRICES ..........
DO 500 I = 2, N , 2
L = I - 1
Z(N,L) = Z(L,L)
Z(L,L) = 1.0E0
h=0.0
if (d(i).ne.0.0)h=1.0/d(i)
C
DO 330 K = 1, L
330 D(K) = Z(K,I) *H
C
do 335 j=1,l
335 yy(j)=0.0
call sxmpy(l,1,yy,l,1,z(1,i),nm,z)
if (i.lt.n)go to 410
do 355 j=1,l
g =yy(j)
do 350 k=1,l
350 z(k,j)=z(k,j)-g*d(k)
355 continue
c DO 360 J = 1, L
c G = 0.0E0
cC
c DO 340 K = 1, L
c 340 G = G + Z(K,I) * Z(K,J)
cC
c DO 360 K = 1, L
c Z(K,J) = Z(K,J) - G * D(K)
c 360 CONTINUE
C
do 400 k=1,l
400 z(k,i)=0.0e0
go to 500
410 z(n,i)=z(i,i)
ip1=i+1
hi=0.0
if (d(i+1).ne.0.0) hi=1.0/d(ip1)
do 430 k=1,i
430 d2(k)=z(k,ip1)*hi
do 435 j=1,i
z(j,i)=0.0
435 y2(j)=0.0
z(i,i)=1.0e0
call sxmpy(i,1,y2,i,1,z(1,ip1),nm,z)
vtw=0.0
do 445 j=1,l
445 vtw=vtw+d(j)*d2(j)
vtw=vtw*hi
if (hi.ne.0.0) vtw=vtw/(hi*hi)
yy(i)=0.0
d(i)=0.0
do 455 j=1,i
g=-yy(j)
g2=-y2(j)+vtw*yy(j)
do 450 k=1,i
450 z(k,j)=z(k,j)+g*d(k)+g2*d2(k)
455 continue
do 460 k=1,i
460 z(k,ip1)=0.0
C
500 CONTINUE
C
510 DO 520 I = 1, N
D(I) = Z(N,I)
Z(N,I) = 0.0E0
520 CONTINUE
C
Z(N,N) = 1.0E0
E(1) = 0.0E0
RETURN
END
SUBROUTINE TRED2(NM,N,A,D,E,Z)
C
INTEGER I,J,K,L,N,II,NM,JP1
REAL A(NM,N),D(N),E(N),Z(NM,N)
REAL F,G,H,HH,SCALE
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
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 TO A
C SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
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 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 Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
C PRODUCED IN THE REDUCTION.
C
C A AND Z MAY COINCIDE. IF DISTINCT, A IS UNALTERED.
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 APRIL 1983.
C
C ------------------------------------------------------------------
C
DO 100 I = 1, N
C
DO 80 J = I, N
80 Z(J,I) = A(J,I)
C
D(I) = A(N,I)
100 CONTINUE
C
IF (N .EQ. 1) GO TO 510
C .......... FOR I=N STEP -1 UNTIL 2 DO -- ..........
DO 300 II = 2, N
I = N + 2 - II
L = I - 1
H = 0.0E0
SCALE = 0.0E0
IF (L .LT. 2) GO TO 130
C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
DO 120 K = 1, L
120 SCALE = SCALE + ABS(D(K))
C
IF (SCALE .NE. 0.0E0) GO TO 140
130 E(I) = D(L)
C
DO 135 J = 1, L
D(J) = Z(L,J)
Z(I,J) = 0.0E0
Z(J,I) = 0.0E0
135 CONTINUE
C
GO TO 290
C
140 DO 150 K = 1, L
D(K) = D(K) / SCALE
H = H + D(K) * D(K)
150 CONTINUE
C
F = D(L)
G = -SIGN(SQRT(H),F)
E(I) = SCALE * G
H = H - F * G
D(L) = F - G
C .......... FORM A*U ..........
DO 170 J = 1, L
170 E(J) = 0.0E0
C
DO 240 J = 1, L
F = D(J)
Z(J,I) = F
G = E(J) + Z(J,J) * F
JP1 = J + 1
IF (L .LT. JP1) GO TO 220
C
DO 200 K = JP1, L
G = G + Z(K,J) * D(K)
E(K) = E(K) + Z(K,J) * F
200 CONTINUE
C
220 E(J) = G
240 CONTINUE
C .......... FORM P ..........
F = 0.0E0
C
DO 245 J = 1, L
E(J) = E(J) / H
F = F + E(J) * D(J)
245 CONTINUE
C
HH = F / (H + H)
C .......... FORM Q ..........
DO 250 J = 1, L
250 E(J) = E(J) - HH * D(J)
C .......... FORM REDUCED A ..........
DO 280 J = 1, L
F = D(J)
G = E(J)
C
DO 260 K = J, L
260 Z(K,J) = Z(K,J) - F * E(K) - G * D(K)
C
D(J) = Z(L,J)
Z(I,J) = 0.0E0
280 CONTINUE
C
290 D(I) = H
300 CONTINUE
C .......... ACCUMULATION OF TRANSFORMATION MATRICES ..........
DO 500 I = 2, N
L = I - 1
Z(N,L) = Z(L,L)
Z(L,L) = 1.0E0
H = D(I)
IF (H .EQ. 0.0E0) GO TO 380
C
DO 330 K = 1, L
330 D(K) = Z(K,I) / H
C
DO 360 J = 1, L
G = 0.0E0
C
DO 340 K = 1, L
340 G = G + Z(K,I) * Z(K,J)
C
DO 360 K = 1, L
Z(K,J) = Z(K,J) - G * D(K)
360 CONTINUE
C
380 DO 400 K = 1, L
400 Z(K,I) = 0.0E0
C
500 CONTINUE
C
510 DO 520 I = 1, N
D(I) = Z(N,I)
Z(N,I) = 0.0E0
520 CONTINUE
C
Z(N,N) = 1.0E0
E(1) = 0.0E0
RETURN
END
SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR)
C
INTEGER I,J,K,L,M,N,II,NM,MML,IERR
REAL D(N),E(N),Z(NM,N)
REAL B,C,F,G,P,R,S,TST1,TST2,PYTHAG
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 CALLS PYTHAG FOR SQRT(A*A + B*B) .
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 APRIL 1983.
C
C ------------------------------------------------------------------
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.0E0
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
TST1 = ABS(D(M)) + ABS(D(M+1))
TST2 = TST1 + ABS(E(M))
IF (TST2 .EQ. TST1) GO TO 120
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 = (D(L+1) - P) / (2.0E0 * E(L))
R = PYTHAG(G,1.0E0)
G = D(M) - P + E(L) / (G + SIGN(R,G))
S = 1.0E0
C = 1.0E0
P = 0.0E0
MML = M - L
C .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
DO 200 II = 1, MML
I = M - II
F = S * E(I)
B = C * E(I)
R = PYTHAG(F,G)
E(I+1) = R
S = F / R
C = G / R
G = D(I+1) - P
R = (D(I) - G) * S + 2.0E0 * C * B
P = S * R
D(I+1) = G + P
G = C * R - B
C .......... FORM VECTOR ..........
DO 180 K = 1, N
F = Z(K,I+1)
Z(K,I+1) = S * Z(K,I) + C * F
Z(K,I) = C * Z(K,I) - S * F
180 CONTINUE
C
200 CONTINUE
C
D(L) = D(L) - P
E(L) = G
E(M) = 0.0E0
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
DO 280 J = 1, N
P = Z(J,I)
Z(J,I) = Z(J,K)
Z(J,K) = P
280 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 trbakv(nm,n,a,e,m,z,s,s2,s3)
c
integer i,j,k,l,m,n,nm
real a(nm,n),e(n),z(nm,m),s(n),s2(n),s3(n)
c
if (m .eq. 0) go to 200
if (n .eq. 1) go to 200
c
do 140 i = 2, n ,2
l = i - 1
c
if( i .ne. n ) go to 101
if (e(n) .eq. 0.0e0) go to 140
do 11 j = 1,m
s(j) = 0.0
11 continue
call sxmpy(m,1,s,l,nm,a(n,1),nm,z)
c .......... divisor below is negative of h formed in tred1.
c double division avoids possible underflow ..........
t = (1.0 / a(n,l)) / e(n)
do 12 j=1,m
s(j) = s(j)*t
12 continue
c
do 13 j = 1,m
t = s(j)
do 14 k = 1, l
14 z(k,j) = z(k,j) + t * a(n,k)
c
13 continue
go to 140
c
101 continue
c
c form (I - b2*v*v')(I - b1*u*u')Z
c Z - b1*u*u'*Z - b2*v*v'*Z + b1*b2*(v'*u)*v*u'*Z
c s <- -b1*u'*Z
c s2 <- -b2*v'*Z
c s3 <- -b2*(v'*u)*v + u
c Z <- Z + s3*s' + v*s2'
c
c u is of length l in A(i,1...l)
c v is of length l+1 in A(i+1,1...l+1)
c
do 111 j = 1,m
s(j) = 0.0
s2(j) = 0.0
111 continue
call sxmpy(m,1,s,l,nm,a(i,1),nm,z)
call sxmpy(m,1,s2,l+1,nm,a(i+1,1),nm,z)
c .......... divisor below is negative of h formed in tred1.
c double division avoids possible underflow ..........
if (e(i) .eq. 0.0e0) then
t = 0.0
else
t = (1.0 / a(i,l)) / e(i)
endif
if (e(i+1) .eq. 0.0e0) then
t2 = 0.0
else
t2 = (1.0 / a(i+1,l+1)) / e(i+1)
endif
do 121 j=1,m
s(j) = s(j)*t
s2(j) = s2(j)*t2
121 continue
t3 = 0.0
do 122 j = 1,l
t3 = t3 + a(i,j)*a(i+1,j)
122 continue
t3 = t3*t2
do 123 j = 1,l
s3(j) = t3*a(i+1,j) + a(i,j)
123 continue
s3(l+1) = t3*a(i+1,l+1)
c
do 135 j = 1,m
t = s(j)
t2 = s2(j)
do 120 k = 1, l+1
120 z(k,j) = (z(k,j) + t * s3(k)) + t2*a(i+1,k)
c
135 continue
c
140 continue
c
200 return
end
subroutine trbk1(nm,n,a,e,m,z,s)
c
integer i,j,k,l,m,n,nm
real a(nm,n),e(n),z(nm,m),s(n)
c
if (m .eq. 0) go to 200
if (n .eq. 1) go to 200
c
do 140 i = 2, n
l = i - 1
if (e(i) .eq. 0.0e0) go to 140
c
do 111 j = 1,m
s(j) = 0.0
111 continue
call sxmpy(m,1,s,l,nm,a(i,1),nm,z)
c .......... divisor below is negative of h formed in tred1.
c double division avoids possible underflow ..........
t = (1.0 / a(i,l)) / e(i)
do 121 j=1,m
s(j) = s(j)*t
121 continue
c
do 135 j = 1,m
t = s(j)
do 120 k = 1, l
120 z(k,j) = z(k,j) + t * a(i,k)
c
135 continue
c
140 continue
c
200 return
end
REAL FUNCTION PYTHAG(A,B)
REAL A,B
C
C FINDS SQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
C
REAL P,Q,R
P = AMAX1(ABS(A),ABS(B))
Q = AMIN1(ABS(A),ABS(B))
IF (Q .EQ. 0.0E0) GO TO 20
R = (Q/P+1)**2
PYTHAG=p*sqrt(r)
return
20 PYTHAG = P
RETURN
END
.