[CONTACT]

[ABOUT]

[POLICY]

[ADVERTISE]

real real real real do

Found at: ftp.icm.edu.pl:70/packages/netlib/misc/tred2v.f

           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 

		

		

		
.

NEW PAGES:

[ODDNUGGET]

[GOPHER]