ortbak.f 3.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110
  1. *DECK ORTBAK
  2. SUBROUTINE ORTBAK (NM, LOW, IGH, A, ORT, M, Z)
  3. C***BEGIN PROLOGUE ORTBAK
  4. C***PURPOSE Form the eigenvectors of a general real matrix from the
  5. C eigenvectors of the upper Hessenberg matrix output from
  6. C ORTHES.
  7. C***LIBRARY SLATEC (EISPACK)
  8. C***CATEGORY D4C4
  9. C***TYPE SINGLE PRECISION (ORTBAK-S, CORTB-C)
  10. C***KEYWORDS EIGENVALUES, EIGENVECTORS, 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 ORTBAK,
  15. C NUM. MATH. 12, 349-368(1968) by Martin and Wilkinson.
  16. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
  17. C
  18. C This subroutine forms the eigenvectors of a REAL GENERAL
  19. C matrix by back transforming those of the corresponding
  20. C upper Hessenberg matrix determined by ORTHES.
  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 parameters, A and Z, as declared in the calling
  26. C program dimension statement. NM is an INTEGER variable.
  27. C
  28. C LOW and IGH are two INTEGER variables determined by the
  29. C balancing subroutine BALANC. If BALANC has not been
  30. C used, set LOW=1 and IGH equal to the order of the matrix.
  31. C
  32. C A contains some information about the orthogonal trans-
  33. C formations used in the reduction to Hessenberg form by
  34. C ORTHES in its strict lower triangle. A is a two-dimensional
  35. C REAL array, dimensioned A(NM,IGH).
  36. C
  37. C ORT contains further information about the orthogonal trans-
  38. C formations used in the reduction by ORTHES. Only elements
  39. C LOW through IGH are used. ORT is a one-dimensional REAL
  40. C array, dimensioned ORT(IGH).
  41. C
  42. C M is the number of columns of Z to be back transformed.
  43. C M is an INTEGER variable.
  44. C
  45. C Z contains the real and imaginary parts of the eigenvectors to
  46. C be back transformed in its first M columns. Z is a two-
  47. C dimensional REAL array, dimensioned Z(NM,M).
  48. C
  49. C On OUTPUT
  50. C
  51. C Z contains the real and imaginary parts of the transformed
  52. C eigenvectors in its first M columns.
  53. C
  54. C ORT has been used for temporary storage as is not restored.
  55. C
  56. C NOTE that ORTBAK preserves vector Euclidean norms.
  57. C
  58. C Questions and comments should be directed to B. S. Garbow,
  59. C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  60. C ------------------------------------------------------------------
  61. C
  62. C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
  63. C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
  64. C system Routines - EISPACK Guide, Springer-Verlag,
  65. C 1976.
  66. C***ROUTINES CALLED (NONE)
  67. C***REVISION HISTORY (YYMMDD)
  68. C 760101 DATE WRITTEN
  69. C 890831 Modified array declarations. (WRB)
  70. C 890831 REVISION DATE from Version 3.2
  71. C 891214 Prologue converted to Version 4.0 format. (BAB)
  72. C 920501 Reformatted the REFERENCES section. (WRB)
  73. C***END PROLOGUE ORTBAK
  74. C
  75. INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1
  76. REAL A(NM,*),ORT(*),Z(NM,*)
  77. REAL G
  78. C
  79. C***FIRST EXECUTABLE STATEMENT ORTBAK
  80. IF (M .EQ. 0) GO TO 200
  81. LA = IGH - 1
  82. KP1 = LOW + 1
  83. IF (LA .LT. KP1) GO TO 200
  84. C .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
  85. DO 140 MM = KP1, LA
  86. MP = LOW + IGH - MM
  87. IF (A(MP,MP-1) .EQ. 0.0E0) GO TO 140
  88. MP1 = MP + 1
  89. C
  90. DO 100 I = MP1, IGH
  91. 100 ORT(I) = A(I,MP-1)
  92. C
  93. DO 130 J = 1, M
  94. G = 0.0E0
  95. C
  96. DO 110 I = MP, IGH
  97. 110 G = G + ORT(I) * Z(I,J)
  98. C .......... DIVISOR BELOW IS NEGATIVE OF H FORMED IN ORTHES.
  99. C DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ..........
  100. G = (G / ORT(MP)) / A(MP,MP-1)
  101. C
  102. DO 120 I = MP, IGH
  103. 120 Z(I,J) = Z(I,J) + G * ORT(I)
  104. C
  105. 130 CONTINUE
  106. C
  107. 140 CONTINUE
  108. C
  109. 200 RETURN
  110. END