binaries/src/disembl/Tisean_3.0.1/source_f/slatec/tred2.f

   1 *DECK TRED2
   2       SUBROUTINE TRED2 (NM, N, A, D, E, Z)
   3 C***BEGIN PROLOGUE  TRED2
   4 C***PURPOSE  Reduce a real symmetric matrix to a symmetric tridiagonal
   5 C            matrix using and accumulating orthogonal transformations.
   6 C***LIBRARY   SLATEC (EISPACK)
   7 C***CATEGORY  D4C1B1
   8 C***TYPE      SINGLE PRECISION (TRED2-S)
   9 C***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
  10 C***AUTHOR  Smith, B. T., et al.
  11 C***DESCRIPTION
  12 C
  13 C     This subroutine is a translation of the ALGOL procedure TRED2,
  14 C     NUM. MATH. 11, 181-195(1968) by Martin, Reinsch, and Wilkinson.
  15 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  16 C
  17 C     This subroutine reduces a REAL SYMMETRIC matrix to a
  18 C     symmetric tridiagonal matrix using and accumulating
  19 C     orthogonal similarity transformations.
  20 C
  21 C     On Input
  22 C
  23 C        NM must be set to the row dimension of the two-dimensional
  24 C          array parameters, A and Z, as declared in the calling
  25 C          program dimension statement.  NM is an INTEGER variable.
  26 C
  27 C        N is the order of the matrix A.  N is an INTEGER variable.
  28 C          N must be less than or equal to NM.
  29 C
  30 C        A contains the real symmetric input matrix.  Only the lower
  31 C          triangle of the matrix need be supplied.  A is a two-
  32 C          dimensional REAL array, dimensioned A(NM,N).
  33 C
  34 C     On Output
  35 C
  36 C        D contains the diagonal elements of the symmetric tridiagonal
  37 C          matrix.  D is a one-dimensional REAL array, dimensioned D(N).
  38 C
  39 C        E contains the subdiagonal elements of the symmetric
  40 C          tridiagonal matrix in its last N-1 positions.  E(1) is set
  41 C          to zero.  E is a one-dimensional REAL array, dimensioned
  42 C          E(N).
  43 C
  44 C        Z contains the orthogonal transformation matrix produced in
  45 C          the reduction.  Z is a two-dimensional REAL array,
  46 C          dimensioned Z(NM,N).
  47 C
  48 C        A and Z may coincide.  If distinct, A is unaltered.
  49 C
  50 C     Questions and comments should be directed to B. S. Garbow,
  51 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  52 C     ------------------------------------------------------------------
  53 C
  54 C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
  55 C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
  56 C                 system Routines - EISPACK Guide, Springer-Verlag,
  57 C                 1976.
  58 C***ROUTINES CALLED  (NONE)
  59 C***REVISION HISTORY  (YYMMDD)
  60 C   760101  DATE WRITTEN
  61 C   890831  Modified array declarations.  (WRB)
  62 C   890831  REVISION DATE from Version 3.2
  63 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  64 C   920501  Reformatted the REFERENCES section.  (WRB)
  65 C***END PROLOGUE  TRED2
  66 C
  67       INTEGER I,J,K,L,N,II,NM,JP1
  68       REAL A(NM,*),D(*),E(*),Z(NM,*)
  69       REAL F,G,H,HH,SCALE
  70 C
  71 C***FIRST EXECUTABLE STATEMENT  TRED2
  72       DO 100 I = 1, N
  73 C
  74          DO 100 J = 1, I
  75             Z(I,J) = A(I,J)
  76   100 CONTINUE
  77 C
  78       IF (N .EQ. 1) GO TO 320
  79 C     .......... FOR I=N STEP -1 UNTIL 2 DO -- ..........
  80       DO 300 II = 2, N
  81          I = N + 2 - II
  82          L = I - 1
  83          H = 0.0E0
  84          SCALE = 0.0E0
  85          IF (L .LT. 2) GO TO 130
  86 C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
  87          DO 120 K = 1, L
  88   120    SCALE = SCALE + ABS(Z(I,K))
  89 C
  90          IF (SCALE .NE. 0.0E0) GO TO 140
  91   130    E(I) = Z(I,L)
  92          GO TO 290
  93 C
  94   140    DO 150 K = 1, L
  95             Z(I,K) = Z(I,K) / SCALE
  96             H = H + Z(I,K) * Z(I,K)
  97   150    CONTINUE
  98 C
  99          F = Z(I,L)
 100          G = -SIGN(SQRT(H),F)
 101          E(I) = SCALE * G
 102          H = H - F * G
 103          Z(I,L) = F - G
 104          F = 0.0E0
 105 C
 106          DO 240 J = 1, L
 107             Z(J,I) = Z(I,J) / H
 108             G = 0.0E0
 109 C     .......... FORM ELEMENT OF A*U ..........
 110             DO 180 K = 1, J
 111   180       G = G + Z(J,K) * Z(I,K)
 112 C
 113             JP1 = J + 1
 114             IF (L .LT. JP1) GO TO 220
 115 C
 116             DO 200 K = JP1, L
 117   200       G = G + Z(K,J) * Z(I,K)
 118 C     .......... FORM ELEMENT OF P ..........
 119   220       E(J) = G / H
 120             F = F + E(J) * Z(I,J)
 121   240    CONTINUE
 122 C
 123          HH = F / (H + H)
 124 C     .......... FORM REDUCED A ..........
 125          DO 260 J = 1, L
 126             F = Z(I,J)
 127             G = E(J) - HH * F
 128             E(J) = G
 129 C
 130             DO 260 K = 1, J
 131                Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
 132   260    CONTINUE
 133 C
 134   290    D(I) = H
 135   300 CONTINUE
 136 C
 137   320 D(1) = 0.0E0
 138       E(1) = 0.0E0
 139 C     .......... ACCUMULATION OF TRANSFORMATION MATRICES ..........
 140       DO 500 I = 1, N
 141          L = I - 1
 142          IF (D(I) .EQ. 0.0E0) GO TO 380
 143 C
 144          DO 360 J = 1, L
 145             G = 0.0E0
 146 C
 147             DO 340 K = 1, L
 148   340       G = G + Z(I,K) * Z(K,J)
 149 C
 150             DO 360 K = 1, L
 151                Z(K,J) = Z(K,J) - G * Z(K,I)
 152   360    CONTINUE
 153 C
 154   380    D(I) = Z(I,I)
 155          Z(I,I) = 1.0E0
 156          IF (L .LT. 1) GO TO 500
 157 C
 158          DO 400 J = 1, L
 159             Z(I,J) = 0.0E0
 160             Z(J,I) = 0.0E0
 161   400    CONTINUE
 162 C
 163   500 CONTINUE
 164 C
 165       RETURN
 166       END