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