ortran.f 3.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111
  1. *DECK ORTRAN
  2. SUBROUTINE ORTRAN (NM, N, LOW, IGH, A, ORT, Z)
  3. C***BEGIN PROLOGUE ORTRAN
  4. C***PURPOSE Accumulate orthogonal similarity transformations in the
  5. C reduction of real general matrix by ORTHES.
  6. C***LIBRARY SLATEC (EISPACK)
  7. C***CATEGORY D4C4
  8. C***TYPE SINGLE PRECISION (ORTRAN-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 ORTRANS,
  14. C NUM. MATH. 16, 181-204(1970) by Peters and Wilkinson.
  15. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
  16. C
  17. C This subroutine accumulates the orthogonal similarity
  18. C transformations used in the reduction of a REAL GENERAL
  19. C matrix to upper Hessenberg form by ORTHES.
  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 LOW and IGH are two INTEGER variables determined by the
  31. C balancing subroutine BALANC. If BALANC has not been
  32. C used, set LOW=1 and IGH equal to the order of the matrix, N.
  33. C
  34. C A contains some information about the orthogonal trans-
  35. C formations used in the reduction to Hessenberg form by
  36. C ORTHES in its strict lower triangle. A is a two-dimensional
  37. C REAL array, dimensioned A(NM,IGH).
  38. C
  39. C ORT contains further information about the orthogonal trans-
  40. C formations used in the reduction by ORTHES. Only elements
  41. C LOW through IGH are used. ORT is a one-dimensional REAL
  42. C array, dimensioned ORT(IGH).
  43. C
  44. C On OUTPUT
  45. C
  46. C Z contains the transformation matrix produced in the reduction
  47. C by ORTHES to the upper Hessenberg form. Z is a two-
  48. C dimensional REAL array, dimensioned Z(NM,N).
  49. C
  50. C ORT has been used for temporary storage as is not restored.
  51. C
  52. C Questions and comments should be directed to B. S. Garbow,
  53. C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  54. C ------------------------------------------------------------------
  55. C
  56. C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
  57. C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
  58. C system Routines - EISPACK Guide, Springer-Verlag,
  59. C 1976.
  60. C***ROUTINES CALLED (NONE)
  61. C***REVISION HISTORY (YYMMDD)
  62. C 760101 DATE WRITTEN
  63. C 890831 Modified array declarations. (WRB)
  64. C 890831 REVISION DATE from Version 3.2
  65. C 891214 Prologue converted to Version 4.0 format. (BAB)
  66. C 920501 Reformatted the REFERENCES section. (WRB)
  67. C***END PROLOGUE ORTRAN
  68. C
  69. INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1
  70. REAL A(NM,*),ORT(*),Z(NM,*)
  71. REAL G
  72. C
  73. C .......... INITIALIZE Z TO IDENTITY MATRIX ..........
  74. C***FIRST EXECUTABLE STATEMENT ORTRAN
  75. DO 80 I = 1, N
  76. C
  77. DO 60 J = 1, N
  78. 60 Z(I,J) = 0.0E0
  79. C
  80. Z(I,I) = 1.0E0
  81. 80 CONTINUE
  82. C
  83. KL = IGH - LOW - 1
  84. IF (KL .LT. 1) GO TO 200
  85. C .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
  86. DO 140 MM = 1, KL
  87. MP = IGH - MM
  88. IF (A(MP,MP-1) .EQ. 0.0E0) GO TO 140
  89. MP1 = MP + 1
  90. C
  91. DO 100 I = MP1, IGH
  92. 100 ORT(I) = A(I,MP-1)
  93. C
  94. DO 130 J = MP, IGH
  95. G = 0.0E0
  96. C
  97. DO 110 I = MP, IGH
  98. 110 G = G + ORT(I) * Z(I,J)
  99. C .......... DIVISOR BELOW IS NEGATIVE OF H FORMED IN ORTHES.
  100. C DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ..........
  101. G = (G / ORT(MP)) / A(MP,MP-1)
  102. C
  103. DO 120 I = MP, IGH
  104. 120 Z(I,J) = Z(I,J) + G * ORT(I)
  105. C
  106. 130 CONTINUE
  107. C
  108. 140 CONTINUE
  109. C
  110. 200 RETURN
  111. END