cortb.f 4.4 KB

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