trbak3.f 3.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107
  1. *DECK TRBAK3
  2. SUBROUTINE TRBAK3 (NM, N, NV, A, M, Z)
  3. C***BEGIN PROLOGUE TRBAK3
  4. C***PURPOSE Form the eigenvectors of a real symmetric matrix from the
  5. C eigenvectors of a symmetric tridiagonal matrix formed
  6. C by TRED3.
  7. C***LIBRARY SLATEC (EISPACK)
  8. C***CATEGORY D4C4
  9. C***TYPE SINGLE PRECISION (TRBAK3-S)
  10. C***KEYWORDS EIGENVECTORS OF A REAL SYMMETRIC MATRIX, EISPACK
  11. C***AUTHOR Smith, B. T., et al.
  12. C***DESCRIPTION
  13. C
  14. C This subroutine is a translation of the ALGOL procedure TRBAK3,
  15. C NUM. MATH. 11, 181-195(1968) by Martin, Reinsch, and Wilkinson.
  16. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  17. C
  18. C This subroutine forms the eigenvectors of a REAL SYMMETRIC
  19. C matrix by back transforming those of the corresponding
  20. C symmetric tridiagonal matrix determined by TRED3.
  21. C
  22. C On Input
  23. C
  24. C NM must be set to the row dimension of the two-dimensional
  25. C array parameter, Z, as declared in the calling program
  26. C dimension statement. NM is an INTEGER variable.
  27. C
  28. C N is the order of the matrix. N is an INTEGER variable.
  29. C N must be less than or equal to NM.
  30. C
  31. C NV is an INTEGER variable set equal to the dimension of the
  32. C array A as specified in the calling program. NV must not
  33. C be less than N*(N+1)/2.
  34. C
  35. C A contains information about the orthogonal transformations
  36. C used in the reduction by TRED3 in its first N*(N+1)/2
  37. C positions. A is a one-dimensional REAL array, dimensioned
  38. C A(NV).
  39. C
  40. C M is the number of columns of Z to be back transformed.
  41. C M is an INTEGER variable.
  42. C
  43. C Z contains the eigenvectors to be back transformed in its
  44. C first M columns. Z is a two-dimensional REAL array,
  45. C dimensioned Z(NM,M).
  46. C
  47. C On Output
  48. C
  49. C Z contains the transformed eigenvectors in its first M columns.
  50. C
  51. C Note that TRBAK3 preserves vector Euclidean norms.
  52. C
  53. C Questions and comments should be directed to b. s. Garbow,
  54. C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  55. C ------------------------------------------------------------------
  56. C
  57. C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
  58. C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
  59. C system Routines - EISPACK Guide, Springer-Verlag,
  60. C 1976.
  61. C***ROUTINES CALLED (NONE)
  62. C***REVISION HISTORY (YYMMDD)
  63. C 760101 DATE WRITTEN
  64. C 890831 Modified array declarations. (WRB)
  65. C 890831 REVISION DATE from Version 3.2
  66. C 891214 Prologue converted to Version 4.0 format. (BAB)
  67. C 920501 Reformatted the REFERENCES section. (WRB)
  68. C***END PROLOGUE TRBAK3
  69. C
  70. INTEGER I,J,K,L,M,N,IK,IZ,NM,NV
  71. REAL A(*),Z(NM,*)
  72. REAL H,S
  73. C
  74. C***FIRST EXECUTABLE STATEMENT TRBAK3
  75. IF (M .EQ. 0) GO TO 200
  76. IF (N .EQ. 1) GO TO 200
  77. C
  78. DO 140 I = 2, N
  79. L = I - 1
  80. IZ = (I * L) / 2
  81. IK = IZ + I
  82. H = A(IK)
  83. IF (H .EQ. 0.0E0) GO TO 140
  84. C
  85. DO 130 J = 1, M
  86. S = 0.0E0
  87. IK = IZ
  88. C
  89. DO 110 K = 1, L
  90. IK = IK + 1
  91. S = S + A(IK) * Z(K,J)
  92. 110 CONTINUE
  93. C .......... DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ..........
  94. S = (S / H) / H
  95. IK = IZ
  96. C
  97. DO 120 K = 1, L
  98. IK = IK + 1
  99. Z(K,J) = Z(K,J) - S * A(IK)
  100. 120 CONTINUE
  101. C
  102. 130 CONTINUE
  103. C
  104. 140 CONTINUE
  105. C
  106. 200 RETURN
  107. END