home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / single / comments next >
Text File  |  1996-09-28  |  156KB  |  4,144 lines
  1.       SUBROUTINE CDIV(AR,AI,BR,BI,CR,CI)
  2.       REAL AR,AI,BR,BI,CR,CI
  3. C
  4. C     COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI)
  5. C
  6.       SUBROUTINE CSROOT(XR,XI,YR,YI)
  7.       REAL XR,XI,YR,YI
  8. C
  9. C     (YR,YI) = COMPLEX SQRT(XR,XI)
  10. C     BRANCH CHOSEN SO THAT YR .GE. 0.0 AND SIGN(YI) .EQ. SIGN(XI)
  11. C
  12.       REAL FUNCTION EPSLON (X)
  13.       REAL X
  14. C
  15. C     ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
  16. C
  17.       REAL FUNCTION PYTHAG(A,B)
  18.       REAL A,B
  19. C
  20. C     FINDS SQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
  21. C
  22.       SUBROUTINE BAKVEC(NM,N,T,E,M,Z,IERR)
  23. C
  24.       INTEGER I,J,M,N,NM,IERR
  25.       REAL T(NM,3),E(N),Z(NM,M)
  26. C
  27. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A NONSYMMETRIC
  28. C     TRIDIAGONAL MATRIX BY BACK TRANSFORMING THOSE OF THE
  29. C     CORRESPONDING SYMMETRIC MATRIX DETERMINED BY  FIGI.
  30. C
  31. C     ON INPUT
  32. C
  33. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  34. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  35. C          DIMENSION STATEMENT.
  36. C
  37. C        N IS THE ORDER OF THE MATRIX.
  38. C
  39. C        T CONTAINS THE NONSYMMETRIC MATRIX.  ITS SUBDIAGONAL IS
  40. C          STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN,
  41. C          ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN,
  42. C          AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF
  43. C          THE THIRD COLUMN.  T(1,1) AND T(N,3) ARE ARBITRARY.
  44. C
  45. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC
  46. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  47. C
  48. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
  49. C
  50. C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
  51. C          IN ITS FIRST M COLUMNS.
  52. C
  53. C     ON OUTPUT
  54. C
  55. C        T IS UNALTERED.
  56. C
  57. C        E IS DESTROYED.
  58. C
  59. C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
  60. C          IN ITS FIRST M COLUMNS.
  61. C
  62. C        IERR IS SET TO
  63. C          ZERO       FOR NORMAL RETURN,
  64. C          2*N+I      IF E(I) IS ZERO WITH T(I,1) OR T(I-1,3) NON-ZERO.
  65. C                     IN THIS CASE, THE SYMMETRIC MATRIX IS NOT SIMILAR
  66. C                     TO THE ORIGINAL MATRIX, AND THE EIGENVECTORS
  67. C                     CANNOT BE FOUND BY THIS PROGRAM.
  68. C
  69. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  70. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  71. C
  72. C     THIS VERSION DATED AUGUST 1983.
  73. C
  74. C     ------------------------------------------------------------------
  75. C
  76.       SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE)
  77. C
  78.       INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC
  79.       REAL A(NM,N),SCALE(N)
  80.       REAL C,F,G,R,S,B2,RADIX
  81.       LOGICAL NOCONV
  82. C
  83. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE,
  84. C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
  85. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
  86. C
  87. C     THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES
  88. C     EIGENVALUES WHENEVER POSSIBLE.
  89. C
  90. C     ON INPUT
  91. C
  92. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  93. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  94. C          DIMENSION STATEMENT.
  95. C
  96. C        N IS THE ORDER OF THE MATRIX.
  97. C
  98. C        A CONTAINS THE INPUT MATRIX TO BE BALANCED.
  99. C
  100. C     ON OUTPUT
  101. C
  102. C        A CONTAINS THE BALANCED MATRIX.
  103. C
  104. C        LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J)
  105. C          IS EQUAL TO ZERO IF
  106. C           (1) I IS GREATER THAN J AND
  107. C           (2) J=1,...,LOW-1 OR I=IGH+1,...,N.
  108. C
  109. C        SCALE CONTAINS INFORMATION DETERMINING THE
  110. C           PERMUTATIONS AND SCALING FACTORS USED.
  111. C
  112. C     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH
  113. C     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED
  114. C     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS
  115. C     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN
  116. C        SCALE(J) = P(J),    FOR J = 1,...,LOW-1
  117. C                 = D(J,J),      J = LOW,...,IGH
  118. C                 = P(J)         J = IGH+1,...,N.
  119. C     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1,
  120. C     THEN 1 TO LOW-1.
  121. C
  122. C     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY.
  123. C
  124. C     THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN
  125. C     BALANC  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS
  126. C     K,L HAVE BEEN REVERSED.)
  127. C
  128. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  129. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  130. C
  131. C     THIS VERSION DATED AUGUST 1983.
  132. C
  133. C     ------------------------------------------------------------------
  134. C
  135.       SUBROUTINE BALBAK(NM,N,LOW,IGH,SCALE,M,Z)
  136. C
  137.       INTEGER I,J,K,M,N,II,NM,IGH,LOW
  138.       REAL SCALE(N),Z(NM,M)
  139.       REAL S
  140. C
  141. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK,
  142. C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
  143. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
  144. C
  145. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL
  146. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  147. C     BALANCED MATRIX DETERMINED BY  BALANC.
  148. C
  149. C     ON INPUT
  150. C
  151. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  152. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  153. C          DIMENSION STATEMENT.
  154. C
  155. C        N IS THE ORDER OF THE MATRIX.
  156. C
  157. C        LOW AND IGH ARE INTEGERS DETERMINED BY  BALANC.
  158. C
  159. C        SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS
  160. C          AND SCALING FACTORS USED BY  BALANC.
  161. C
  162. C        M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED.
  163. C
  164. C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN-
  165. C          VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS.
  166. C
  167. C     ON OUTPUT
  168. C
  169. C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE
  170. C          TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS.
  171. C
  172. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  173. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  174. C
  175. C     THIS VERSION DATED AUGUST 1983.
  176. C
  177. C     ------------------------------------------------------------------
  178. C
  179.       SUBROUTINE BANDR(NM,N,MB,A,D,E,E2,MATZ,Z)
  180. C
  181.       INTEGER J,K,L,N,R,I1,I2,J1,J2,KR,MB,MR,M1,NM,N2,R1,UGL,MAXL,MAXR
  182.       REAL A(NM,MB),D(N),E(N),E2(N),Z(NM,N)
  183.       REAL G,U,B1,B2,C2,F1,F2,S2,DMIN,DMINRT
  184.       LOGICAL MATZ
  185. C
  186. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BANDRD,
  187. C     NUM. MATH. 12, 231-241(1968) BY SCHWARZ.
  188. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 273-283(1971).
  189. C
  190. C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC BAND MATRIX
  191. C     TO A SYMMETRIC TRIDIAGONAL MATRIX USING AND OPTIONALLY
  192. C     ACCUMULATING ORTHOGONAL SIMILARITY TRANSFORMATIONS.
  193. C
  194. C     ON INPUT
  195. C
  196. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  197. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  198. C          DIMENSION STATEMENT.
  199. C
  200. C        N IS THE ORDER OF THE MATRIX.
  201. C
  202. C        MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE
  203. C          NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL
  204. C          DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE
  205. C          LOWER TRIANGLE OF THE MATRIX.
  206. C
  207. C        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT
  208. C          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL
  209. C          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN,
  210. C          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE
  211. C          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY
  212. C          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
  213. C          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY.
  214. C
  215. C        MATZ SHOULD BE SET TO .TRUE. IF THE TRANSFORMATION MATRIX IS
  216. C          TO BE ACCUMULATED, AND TO .FALSE. OTHERWISE.
  217. C
  218. C     ON OUTPUT
  219. C
  220. C        A HAS BEEN DESTROYED, EXCEPT FOR ITS LAST TWO COLUMNS WHICH
  221. C          CONTAIN A COPY OF THE TRIDIAGONAL MATRIX.
  222. C
  223. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
  224. C
  225. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  226. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
  227. C
  228. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  229. C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
  230. C
  231. C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX PRODUCED IN
  232. C          THE REDUCTION IF MATZ HAS BEEN SET TO .TRUE.  OTHERWISE, Z
  233. C          IS NOT REFERENCED.
  234. C
  235. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  236. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  237. C
  238. C     THIS VERSION DATED AUGUST 1983.
  239. C
  240. C     ------------------------------------------------------------------
  241. C
  242.       SUBROUTINE BANDV(NM,N,MBW,A,E21,M,W,Z,IERR,NV,RV,RV6)
  243. C
  244.       INTEGER I,J,K,M,N,R,II,IJ,JJ,KJ,MB,M1,NM,NV,IJ1,ITS,KJ1,MBW,M21,
  245.      X        IERR,MAXJ,MAXK,GROUP
  246.       REAL A(NM,MBW),W(M),Z(NM,M),RV(NV),RV6(N)
  247.       REAL U,V,UK,XU,X0,X1,E21,EPS2,EPS3,EPS4,NORM,ORDER,
  248.      X       EPSLON,PYTHAG
  249. C
  250. C     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL SYMMETRIC
  251. C     BAND MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, USING INVERSE
  252. C     ITERATION.  THE SUBROUTINE MAY ALSO BE USED TO SOLVE SYSTEMS
  253. C     OF LINEAR EQUATIONS WITH A SYMMETRIC OR NON-SYMMETRIC BAND
  254. C     COEFFICIENT MATRIX.
  255. C
  256. C     ON INPUT
  257. C
  258. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  259. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  260. C          DIMENSION STATEMENT.
  261. C
  262. C        N IS THE ORDER OF THE MATRIX.
  263. C
  264. C        MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE
  265. C          BAND MATRIX.  IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF)
  266. C          BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT
  267. C          DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO
  268. C          SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE
  269. C          MATRIX.  IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS
  270. C          OF LINEAR EQUATIONS AND THE COEFFICIENT MATRIX IS NOT
  271. C          SYMMETRIC, IT MUST HOWEVER HAVE THE SAME NUMBER OF ADJACENT
  272. C          DIAGONALS ABOVE THE MAIN DIAGONAL AS BELOW, AND IN THIS
  273. C          CASE, MBW=2*MB-1.
  274. C
  275. C        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT
  276. C          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL
  277. C          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN,
  278. C          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE
  279. C          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY
  280. C          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF COLUMN MB.
  281. C          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR
  282. C          EQUATIONS AND THE COEFFICIENT MATRIX IS NOT SYMMETRIC, A IS
  283. C          N BY 2*MB-1 INSTEAD WITH LOWER TRIANGLE AS ABOVE AND WITH
  284. C          ITS FIRST SUPERDIAGONAL STORED IN THE FIRST N-1 POSITIONS OF
  285. C          COLUMN MB+1, ITS SECOND SUPERDIAGONAL IN THE FIRST N-2
  286. C          POSITIONS OF COLUMN MB+2, FURTHER SUPERDIAGONALS SIMILARLY,
  287. C          AND FINALLY ITS HIGHEST SUPERDIAGONAL IN THE FIRST N+1-MB
  288. C          POSITIONS OF THE LAST COLUMN.
  289. C          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY.
  290. C
  291. C        E21 SPECIFIES THE ORDERING OF THE EIGENVALUES AND CONTAINS
  292. C            0.0E0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR
  293. C            2.0E0 IF THE EIGENVALUES ARE IN DESCENDING ORDER.
  294. C          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR
  295. C          EQUATIONS, E21 SHOULD BE SET TO 1.0E0 IF THE COEFFICIENT
  296. C          MATRIX IS SYMMETRIC AND TO -1.0E0 IF NOT.
  297. C
  298. C        M IS THE NUMBER OF SPECIFIED EIGENVALUES OR THE NUMBER OF
  299. C          SYSTEMS OF LINEAR EQUATIONS.
  300. C
  301. C        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
  302. C          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR
  303. C          EQUATIONS (A-W(R)*I)*X(R)=B(R), WHERE I IS THE IDENTITY
  304. C          MATRIX, W(R) SHOULD BE SET ACCORDINGLY, FOR R=1,2,...,M.
  305. C
  306. C        Z CONTAINS THE CONSTANT MATRIX COLUMNS (B(R),R=1,2,...,M), IF
  307. C          THE SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS.
  308. C
  309. C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV
  310. C          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  311. C
  312. C     ON OUTPUT
  313. C
  314. C        A AND W ARE UNALTERED.
  315. C
  316. C        Z CONTAINS THE ASSOCIATED SET OF ORTHOGONAL EIGENVECTORS.
  317. C          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO.  IF THE
  318. C          SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS,
  319. C          Z CONTAINS THE SOLUTION MATRIX COLUMNS (X(R),R=1,2,...,M).
  320. C
  321. C        IERR IS SET TO
  322. C          ZERO       FOR NORMAL RETURN,
  323. C          -R         IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
  324. C                     EIGENVALUE FAILS TO CONVERGE, OR IF THE R-TH
  325. C                     SYSTEM OF LINEAR EQUATIONS IS NEARLY SINGULAR.
  326. C
  327. C        RV AND RV6 ARE TEMPORARY STORAGE ARRAYS.  NOTE THAT RV IS
  328. C          OF DIMENSION AT LEAST N*(2*MB-1).  IF THE SUBROUTINE
  329. C          IS BEING USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, THE
  330. C          DETERMINANT (UP TO SIGN) OF A-W(M)*I IS AVAILABLE, UPON
  331. C          RETURN, AS THE PRODUCT OF THE FIRST N ELEMENTS OF RV.
  332. C
  333. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  334. C
  335. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  336. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  337. C
  338. C     THIS VERSION DATED AUGUST 1983.
  339. C
  340. C     ------------------------------------------------------------------
  341. C
  342.       SUBROUTINE BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,IERR,RV4,RV5)
  343. C
  344.       INTEGER I,J,K,L,M,N,P,Q,R,S,II,MM,M1,M2,TAG,IERR,ISTURM
  345.       REAL D(N),E(N),E2(N),W(MM),RV4(N),RV5(N)
  346.       REAL U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,TST1,TST2,EPSLON
  347.       INTEGER IND(MM)
  348. C
  349. C     THIS SUBROUTINE IS A TRANSLATION OF THE BISECTION TECHNIQUE
  350. C     IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
  351. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
  352. C
  353. C     THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL
  354. C     SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL,
  355. C     USING BISECTION.
  356. C
  357. C     ON INPUT
  358. C
  359. C        N IS THE ORDER OF THE MATRIX.
  360. C
  361. C        EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED
  362. C          EIGENVALUES.  IF THE INPUT EPS1 IS NON-POSITIVE,
  363. C          IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE,
  364. C          NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE
  365. C          PRECISION AND THE 1-NORM OF THE SUBMATRIX.
  366. C
  367. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  368. C
  369. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  370. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  371. C
  372. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  373. C          E2(1) IS ARBITRARY.
  374. C
  375. C        LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES.
  376. C          IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND.
  377. C
  378. C        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF
  379. C          EIGENVALUES IN THE INTERVAL.  WARNING. IF MORE THAN
  380. C          MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL,
  381. C          AN ERROR RETURN IS MADE WITH NO EIGENVALUES FOUND.
  382. C
  383. C     ON OUTPUT
  384. C
  385. C        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS
  386. C          (LAST) DEFAULT VALUE.
  387. C
  388. C        D AND E ARE UNALTERED.
  389. C
  390. C        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
  391. C          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
  392. C          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
  393. C          E2(1) IS ALSO SET TO ZERO.
  394. C
  395. C        M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB).
  396. C
  397. C        W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER.
  398. C
  399. C        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
  400. C          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
  401. C          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
  402. C          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC..
  403. C
  404. C        IERR IS SET TO
  405. C          ZERO       FOR NORMAL RETURN,
  406. C          3*N+1      IF M EXCEEDS MM.
  407. C
  408. C        RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS.
  409. C
  410. C     THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM
  411. C     APPEARS IN BISECT IN-LINE.
  412. C
  413. C     NOTE THAT SUBROUTINE TQL1 OR IMTQL1 IS GENERALLY FASTER THAN
  414. C     BISECT, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND.
  415. C
  416. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  417. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  418. C
  419. C     THIS VERSION DATED AUGUST 1983.
  420. C
  421. C     ------------------------------------------------------------------
  422. C
  423.       SUBROUTINE BQR(NM,N,MB,A,T,R,IERR,NV,RV)
  424. C
  425.       INTEGER I,J,K,L,M,N,II,IK,JK,JM,KJ,KK,KM,LL,MB,MK,MN,MZ,
  426.      X        M1,M2,M3,M4,NI,NM,NV,ITS,KJ1,M21,M31,IERR,IMULT
  427.       REAL A(NM,MB),RV(NV)
  428.       REAL F,G,Q,R,S,T,TST1,TST2,SCALE,PYTHAG
  429. C
  430. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BQR,
  431. C     NUM. MATH. 16, 85-92(1970) BY MARTIN, REINSCH, AND WILKINSON.
  432. C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971).
  433. C
  434. C     THIS SUBROUTINE FINDS THE EIGENVALUE OF SMALLEST (USUALLY)
  435. C     MAGNITUDE OF A REAL SYMMETRIC BAND MATRIX USING THE
  436. C     QR ALGORITHM WITH SHIFTS OF ORIGIN.  CONSECUTIVE CALLS
  437. C     CAN BE MADE TO FIND FURTHER EIGENVALUES.
  438. C
  439. C     ON INPUT
  440. C
  441. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  442. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  443. C          DIMENSION STATEMENT.
  444. C
  445. C        N IS THE ORDER OF THE MATRIX.
  446. C
  447. C        MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE
  448. C          NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL
  449. C          DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE
  450. C          LOWER TRIANGLE OF THE MATRIX.
  451. C
  452. C        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT
  453. C          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL
  454. C          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN,
  455. C          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE
  456. C          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY
  457. C          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
  458. C          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY.
  459. C          ON A SUBSEQUENT CALL, ITS OUTPUT CONTENTS FROM THE PREVIOUS
  460. C          CALL SHOULD BE PASSED.
  461. C
  462. C        T SPECIFIES THE SHIFT (OF EIGENVALUES) APPLIED TO THE DIAGONAL
  463. C          OF A IN FORMING THE INPUT MATRIX. WHAT IS ACTUALLY DETERMINED
  464. C          IS THE EIGENVALUE OF A+TI (I IS THE IDENTITY MATRIX) NEAREST
  465. C          TO T.  ON A SUBSEQUENT CALL, THE OUTPUT VALUE OF T FROM THE
  466. C          PREVIOUS CALL SHOULD BE PASSED IF THE NEXT NEAREST EIGENVALUE
  467. C          IS SOUGHT.
  468. C
  469. C        R SHOULD BE SPECIFIED AS ZERO ON THE FIRST CALL, AND AS ITS
  470. C          OUTPUT VALUE FROM THE PREVIOUS CALL ON A SUBSEQUENT CALL.
  471. C          IT IS USED TO DETERMINE WHEN THE LAST ROW AND COLUMN OF
  472. C          THE TRANSFORMED BAND MATRIX CAN BE REGARDED AS NEGLIGIBLE.
  473. C
  474. C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV
  475. C          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  476. C
  477. C     ON OUTPUT
  478. C
  479. C        A CONTAINS THE TRANSFORMED BAND MATRIX.  THE MATRIX A+TI
  480. C          DERIVED FROM THE OUTPUT PARAMETERS IS SIMILAR TO THE
  481. C          INPUT A+TI TO WITHIN ROUNDING ERRORS.  ITS LAST ROW AND
  482. C          COLUMN ARE NULL (IF IERR IS ZERO).
  483. C
  484. C        T CONTAINS THE COMPUTED EIGENVALUE OF A+TI (IF IERR IS ZERO).
  485. C
  486. C        R CONTAINS THE MAXIMUM OF ITS INPUT VALUE AND THE NORM OF THE
  487. C          LAST COLUMN OF THE INPUT MATRIX A.
  488. C
  489. C        IERR IS SET TO
  490. C          ZERO       FOR NORMAL RETURN,
  491. C          N          IF THE EIGENVALUE HAS NOT BEEN
  492. C                     DETERMINED AFTER 30 ITERATIONS.
  493. C
  494. C        RV IS A TEMPORARY STORAGE ARRAY OF DIMENSION AT LEAST
  495. C          (2*MB**2+4*MB-3).  THE FIRST (3*MB-2) LOCATIONS CORRESPOND
  496. C          TO THE ALGOL ARRAY B, THE NEXT (2*MB-1) LOCATIONS CORRESPOND
  497. C          TO THE ALGOL ARRAY H, AND THE FINAL (2*MB**2-MB) LOCATIONS
  498. C          CORRESPOND TO THE MB BY (2*MB-1) ALGOL ARRAY U.
  499. C
  500. C     NOTE. FOR A SUBSEQUENT CALL, N SHOULD BE REPLACED BY N-1, BUT
  501. C     MB SHOULD NOT BE ALTERED EVEN WHEN IT EXCEEDS THE CURRENT N.
  502. C
  503. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  504. C
  505. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  506. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  507. C
  508. C     THIS VERSION DATED AUGUST 1983.
  509. C
  510. C     ------------------------------------------------------------------
  511. C
  512.       SUBROUTINE CBABK2(NM,N,LOW,IGH,SCALE,M,ZR,ZI)
  513. C
  514.       INTEGER I,J,K,M,N,II,NM,IGH,LOW
  515.       REAL SCALE(N),ZR(NM,M),ZI(NM,M)
  516.       REAL S
  517. C
  518. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE
  519. C     CBABK2, WHICH IS A COMPLEX VERSION OF BALBAK,
  520. C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
  521. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
  522. C
  523. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL
  524. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  525. C     BALANCED MATRIX DETERMINED BY  CBAL.
  526. C
  527. C     ON INPUT
  528. C
  529. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  530. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  531. C          DIMENSION STATEMENT.
  532. C
  533. C        N IS THE ORDER OF THE MATRIX.
  534. C
  535. C        LOW AND IGH ARE INTEGERS DETERMINED BY  CBAL.
  536. C
  537. C        SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS
  538. C          AND SCALING FACTORS USED BY  CBAL.
  539. C
  540. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
  541. C
  542. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  543. C          RESPECTIVELY, OF THE EIGENVECTORS TO BE
  544. C          BACK TRANSFORMED IN THEIR FIRST M COLUMNS.
  545. C
  546. C     ON OUTPUT
  547. C
  548. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  549. C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
  550. C          IN THEIR FIRST M COLUMNS.
  551. C
  552. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  553. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  554. C
  555. C     THIS VERSION DATED AUGUST 1983.
  556. C
  557. C     ------------------------------------------------------------------
  558. C
  559.       SUBROUTINE CBAL(NM,N,AR,AI,LOW,IGH,SCALE)
  560. C
  561.       INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC
  562.       REAL AR(NM,N),AI(NM,N),SCALE(N)
  563.       REAL C,F,G,R,S,B2,RADIX
  564.       LOGICAL NOCONV
  565. C
  566. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE
  567. C     CBALANCE, WHICH IS A COMPLEX VERSION OF BALANCE,
  568. C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
  569. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
  570. C
  571. C     THIS SUBROUTINE BALANCES A COMPLEX MATRIX AND ISOLATES
  572. C     EIGENVALUES WHENEVER POSSIBLE.
  573. C
  574. C     ON INPUT
  575. C
  576. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  577. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  578. C          DIMENSION STATEMENT.
  579. C
  580. C        N IS THE ORDER OF THE MATRIX.
  581. C
  582. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
  583. C          RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED.
  584. C
  585. C     ON OUTPUT
  586. C
  587. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
  588. C          RESPECTIVELY, OF THE BALANCED MATRIX.
  589. C
  590. C        LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J)
  591. C          ARE EQUAL TO ZERO IF
  592. C           (1) I IS GREATER THAN J AND
  593. C           (2) J=1,...,LOW-1 OR I=IGH+1,...,N.
  594. C
  595. C        SCALE CONTAINS INFORMATION DETERMINING THE
  596. C           PERMUTATIONS AND SCALING FACTORS USED.
  597. C
  598. C     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH
  599. C     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED
  600. C     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS
  601. C     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN
  602. C        SCALE(J) = P(J),    FOR J = 1,...,LOW-1
  603. C                 = D(J,J)       J = LOW,...,IGH
  604. C                 = P(J)         J = IGH+1,...,N.
  605. C     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1,
  606. C     THEN 1 TO LOW-1.
  607. C
  608. C     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY.
  609. C
  610. C     THE ALGOL PROCEDURE EXC CONTAINED IN CBALANCE APPEARS IN
  611. C     CBAL  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS
  612. C     K,L HAVE BEEN REVERSED.)
  613. C
  614. C     ARITHMETIC IS REAL THROUGHOUT.
  615. C
  616. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  617. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  618. C
  619. C     THIS VERSION DATED AUGUST 1983.
  620. C
  621. C     ------------------------------------------------------------------
  622. C
  623.       SUBROUTINE CG(NM,N,AR,AI,WR,WI,MATZ,ZR,ZI,FV1,FV2,FV3,IERR)
  624. C
  625.       INTEGER N,NM,IS1,IS2,IERR,MATZ
  626.       REAL AR(NM,N),AI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N),
  627.      X       FV1(N),FV2(N),FV3(N)
  628. C
  629. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  630. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  631. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  632. C     OF A COMPLEX GENERAL MATRIX.
  633. C
  634. C     ON INPUT
  635. C
  636. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  637. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  638. C        DIMENSION STATEMENT.
  639. C
  640. C        N  IS THE ORDER OF THE MATRIX  A=(AR,AI).
  641. C
  642. C        AR  AND  AI  CONTAIN THE REAL AND IMAGINARY PARTS,
  643. C        RESPECTIVELY, OF THE COMPLEX GENERAL MATRIX.
  644. C
  645. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  646. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  647. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  648. C
  649. C     ON OUTPUT
  650. C
  651. C        WR  AND  WI  CONTAIN THE REAL AND IMAGINARY PARTS,
  652. C        RESPECTIVELY, OF THE EIGENVALUES.
  653. C
  654. C        ZR  AND  ZI  CONTAIN THE REAL AND IMAGINARY PARTS,
  655. C        RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO.
  656. C
  657. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  658. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR COMQR
  659. C           AND COMQR2.  THE NORMAL COMPLETION CODE IS ZERO.
  660. C
  661. C        FV1, FV2, AND  FV3  ARE TEMPORARY STORAGE ARRAYS.
  662. C
  663. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  664. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  665. C
  666. C     THIS VERSION DATED AUGUST 1983.
  667. C
  668. C     ------------------------------------------------------------------
  669. C
  670.       SUBROUTINE CH(NM,N,AR,AI,W,MATZ,ZR,ZI,FV1,FV2,FM1,IERR)
  671. C
  672.       INTEGER I,J,N,NM,IERR,MATZ
  673.       REAL AR(NM,N),AI(NM,N),W(N),ZR(NM,N),ZI(NM,N),
  674.      X       FV1(N),FV2(N),FM1(2,N)
  675. C
  676. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  677. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  678. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  679. C     OF A COMPLEX HERMITIAN MATRIX.
  680. C
  681. C     ON INPUT
  682. C
  683. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  684. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  685. C        DIMENSION STATEMENT.
  686. C
  687. C        N  IS THE ORDER OF THE MATRIX  A=(AR,AI).
  688. C
  689. C        AR  AND  AI  CONTAIN THE REAL AND IMAGINARY PARTS,
  690. C        RESPECTIVELY, OF THE COMPLEX HERMITIAN MATRIX.
  691. C
  692. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  693. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  694. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  695. C
  696. C     ON OUTPUT
  697. C
  698. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
  699. C
  700. C        ZR  AND  ZI  CONTAIN THE REAL AND IMAGINARY PARTS,
  701. C        RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO.
  702. C
  703. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  704. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT
  705. C           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO.
  706. C
  707. C        FV1, FV2, AND  FM1  ARE TEMPORARY STORAGE ARRAYS.
  708. C
  709. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  710. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  711. C
  712. C     THIS VERSION DATED AUGUST 1983.
  713. C
  714. C     ------------------------------------------------------------------
  715. C
  716.       SUBROUTINE CINVIT(NM,N,AR,AI,WR,WI,SELECT,MM,M,ZR,ZI,
  717.      X                  IERR,RM1,RM2,RV1,RV2)
  718. C
  719.       INTEGER I,J,K,M,N,S,II,MM,MP,NM,UK,IP1,ITS,KM1,IERR
  720.       REAL AR(NM,N),AI(NM,N),WR(N),WI(N),ZR(NM,MM),
  721.      X       ZI(NM,MM),RM1(N,N),RM2(N,N),RV1(N),RV2(N)
  722.       REAL X,Y,EPS3,NORM,NORMV,EPSLON,GROWTO,ILAMBD,PYTHAG,
  723.      X       RLAMBD,UKROOT
  724.       LOGICAL SELECT(N)
  725. C
  726. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE CX INVIT
  727. C     BY PETERS AND WILKINSON.
  728. C     HANDBOOK FOR AUTO. COMP. VOL.II-LINEAR ALGEBRA, 418-439(1971).
  729. C
  730. C     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A COMPLEX UPPER
  731. C     HESSENBERG MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
  732. C     USING INVERSE ITERATION.
  733. C
  734. C     ON INPUT
  735. C
  736. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  737. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  738. C          DIMENSION STATEMENT.
  739. C
  740. C        N IS THE ORDER OF THE MATRIX.
  741. C
  742. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
  743. C          RESPECTIVELY, OF THE HESSENBERG MATRIX.
  744. C
  745. C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY,
  746. C          OF THE EIGENVALUES OF THE MATRIX.  THE EIGENVALUES MUST BE
  747. C          STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE  COMLR,
  748. C          WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX.
  749. C
  750. C        SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND.  THE
  751. C          EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS
  752. C          SPECIFIED BY SETTING SELECT(J) TO .TRUE..
  753. C
  754. C        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF
  755. C          EIGENVECTORS TO BE FOUND.
  756. C
  757. C     ON OUTPUT
  758. C
  759. C        AR, AI, WI, AND SELECT ARE UNALTERED.
  760. C
  761. C        WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED
  762. C          SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS.
  763. C
  764. C        M IS THE NUMBER OF EIGENVECTORS ACTUALLY FOUND.
  765. C
  766. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY,
  767. C          OF THE EIGENVECTORS.  THE EIGENVECTORS ARE NORMALIZED
  768. C          SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1.
  769. C          ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO.
  770. C
  771. C        IERR IS SET TO
  772. C          ZERO       FOR NORMAL RETURN,
  773. C          -(2*N+1)   IF MORE THAN MM EIGENVECTORS HAVE BEEN SPECIFIED,
  774. C          -K         IF THE ITERATION CORRESPONDING TO THE K-TH
  775. C                     VALUE FAILS,
  776. C          -(N+K)     IF BOTH ERROR SITUATIONS OCCUR.
  777. C
  778. C        RM1, RM2, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS.
  779. C
  780. C     THE ALGOL PROCEDURE GUESSVEC APPEARS IN CINVIT IN LINE.
  781. C
  782. C     CALLS CDIV FOR COMPLEX DIVISION.
  783. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  784. C
  785. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  786. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  787. C
  788. C     THIS VERSION DATED AUGUST 1983.
  789. C
  790. C     ------------------------------------------------------------------
  791. C
  792.       SUBROUTINE COMBAK(NM,LOW,IGH,AR,AI,INT,M,ZR,ZI)
  793. C
  794.       INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1
  795.       REAL AR(NM,IGH),AI(NM,IGH),ZR(NM,M),ZI(NM,M)
  796.       REAL XR,XI
  797.       INTEGER INT(IGH)
  798. C
  799. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMBAK,
  800. C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
  801. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
  802. C
  803. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL
  804. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  805. C     UPPER HESSENBERG MATRIX DETERMINED BY  COMHES.
  806. C
  807. C     ON INPUT
  808. C
  809. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  810. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  811. C          DIMENSION STATEMENT.
  812. C
  813. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  814. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
  815. C          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX.
  816. C
  817. C        AR AND AI CONTAIN THE MULTIPLIERS WHICH WERE USED IN THE
  818. C          REDUCTION BY  COMHES  IN THEIR LOWER TRIANGLES
  819. C          BELOW THE SUBDIAGONAL.
  820. C
  821. C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
  822. C          INTERCHANGED IN THE REDUCTION BY  COMHES.
  823. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  824. C
  825. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
  826. C
  827. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  828. C          RESPECTIVELY, OF THE EIGENVECTORS TO BE
  829. C          BACK TRANSFORMED IN THEIR FIRST M COLUMNS.
  830. C
  831. C     ON OUTPUT
  832. C
  833. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  834. C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
  835. C          IN THEIR FIRST M COLUMNS.
  836. C
  837. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  838. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  839. C
  840. C     THIS VERSION DATED AUGUST 1983.
  841. C
  842. C     ------------------------------------------------------------------
  843. C
  844.       SUBROUTINE COMHES(NM,N,LOW,IGH,AR,AI,INT)
  845. C
  846.       INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1
  847.       REAL AR(NM,N),AI(NM,N)
  848.       REAL XR,XI,YR,YI
  849.       INTEGER INT(IGH)
  850. C
  851. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMHES,
  852. C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
  853. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
  854. C
  855. C     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE
  856. C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
  857. C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
  858. C     STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS.
  859. C
  860. C     ON INPUT
  861. C
  862. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  863. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  864. C          DIMENSION STATEMENT.
  865. C
  866. C        N IS THE ORDER OF THE MATRIX.
  867. C
  868. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  869. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
  870. C          SET LOW=1, IGH=N.
  871. C
  872. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
  873. C          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX.
  874. C
  875. C     ON OUTPUT
  876. C
  877. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
  878. C          RESPECTIVELY, OF THE HESSENBERG MATRIX.  THE
  879. C          MULTIPLIERS WHICH WERE USED IN THE REDUCTION
  880. C          ARE STORED IN THE REMAINING TRIANGLES UNDER THE
  881. C          HESSENBERG MATRIX.
  882. C
  883. C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
  884. C          INTERCHANGED IN THE REDUCTION.
  885. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  886. C
  887. C     CALLS CDIV FOR COMPLEX DIVISION.
  888. C
  889. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  890. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  891. C
  892. C     THIS VERSION DATED AUGUST 1983.
  893. C
  894. C     ------------------------------------------------------------------
  895. C
  896.       SUBROUTINE COMLR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR)
  897. C
  898.       INTEGER I,J,L,M,N,EN,LL,MM,NM,IGH,IM1,ITN,ITS,LOW,MP1,ENM1,IERR
  899.       REAL HR(NM,N),HI(NM,N),WR(N),WI(N)
  900.       REAL SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,TST1,TST2
  901. C
  902. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR,
  903. C     NUM. MATH. 12, 369-376(1968) BY MARTIN AND WILKINSON.
  904. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971).
  905. C
  906. C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX
  907. C     UPPER HESSENBERG MATRIX BY THE MODIFIED LR METHOD.
  908. C
  909. C     ON INPUT
  910. C
  911. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  912. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  913. C          DIMENSION STATEMENT.
  914. C
  915. C        N IS THE ORDER OF THE MATRIX.
  916. C
  917. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  918. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
  919. C          SET LOW=1, IGH=N.
  920. C
  921. C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
  922. C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
  923. C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE
  924. C          MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY  COMHES,
  925. C          IF PERFORMED.
  926. C
  927. C     ON OUTPUT
  928. C
  929. C        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN
  930. C          DESTROYED.  THEREFORE, THEY MUST BE SAVED BEFORE
  931. C          CALLING  COMLR  IF SUBSEQUENT CALCULATION OF
  932. C          EIGENVECTORS IS TO BE PERFORMED.
  933. C
  934. C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
  935. C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
  936. C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
  937. C          FOR INDICES IERR+1,...,N.
  938. C
  939. C        IERR IS SET TO
  940. C          ZERO       FOR NORMAL RETURN,
  941. C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
  942. C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
  943. C
  944. C     CALLS CDIV FOR COMPLEX DIVISION.
  945. C     CALLS CSROOT FOR COMPLEX SQUARE ROOT.
  946. C
  947. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  948. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  949. C
  950. C     THIS VERSION DATED AUGUST 1983.
  951. C
  952. C     ------------------------------------------------------------------
  953. C
  954.       SUBROUTINE COMLR2(NM,N,LOW,IGH,INT,HR,HI,WR,WI,ZR,ZI,IERR)
  955. C
  956.       INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NM,NN,IGH,IM1,IP1,
  957.      X        ITN,ITS,LOW,MP1,ENM1,IEND,IERR
  958.       REAL HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N)
  959.       REAL SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2
  960.       INTEGER INT(IGH)
  961. C
  962. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR2,
  963. C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
  964. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
  965. C
  966. C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
  967. C     OF A COMPLEX UPPER HESSENBERG MATRIX BY THE MODIFIED LR
  968. C     METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX
  969. C     CAN ALSO BE FOUND IF  COMHES  HAS BEEN USED TO REDUCE
  970. C     THIS GENERAL MATRIX TO HESSENBERG FORM.
  971. C
  972. C     ON INPUT
  973. C
  974. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  975. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  976. C          DIMENSION STATEMENT.
  977. C
  978. C        N IS THE ORDER OF THE MATRIX.
  979. C
  980. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  981. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
  982. C          SET LOW=1, IGH=N.
  983. C
  984. C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS INTERCHANGED
  985. C          IN THE REDUCTION BY  COMHES, IF PERFORMED.  ONLY ELEMENTS
  986. C          LOW THROUGH IGH ARE USED.  IF THE EIGENVECTORS OF THE HESSEN-
  987. C          BERG MATRIX ARE DESIRED, SET INT(J)=J FOR THESE ELEMENTS.
  988. C
  989. C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
  990. C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
  991. C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE
  992. C          MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY  COMHES,
  993. C          IF PERFORMED.  IF THE EIGENVECTORS OF THE HESSENBERG
  994. C          MATRIX ARE DESIRED, THESE ELEMENTS MUST BE SET TO ZERO.
  995. C
  996. C     ON OUTPUT
  997. C
  998. C        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN
  999. C          DESTROYED, BUT THE LOCATION HR(1,1) CONTAINS THE NORM
  1000. C          OF THE TRIANGULARIZED MATRIX.
  1001. C
  1002. C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
  1003. C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
  1004. C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
  1005. C          FOR INDICES IERR+1,...,N.
  1006. C
  1007. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  1008. C          RESPECTIVELY, OF THE EIGENVECTORS.  THE EIGENVECTORS
  1009. C          ARE UNNORMALIZED.  IF AN ERROR EXIT IS MADE, NONE OF
  1010. C          THE EIGENVECTORS HAS BEEN FOUND.
  1011. C
  1012. C        IERR IS SET TO
  1013. C          ZERO       FOR NORMAL RETURN,
  1014. C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
  1015. C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
  1016. C
  1017. C
  1018. C     CALLS CDIV FOR COMPLEX DIVISION.
  1019. C     CALLS CSROOT FOR COMPLEX SQUARE ROOT.
  1020. C
  1021. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1022. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1023. C
  1024. C     THIS VERSION DATED AUGUST 1983.
  1025. C
  1026. C     ------------------------------------------------------------------
  1027. C
  1028.       SUBROUTINE COMQR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR)
  1029. C
  1030.       INTEGER I,J,L,N,EN,LL,NM,IGH,ITN,ITS,LOW,LP1,ENM1,IERR
  1031.       REAL HR(NM,N),HI(NM,N),WR(N),WI(N)
  1032.       REAL SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2,
  1033.      X       PYTHAG
  1034. C
  1035. C     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
  1036. C     ALGOL PROCEDURE  COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN
  1037. C     AND WILKINSON.
  1038. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971).
  1039. C     THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
  1040. C     (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.
  1041. C
  1042. C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX
  1043. C     UPPER HESSENBERG MATRIX BY THE QR METHOD.
  1044. C
  1045. C     ON INPUT
  1046. C
  1047. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1048. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1049. C          DIMENSION STATEMENT.
  1050. C
  1051. C        N IS THE ORDER OF THE MATRIX.
  1052. C
  1053. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  1054. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
  1055. C          SET LOW=1, IGH=N.
  1056. C
  1057. C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
  1058. C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
  1059. C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN
  1060. C          INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN
  1061. C          THE REDUCTION BY  CORTH, IF PERFORMED.
  1062. C
  1063. C     ON OUTPUT
  1064. C
  1065. C        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN
  1066. C          DESTROYED.  THEREFORE, THEY MUST BE SAVED BEFORE
  1067. C          CALLING  COMQR  IF SUBSEQUENT CALCULATION OF
  1068. C          EIGENVECTORS IS TO BE PERFORMED.
  1069. C
  1070. C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
  1071. C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
  1072. C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
  1073. C          FOR INDICES IERR+1,...,N.
  1074. C
  1075. C        IERR IS SET TO
  1076. C          ZERO       FOR NORMAL RETURN,
  1077. C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
  1078. C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
  1079. C
  1080. C     CALLS CDIV FOR COMPLEX DIVISION.
  1081. C     CALLS CSROOT FOR COMPLEX SQUARE ROOT.
  1082. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  1083. C
  1084. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1085. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1086. C
  1087. C     THIS VERSION DATED AUGUST 1983.
  1088. C
  1089. C     ------------------------------------------------------------------
  1090. C
  1091.       SUBROUTINE COMQR2(NM,N,LOW,IGH,ORTR,ORTI,HR,HI,WR,WI,ZR,ZI,IERR)
  1092. C
  1093.       INTEGER I,J,K,L,M,N,EN,II,JJ,LL,NM,NN,IGH,IP1,
  1094.      X        ITN,ITS,LOW,LP1,ENM1,IEND,IERR
  1095.       REAL HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N),
  1096.      X       ORTR(IGH),ORTI(IGH)
  1097.       REAL SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2,
  1098.      X       PYTHAG
  1099. C
  1100. C     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
  1101. C     ALGOL PROCEDURE  COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS
  1102. C     AND WILKINSON.
  1103. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
  1104. C     THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
  1105. C     (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.
  1106. C
  1107. C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
  1108. C     OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR
  1109. C     METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX
  1110. C     CAN ALSO BE FOUND IF  CORTH  HAS BEEN USED TO REDUCE
  1111. C     THIS GENERAL MATRIX TO HESSENBERG FORM.
  1112. C
  1113. C     ON INPUT
  1114. C
  1115. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1116. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1117. C          DIMENSION STATEMENT.
  1118. C
  1119. C        N IS THE ORDER OF THE MATRIX.
  1120. C
  1121. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  1122. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
  1123. C          SET LOW=1, IGH=N.
  1124. C
  1125. C        ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
  1126. C          FORMATIONS USED IN THE REDUCTION BY  CORTH, IF PERFORMED.
  1127. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.  IF THE EIGENVECTORS
  1128. C          OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND
  1129. C          ORTI(J) TO 0.0E0 FOR THESE ELEMENTS.
  1130. C
  1131. C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
  1132. C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
  1133. C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER
  1134. C          INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE
  1135. C          REDUCTION BY  CORTH, IF PERFORMED.  IF THE EIGENVECTORS OF
  1136. C          THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE
  1137. C          ARBITRARY.
  1138. C
  1139. C     ON OUTPUT
  1140. C
  1141. C        ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI
  1142. C          HAVE BEEN DESTROYED.
  1143. C
  1144. C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
  1145. C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
  1146. C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
  1147. C          FOR INDICES IERR+1,...,N.
  1148. C
  1149. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  1150. C          RESPECTIVELY, OF THE EIGENVECTORS.  THE EIGENVECTORS
  1151. C          ARE UNNORMALIZED.  IF AN ERROR EXIT IS MADE, NONE OF
  1152. C          THE EIGENVECTORS HAS BEEN FOUND.
  1153. C
  1154. C        IERR IS SET TO
  1155. C          ZERO       FOR NORMAL RETURN,
  1156. C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
  1157. C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
  1158. C
  1159. C     CALLS CDIV FOR COMPLEX DIVISION.
  1160. C     CALLS CSROOT FOR COMPLEX SQUARE ROOT.
  1161. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  1162. C
  1163. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1164. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1165. C
  1166. C     THIS VERSION DATED AUGUST 1983.
  1167. C
  1168. C     ------------------------------------------------------------------
  1169. C
  1170.       SUBROUTINE CORTB(NM,LOW,IGH,AR,AI,ORTR,ORTI,M,ZR,ZI)
  1171. C
  1172.       INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1
  1173.       REAL AR(NM,IGH),AI(NM,IGH),ORTR(IGH),ORTI(IGH),
  1174.      X       ZR(NM,M),ZI(NM,M)
  1175.       REAL H,GI,GR
  1176. C
  1177. C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
  1178. C     THE ALGOL PROCEDURE ORTBAK, NUM. MATH. 12, 349-368(1968)
  1179. C     BY MARTIN AND WILKINSON.
  1180. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
  1181. C
  1182. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL
  1183. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  1184. C     UPPER HESSENBERG MATRIX DETERMINED BY  CORTH.
  1185. C
  1186. C     ON INPUT
  1187. C
  1188. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1189. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1190. C          DIMENSION STATEMENT.
  1191. C
  1192. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  1193. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
  1194. C          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX.
  1195. C
  1196. C        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY
  1197. C          TRANSFORMATIONS USED IN THE REDUCTION BY  CORTH
  1198. C          IN THEIR STRICT LOWER TRIANGLES.
  1199. C
  1200. C        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE
  1201. C          TRANSFORMATIONS USED IN THE REDUCTION BY  CORTH.
  1202. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  1203. C
  1204. C        M IS THE NUMBER OF COLUMNS OF ZR AND ZI TO BE BACK TRANSFORMED.
  1205. C
  1206. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  1207. C          RESPECTIVELY, OF THE EIGENVECTORS TO BE
  1208. C          BACK TRANSFORMED IN THEIR FIRST M COLUMNS.
  1209. C
  1210. C     ON OUTPUT
  1211. C
  1212. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  1213. C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
  1214. C          IN THEIR FIRST M COLUMNS.
  1215. C
  1216. C        ORTR AND ORTI HAVE BEEN ALTERED.
  1217. C
  1218. C     NOTE THAT CORTB PRESERVES VECTOR EUCLIDEAN NORMS.
  1219. C
  1220. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1221. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1222. C
  1223. C     THIS VERSION DATED AUGUST 1983.
  1224. C
  1225. C     ------------------------------------------------------------------
  1226. C
  1227.       SUBROUTINE CORTH(NM,N,LOW,IGH,AR,AI,ORTR,ORTI)
  1228. C
  1229.       INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW
  1230.       REAL AR(NM,N),AI(NM,N),ORTR(IGH),ORTI(IGH)
  1231.       REAL F,G,H,FI,FR,SCALE,PYTHAG
  1232. C
  1233. C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
  1234. C     THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968)
  1235. C     BY MARTIN AND WILKINSON.
  1236. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
  1237. C
  1238. C     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE
  1239. C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
  1240. C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
  1241. C     UNITARY SIMILARITY TRANSFORMATIONS.
  1242. C
  1243. C     ON INPUT
  1244. C
  1245. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1246. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1247. C          DIMENSION STATEMENT.
  1248. C
  1249. C        N IS THE ORDER OF THE MATRIX.
  1250. C
  1251. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  1252. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
  1253. C          SET LOW=1, IGH=N.
  1254. C
  1255. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
  1256. C          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX.
  1257. C
  1258. C     ON OUTPUT
  1259. C
  1260. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
  1261. C          RESPECTIVELY, OF THE HESSENBERG MATRIX.  INFORMATION
  1262. C          ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION
  1263. C          IS STORED IN THE REMAINING TRIANGLES UNDER THE
  1264. C          HESSENBERG MATRIX.
  1265. C
  1266. C        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE
  1267. C          TRANSFORMATIONS.  ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  1268. C
  1269. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  1270. C
  1271. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1272. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1273. C
  1274. C     THIS VERSION DATED AUGUST 1983.
  1275. C
  1276. C     ------------------------------------------------------------------
  1277. C
  1278.       SUBROUTINE ELMBAK(NM,LOW,IGH,A,INT,M,Z)
  1279. C
  1280.       INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1
  1281.       REAL A(NM,IGH),Z(NM,M)
  1282.       REAL X
  1283.       INTEGER INT(IGH)
  1284. C
  1285. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMBAK,
  1286. C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
  1287. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
  1288. C
  1289. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL
  1290. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  1291. C     UPPER HESSENBERG MATRIX DETERMINED BY  ELMHES.
  1292. C
  1293. C     ON INPUT
  1294. C
  1295. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1296. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1297. C          DIMENSION STATEMENT.
  1298. C
  1299. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  1300. C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
  1301. C          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX.
  1302. C
  1303. C        A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE
  1304. C          REDUCTION BY  ELMHES  IN ITS LOWER TRIANGLE
  1305. C          BELOW THE SUBDIAGONAL.
  1306. C
  1307. C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
  1308. C          INTERCHANGED IN THE REDUCTION BY  ELMHES.
  1309. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  1310. C
  1311. C        M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED.
  1312. C
  1313. C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN-
  1314. C          VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS.
  1315. C
  1316. C     ON OUTPUT
  1317. C
  1318. C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE
  1319. C          TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS.
  1320. C
  1321. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1322. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1323. C
  1324. C     THIS VERSION DATED AUGUST 1983.
  1325. C
  1326. C     ------------------------------------------------------------------
  1327. C
  1328.       SUBROUTINE ELMHES(NM,N,LOW,IGH,A,INT)
  1329. C
  1330.       INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1
  1331.       REAL A(NM,N)
  1332.       REAL X,Y
  1333.       INTEGER INT(IGH)
  1334. C
  1335. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES,
  1336. C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
  1337. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
  1338. C
  1339. C     GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE
  1340. C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
  1341. C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
  1342. C     STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS.
  1343. C
  1344. C     ON INPUT
  1345. C
  1346. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1347. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1348. C          DIMENSION STATEMENT.
  1349. C
  1350. C        N IS THE ORDER OF THE MATRIX.
  1351. C
  1352. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  1353. C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
  1354. C          SET LOW=1, IGH=N.
  1355. C
  1356. C        A CONTAINS THE INPUT MATRIX.
  1357. C
  1358. C     ON OUTPUT
  1359. C
  1360. C        A CONTAINS THE HESSENBERG MATRIX.  THE MULTIPLIERS
  1361. C          WHICH WERE USED IN THE REDUCTION ARE STORED IN THE
  1362. C          REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX.
  1363. C
  1364. C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
  1365. C          INTERCHANGED IN THE REDUCTION.
  1366. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  1367. C
  1368. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1369. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1370. C
  1371. C     THIS VERSION DATED AUGUST 1983.
  1372. C
  1373. C     ------------------------------------------------------------------
  1374. C
  1375.       SUBROUTINE ELTRAN(NM,N,LOW,IGH,A,INT,Z)
  1376. C
  1377.       INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1
  1378.       REAL A(NM,IGH),Z(NM,N)
  1379.       INTEGER INT(IGH)
  1380. C
  1381. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMTRANS,
  1382. C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
  1383. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
  1384. C
  1385. C     THIS SUBROUTINE ACCUMULATES THE STABILIZED ELEMENTARY
  1386. C     SIMILARITY TRANSFORMATIONS USED IN THE REDUCTION OF A
  1387. C     REAL GENERAL MATRIX TO UPPER HESSENBERG FORM BY  ELMHES.
  1388. C
  1389. C     ON INPUT
  1390. C
  1391. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1392. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1393. C          DIMENSION STATEMENT.
  1394. C
  1395. C        N IS THE ORDER OF THE MATRIX.
  1396. C
  1397. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  1398. C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
  1399. C          SET LOW=1, IGH=N.
  1400. C
  1401. C        A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE
  1402. C          REDUCTION BY  ELMHES  IN ITS LOWER TRIANGLE
  1403. C          BELOW THE SUBDIAGONAL.
  1404. C
  1405. C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
  1406. C          INTERCHANGED IN THE REDUCTION BY  ELMHES.
  1407. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  1408. C
  1409. C     ON OUTPUT
  1410. C
  1411. C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
  1412. C          REDUCTION BY  ELMHES.
  1413. C
  1414. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1415. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1416. C
  1417. C     THIS VERSION DATED AUGUST 1983.
  1418. C
  1419. C     ------------------------------------------------------------------
  1420. C
  1421. C
  1422.       SUBROUTINE FIGI(NM,N,T,D,E,E2,IERR)
  1423. C
  1424.       INTEGER I,N,NM,IERR
  1425.       REAL T(NM,3),D(N),E(N),E2(N)
  1426. C
  1427. C     GIVEN A NONSYMMETRIC TRIDIAGONAL MATRIX SUCH THAT THE PRODUCTS
  1428. C     OF CORRESPONDING PAIRS OF OFF-DIAGONAL ELEMENTS ARE ALL
  1429. C     NON-NEGATIVE, THIS SUBROUTINE REDUCES IT TO A SYMMETRIC
  1430. C     TRIDIAGONAL MATRIX WITH THE SAME EIGENVALUES.  IF, FURTHER,
  1431. C     A ZERO PRODUCT ONLY OCCURS WHEN BOTH FACTORS ARE ZERO,
  1432. C     THE REDUCED MATRIX IS SIMILAR TO THE ORIGINAL MATRIX.
  1433. C
  1434. C     ON INPUT
  1435. C
  1436. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1437. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1438. C          DIMENSION STATEMENT.
  1439. C
  1440. C        N IS THE ORDER OF THE MATRIX.
  1441. C
  1442. C        T CONTAINS THE INPUT MATRIX.  ITS SUBDIAGONAL IS
  1443. C          STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN,
  1444. C          ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN,
  1445. C          AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF
  1446. C          THE THIRD COLUMN.  T(1,1) AND T(N,3) ARE ARBITRARY.
  1447. C
  1448. C     ON OUTPUT
  1449. C
  1450. C        T IS UNALTERED.
  1451. C
  1452. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE SYMMETRIC MATRIX.
  1453. C
  1454. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC
  1455. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS NOT SET.
  1456. C
  1457. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  1458. C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
  1459. C
  1460. C        IERR IS SET TO
  1461. C          ZERO       FOR NORMAL RETURN,
  1462. C          N+I        IF T(I,1)*T(I-1,3) IS NEGATIVE,
  1463. C          -(3*N+I)   IF T(I,1)*T(I-1,3) IS ZERO WITH ONE FACTOR
  1464. C                     NON-ZERO.  IN THIS CASE, THE EIGENVECTORS OF
  1465. C                     THE SYMMETRIC MATRIX ARE NOT SIMPLY RELATED
  1466. C                     TO THOSE OF  T  AND SHOULD NOT BE SOUGHT.
  1467. C
  1468. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1469. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1470. C
  1471. C     THIS VERSION DATED AUGUST 1983.
  1472. C
  1473. C     ------------------------------------------------------------------
  1474. C
  1475.       SUBROUTINE FIGI2(NM,N,T,D,E,Z,IERR)
  1476. C
  1477.       INTEGER I,J,N,NM,IERR
  1478.       REAL T(NM,3),D(N),E(N),Z(NM,N)
  1479.       REAL H
  1480. C
  1481. C     GIVEN A NONSYMMETRIC TRIDIAGONAL MATRIX SUCH THAT THE PRODUCTS
  1482. C     OF CORRESPONDING PAIRS OF OFF-DIAGONAL ELEMENTS ARE ALL
  1483. C     NON-NEGATIVE, AND ZERO ONLY WHEN BOTH FACTORS ARE ZERO, THIS
  1484. C     SUBROUTINE REDUCES IT TO A SYMMETRIC TRIDIAGONAL MATRIX
  1485. C     USING AND ACCUMULATING DIAGONAL SIMILARITY TRANSFORMATIONS.
  1486. C
  1487. C     ON INPUT
  1488. C
  1489. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1490. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1491. C          DIMENSION STATEMENT.
  1492. C
  1493. C        N IS THE ORDER OF THE MATRIX.
  1494. C
  1495. C        T CONTAINS THE INPUT MATRIX.  ITS SUBDIAGONAL IS
  1496. C          STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN,
  1497. C          ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN,
  1498. C          AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF
  1499. C          THE THIRD COLUMN.  T(1,1) AND T(N,3) ARE ARBITRARY.
  1500. C
  1501. C     ON OUTPUT
  1502. C
  1503. C        T IS UNALTERED.
  1504. C
  1505. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE SYMMETRIC MATRIX.
  1506. C
  1507. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC
  1508. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS NOT SET.
  1509. C
  1510. C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN
  1511. C          THE REDUCTION.
  1512. C
  1513. C        IERR IS SET TO
  1514. C          ZERO       FOR NORMAL RETURN,
  1515. C          N+I        IF T(I,1)*T(I-1,3) IS NEGATIVE,
  1516. C          2*N+I      IF T(I,1)*T(I-1,3) IS ZERO WITH
  1517. C                     ONE FACTOR NON-ZERO.
  1518. C
  1519. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1520. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1521. C
  1522. C     THIS VERSION DATED AUGUST 1983.
  1523. C
  1524. C     ------------------------------------------------------------------
  1525. C
  1526.       SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR)
  1527. C
  1528.       INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR
  1529.       REAL H(NM,N),WR(N),WI(N)
  1530.       REAL P,Q,R,S,T,W,X,Y,ZZ,NORM,TST1,TST2
  1531.       LOGICAL NOTLAS
  1532. C
  1533. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR,
  1534. C     NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON.
  1535. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971).
  1536. C
  1537. C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL
  1538. C     UPPER HESSENBERG MATRIX BY THE QR METHOD.
  1539. C
  1540. C     ON INPUT
  1541. C
  1542. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1543. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1544. C          DIMENSION STATEMENT.
  1545. C
  1546. C        N IS THE ORDER OF THE MATRIX.
  1547. C
  1548. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  1549. C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
  1550. C          SET LOW=1, IGH=N.
  1551. C
  1552. C        H CONTAINS THE UPPER HESSENBERG MATRIX.  INFORMATION ABOUT
  1553. C          THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG
  1554. C          FORM BY  ELMHES  OR  ORTHES, IF PERFORMED, IS STORED
  1555. C          IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX.
  1556. C
  1557. C     ON OUTPUT
  1558. C
  1559. C        H HAS BEEN DESTROYED.  THEREFORE, IT MUST BE SAVED
  1560. C          BEFORE CALLING  HQR  IF SUBSEQUENT CALCULATION AND
  1561. C          BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED.
  1562. C
  1563. C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
  1564. C          RESPECTIVELY, OF THE EIGENVALUES.  THE EIGENVALUES
  1565. C          ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS
  1566. C          OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE
  1567. C          HAVING THE POSITIVE IMAGINARY PART FIRST.  IF AN
  1568. C          ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
  1569. C          FOR INDICES IERR+1,...,N.
  1570. C
  1571. C        IERR IS SET TO
  1572. C          ZERO       FOR NORMAL RETURN,
  1573. C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
  1574. C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
  1575. C
  1576. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1577. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1578. C
  1579. C     THIS VERSION DATED AUGUST 1983.
  1580. C
  1581. C     ------------------------------------------------------------------
  1582. C
  1583.       SUBROUTINE HQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR)
  1584. C
  1585.       INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN,
  1586.      X        IGH,ITN,ITS,LOW,MP2,ENM2,IERR
  1587.       REAL H(NM,N),WR(N),WI(N),Z(NM,N)
  1588.       REAL P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,TST1,TST2
  1589.       LOGICAL NOTLAS
  1590. C
  1591. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2,
  1592. C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
  1593. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
  1594. C
  1595. C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
  1596. C     OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD.  THE
  1597. C     EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND
  1598. C     IF  ELMHES  AND  ELTRAN  OR  ORTHES  AND  ORTRAN  HAVE
  1599. C     BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM
  1600. C     AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS.
  1601. C
  1602. C     ON INPUT
  1603. C
  1604. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1605. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1606. C          DIMENSION STATEMENT.
  1607. C
  1608. C        N IS THE ORDER OF THE MATRIX.
  1609. C
  1610. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  1611. C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
  1612. C          SET LOW=1, IGH=N.
  1613. C
  1614. C        H CONTAINS THE UPPER HESSENBERG MATRIX.
  1615. C
  1616. C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY  ELTRAN
  1617. C          AFTER THE REDUCTION BY  ELMHES, OR BY  ORTRAN  AFTER THE
  1618. C          REDUCTION BY  ORTHES, IF PERFORMED.  IF THE EIGENVECTORS
  1619. C          OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE
  1620. C          IDENTITY MATRIX.
  1621. C
  1622. C     ON OUTPUT
  1623. C
  1624. C        H HAS BEEN DESTROYED.
  1625. C
  1626. C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
  1627. C          RESPECTIVELY, OF THE EIGENVALUES.  THE EIGENVALUES
  1628. C          ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS
  1629. C          OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE
  1630. C          HAVING THE POSITIVE IMAGINARY PART FIRST.  IF AN
  1631. C          ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
  1632. C          FOR INDICES IERR+1,...,N.
  1633. C
  1634. C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS.
  1635. C          IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z
  1636. C          CONTAINS ITS EIGENVECTOR.  IF THE I-TH EIGENVALUE IS COMPLEX
  1637. C          WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH
  1638. C          COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS
  1639. C          EIGENVECTOR.  THE EIGENVECTORS ARE UNNORMALIZED.  IF AN
  1640. C          ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND.
  1641. C
  1642. C        IERR IS SET TO
  1643. C          ZERO       FOR NORMAL RETURN,
  1644. C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
  1645. C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
  1646. C
  1647. C     CALLS CDIV FOR COMPLEX DIVISION.
  1648. C
  1649. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1650. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1651. C
  1652. C     THIS VERSION DATED AUGUST 1983.
  1653. C
  1654. C     ------------------------------------------------------------------
  1655. C
  1656.       SUBROUTINE HTRIB3(NM,N,A,TAU,M,ZR,ZI)
  1657. C
  1658.       INTEGER I,J,K,L,M,N,NM
  1659.       REAL A(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M)
  1660.       REAL H,S,SI
  1661. C
  1662. C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
  1663. C     THE ALGOL PROCEDURE TRBAK3, NUM. MATH. 11, 181-195(1968)
  1664. C     BY MARTIN, REINSCH, AND WILKINSON.
  1665. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  1666. C
  1667. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN
  1668. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  1669. C     REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  HTRID3.
  1670. C
  1671. C     ON INPUT
  1672. C
  1673. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1674. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1675. C          DIMENSION STATEMENT.
  1676. C
  1677. C        N IS THE ORDER OF THE MATRIX.
  1678. C
  1679. C        A CONTAINS INFORMATION ABOUT THE UNITARY TRANSFORMATIONS
  1680. C          USED IN THE REDUCTION BY  HTRID3.
  1681. C
  1682. C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
  1683. C
  1684. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
  1685. C
  1686. C        ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
  1687. C          IN ITS FIRST M COLUMNS.
  1688. C
  1689. C     ON OUTPUT
  1690. C
  1691. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  1692. C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
  1693. C          IN THEIR FIRST M COLUMNS.
  1694. C
  1695. C     NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR
  1696. C     IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED.
  1697. C
  1698. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1699. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1700. C
  1701. C     THIS VERSION DATED AUGUST 1983.
  1702. C
  1703. C     ------------------------------------------------------------------
  1704. C
  1705.       SUBROUTINE HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)
  1706. C
  1707.       INTEGER I,J,K,L,M,N,NM
  1708.       REAL AR(NM,N),AI(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M)
  1709.       REAL H,S,SI
  1710. C
  1711. C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
  1712. C     THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968)
  1713. C     BY MARTIN, REINSCH, AND WILKINSON.
  1714. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  1715. C
  1716. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN
  1717. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  1718. C     REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  HTRIDI.
  1719. C
  1720. C     ON INPUT
  1721. C
  1722. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1723. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1724. C          DIMENSION STATEMENT.
  1725. C
  1726. C        N IS THE ORDER OF THE MATRIX.
  1727. C
  1728. C        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
  1729. C          FORMATIONS USED IN THE REDUCTION BY  HTRIDI  IN THEIR
  1730. C          FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR.
  1731. C
  1732. C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
  1733. C
  1734. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
  1735. C
  1736. C        ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
  1737. C          IN ITS FIRST M COLUMNS.
  1738. C
  1739. C     ON OUTPUT
  1740. C
  1741. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
  1742. C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
  1743. C          IN THEIR FIRST M COLUMNS.
  1744. C
  1745. C     NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR
  1746. C     IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED.
  1747. C
  1748. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1749. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1750. C
  1751. C     THIS VERSION DATED AUGUST 1983.
  1752. C
  1753. C     ------------------------------------------------------------------
  1754. C
  1755.       SUBROUTINE HTRID3(NM,N,A,D,E,E2,TAU)
  1756. C
  1757.       INTEGER I,J,K,L,N,II,NM,JM1,JP1
  1758.       REAL A(NM,N),D(N),E(N),E2(N),TAU(2,N)
  1759.       REAL F,G,H,FI,GI,HH,SI,SCALE,PYTHAG
  1760. C
  1761. C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
  1762. C     THE ALGOL PROCEDURE TRED3, NUM. MATH. 11, 181-195(1968)
  1763. C     BY MARTIN, REINSCH, AND WILKINSON.
  1764. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  1765. C
  1766. C     THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX, STORED AS
  1767. C     A SINGLE SQUARE ARRAY, TO A REAL SYMMETRIC TRIDIAGONAL MATRIX
  1768. C     USING UNITARY SIMILARITY TRANSFORMATIONS.
  1769. C
  1770. C     ON INPUT
  1771. C
  1772. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1773. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1774. C          DIMENSION STATEMENT.
  1775. C
  1776. C        N IS THE ORDER OF THE MATRIX.
  1777. C
  1778. C        A CONTAINS THE LOWER TRIANGLE OF THE COMPLEX HERMITIAN INPUT
  1779. C          MATRIX.  THE REAL PARTS OF THE MATRIX ELEMENTS ARE STORED
  1780. C          IN THE FULL LOWER TRIANGLE OF A, AND THE IMAGINARY PARTS
  1781. C          ARE STORED IN THE TRANSPOSED POSITIONS OF THE STRICT UPPER
  1782. C          TRIANGLE OF A.  NO STORAGE IS REQUIRED FOR THE ZERO
  1783. C          IMAGINARY PARTS OF THE DIAGONAL ELEMENTS.
  1784. C
  1785. C     ON OUTPUT
  1786. C
  1787. C        A CONTAINS INFORMATION ABOUT THE UNITARY TRANSFORMATIONS
  1788. C          USED IN THE REDUCTION.
  1789. C
  1790. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX.
  1791. C
  1792. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  1793. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
  1794. C
  1795. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  1796. C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
  1797. C
  1798. C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
  1799. C
  1800. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  1801. C
  1802. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1803. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1804. C
  1805. C     THIS VERSION DATED AUGUST 1983.
  1806. C
  1807. C     ------------------------------------------------------------------
  1808. C
  1809.       SUBROUTINE HTRIDI(NM,N,AR,AI,D,E,E2,TAU)
  1810. C
  1811.       INTEGER I,J,K,L,N,II,NM,JP1
  1812.       REAL AR(NM,N),AI(NM,N),D(N),E(N),E2(N),TAU(2,N)
  1813.       REAL F,G,H,FI,GI,HH,SI,SCALE,PYTHAG
  1814. C
  1815. C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
  1816. C     THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968)
  1817. C     BY MARTIN, REINSCH, AND WILKINSON.
  1818. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  1819. C
  1820. C     THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX
  1821. C     TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING
  1822. C     UNITARY SIMILARITY TRANSFORMATIONS.
  1823. C
  1824. C     ON INPUT
  1825. C
  1826. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1827. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1828. C          DIMENSION STATEMENT.
  1829. C
  1830. C        N IS THE ORDER OF THE MATRIX.
  1831. C
  1832. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
  1833. C          RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX.
  1834. C          ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
  1835. C
  1836. C     ON OUTPUT
  1837. C
  1838. C        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
  1839. C          FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER
  1840. C          TRIANGLES.  THEIR STRICT UPPER TRIANGLES AND THE
  1841. C          DIAGONAL OF AR ARE UNALTERED.
  1842. C
  1843. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX.
  1844. C
  1845. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  1846. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
  1847. C
  1848. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  1849. C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
  1850. C
  1851. C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
  1852. C
  1853. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  1854. C
  1855. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1856. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1857. C
  1858. C     THIS VERSION DATED AUGUST 1983.
  1859. C
  1860. C     ------------------------------------------------------------------
  1861. C
  1862.       SUBROUTINE IMTQL1(N,D,E,IERR)
  1863. C
  1864.       INTEGER I,J,L,M,N,II,MML,IERR
  1865.       REAL D(N),E(N)
  1866.       REAL B,C,F,G,P,R,S,TST1,TST2,PYTHAG
  1867. C
  1868. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL1,
  1869. C     NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON,
  1870. C     AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
  1871. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
  1872. C
  1873. C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
  1874. C     TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD.
  1875. C
  1876. C     ON INPUT
  1877. C
  1878. C        N IS THE ORDER OF THE MATRIX.
  1879. C
  1880. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  1881. C
  1882. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  1883. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  1884. C
  1885. C      ON OUTPUT
  1886. C
  1887. C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
  1888. C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
  1889. C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
  1890. C          THE SMALLEST EIGENVALUES.
  1891. C
  1892. C        E HAS BEEN DESTROYED.
  1893. C
  1894. C        IERR IS SET TO
  1895. C          ZERO       FOR NORMAL RETURN,
  1896. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
  1897. C                     DETERMINED AFTER 30 ITERATIONS.
  1898. C
  1899. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  1900. C
  1901. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1902. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1903. C
  1904. C     THIS VERSION DATED AUGUST 1983.
  1905. C
  1906. C     ------------------------------------------------------------------
  1907. C
  1908.       SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR)
  1909. C
  1910.       INTEGER I,J,K,L,M,N,II,NM,MML,IERR
  1911.       REAL D(N),E(N),Z(NM,N)
  1912.       REAL B,C,F,G,P,R,S,TST1,TST2,PYTHAG
  1913. C
  1914. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2,
  1915. C     NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON,
  1916. C     AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
  1917. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
  1918. C
  1919. C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
  1920. C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD.
  1921. C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
  1922. C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
  1923. C     FULL MATRIX TO TRIDIAGONAL FORM.
  1924. C
  1925. C     ON INPUT
  1926. C
  1927. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  1928. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  1929. C          DIMENSION STATEMENT.
  1930. C
  1931. C        N IS THE ORDER OF THE MATRIX.
  1932. C
  1933. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  1934. C
  1935. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  1936. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  1937. C
  1938. C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
  1939. C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
  1940. C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
  1941. C          THE IDENTITY MATRIX.
  1942. C
  1943. C      ON OUTPUT
  1944. C
  1945. C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
  1946. C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
  1947. C          UNORDERED FOR INDICES 1,2,...,IERR-1.
  1948. C
  1949. C        E HAS BEEN DESTROYED.
  1950. C
  1951. C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
  1952. C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
  1953. C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
  1954. C          EIGENVALUES.
  1955. C
  1956. C        IERR IS SET TO
  1957. C          ZERO       FOR NORMAL RETURN,
  1958. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
  1959. C                     DETERMINED AFTER 30 ITERATIONS.
  1960. C
  1961. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  1962. C
  1963. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  1964. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  1965. C
  1966. C     THIS VERSION DATED AUGUST 1983.
  1967. C
  1968. C     ------------------------------------------------------------------
  1969. C
  1970.       SUBROUTINE IMTQLV(N,D,E,E2,W,IND,IERR,RV1)
  1971. C
  1972.       INTEGER I,J,K,L,M,N,II,MML,TAG,IERR
  1973.       REAL D(N),E(N),E2(N),W(N),RV1(N)
  1974.       REAL B,C,F,G,P,R,S,TST1,TST2,PYTHAG
  1975.       INTEGER IND(N)
  1976. C
  1977. C     THIS SUBROUTINE IS A VARIANT OF  IMTQL1  WHICH IS A TRANSLATION OF
  1978. C     ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND
  1979. C     WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
  1980. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
  1981. C
  1982. C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
  1983. C     MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM
  1984. C     THEIR CORRESPONDING SUBMATRIX INDICES.
  1985. C
  1986. C     ON INPUT
  1987. C
  1988. C        N IS THE ORDER OF THE MATRIX.
  1989. C
  1990. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  1991. C
  1992. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  1993. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  1994. C
  1995. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  1996. C          E2(1) IS ARBITRARY.
  1997. C
  1998. C     ON OUTPUT
  1999. C
  2000. C        D AND E ARE UNALTERED.
  2001. C
  2002. C        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
  2003. C          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
  2004. C          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
  2005. C          E2(1) IS ALSO SET TO ZERO.
  2006. C
  2007. C        W CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
  2008. C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
  2009. C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
  2010. C          THE SMALLEST EIGENVALUES.
  2011. C
  2012. C        IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE
  2013. C          CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES
  2014. C          BELONGING TO THE FIRST SUBMATRIX FROM THE TOP,
  2015. C          2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC..
  2016. C
  2017. C        IERR IS SET TO
  2018. C          ZERO       FOR NORMAL RETURN,
  2019. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
  2020. C                     DETERMINED AFTER 30 ITERATIONS.
  2021. C
  2022. C        RV1 IS A TEMPORARY STORAGE ARRAY.
  2023. C
  2024. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  2025. C
  2026. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2027. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2028. C
  2029. C     THIS VERSION DATED AUGUST 1983.
  2030. C
  2031. C     ------------------------------------------------------------------
  2032. C
  2033.       SUBROUTINE INVIT(NM,N,A,WR,WI,SELECT,MM,M,Z,IERR,RM1,RV1,RV2)
  2034. C
  2035.       INTEGER I,J,K,L,M,N,S,II,IP,MM,MP,NM,NS,N1,UK,IP1,ITS,KM1,IERR
  2036.       REAL A(NM,N),WR(N),WI(N),Z(NM,MM),RM1(N,N),
  2037.      X       RV1(N),RV2(N)
  2038.       REAL T,W,X,Y,EPS3,NORM,NORMV,EPSLON,GROWTO,ILAMBD,
  2039.      X       PYTHAG,RLAMBD,UKROOT
  2040.       LOGICAL SELECT(N)
  2041. C
  2042. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE INVIT
  2043. C     BY PETERS AND WILKINSON.
  2044. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
  2045. C
  2046. C     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL UPPER
  2047. C     HESSENBERG MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
  2048. C     USING INVERSE ITERATION.
  2049. C
  2050. C     ON INPUT
  2051. C
  2052. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2053. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2054. C          DIMENSION STATEMENT.
  2055. C
  2056. C        N IS THE ORDER OF THE MATRIX.
  2057. C
  2058. C        A CONTAINS THE HESSENBERG MATRIX.
  2059. C
  2060. C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY,
  2061. C          OF THE EIGENVALUES OF THE MATRIX.  THE EIGENVALUES MUST BE
  2062. C          STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE  HQR,
  2063. C          WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX.
  2064. C
  2065. C        SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND. THE
  2066. C          EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS
  2067. C          SPECIFIED BY SETTING SELECT(J) TO .TRUE..
  2068. C
  2069. C        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF
  2070. C          COLUMNS REQUIRED TO STORE THE EIGENVECTORS TO BE FOUND.
  2071. C          NOTE THAT TWO COLUMNS ARE REQUIRED TO STORE THE
  2072. C          EIGENVECTOR CORRESPONDING TO A COMPLEX EIGENVALUE.
  2073. C
  2074. C     ON OUTPUT
  2075. C
  2076. C        A AND WI ARE UNALTERED.
  2077. C
  2078. C        WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED
  2079. C          SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS.
  2080. C
  2081. C        SELECT MAY HAVE BEEN ALTERED.  IF THE ELEMENTS CORRESPONDING
  2082. C          TO A PAIR OF CONJUGATE COMPLEX EIGENVALUES WERE EACH
  2083. C          INITIALLY SET TO .TRUE., THE PROGRAM RESETS THE SECOND OF
  2084. C          THE TWO ELEMENTS TO .FALSE..
  2085. C
  2086. C        M IS THE NUMBER OF COLUMNS ACTUALLY USED TO STORE
  2087. C          THE EIGENVECTORS.
  2088. C
  2089. C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS.
  2090. C          IF THE NEXT SELECTED EIGENVALUE IS REAL, THE NEXT COLUMN
  2091. C          OF Z CONTAINS ITS EIGENVECTOR.  IF THE EIGENVALUE IS
  2092. C          COMPLEX, THE NEXT TWO COLUMNS OF Z CONTAIN THE REAL AND
  2093. C          IMAGINARY PARTS OF ITS EIGENVECTOR.  THE EIGENVECTORS ARE
  2094. C          NORMALIZED SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1.
  2095. C          ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO.
  2096. C
  2097. C        IERR IS SET TO
  2098. C          ZERO       FOR NORMAL RETURN,
  2099. C          -(2*N+1)   IF MORE THAN MM COLUMNS OF Z ARE NECESSARY
  2100. C                     TO STORE THE EIGENVECTORS CORRESPONDING TO
  2101. C                     THE SPECIFIED EIGENVALUES.
  2102. C          -K         IF THE ITERATION CORRESPONDING TO THE K-TH
  2103. C                     VALUE FAILS,
  2104. C          -(N+K)     IF BOTH ERROR SITUATIONS OCCUR.
  2105. C
  2106. C        RM1, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS.  NOTE THAT RM1
  2107. C          IS SQUARE OF DIMENSION N BY N AND, AUGMENTED BY TWO COLUMNS
  2108. C          OF Z, IS THE TRANSPOSE OF THE CORRESPONDING ALGOL B ARRAY.
  2109. C
  2110. C     THE ALGOL PROCEDURE GUESSVEC APPEARS IN INVIT IN LINE.
  2111. C
  2112. C     CALLS CDIV FOR COMPLEX DIVISION.
  2113. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  2114. C
  2115. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2116. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2117. C
  2118. C     THIS VERSION DATED AUGUST 1983.
  2119. C
  2120. C     ------------------------------------------------------------------
  2121. C
  2122.       SUBROUTINE MINFIT(NM,M,N,A,W,IP,B,IERR,RV1)
  2123. C
  2124.       INTEGER I,J,K,L,M,N,II,IP,I1,KK,K1,LL,L1,M1,NM,ITS,IERR
  2125.       REAL A(NM,N),W(N),B(NM,IP),RV1(N)
  2126.       REAL C,F,G,H,S,X,Y,Z,TST1,TST2,SCALE,PYTHAG
  2127. C
  2128. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE MINFIT,
  2129. C     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH.
  2130. C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
  2131. C
  2132. C     THIS SUBROUTINE DETERMINES, TOWARDS THE SOLUTION OF THE LINEAR
  2133. C                                                        T
  2134. C     SYSTEM AX=B, THE SINGULAR VALUE DECOMPOSITION A=USV  OF A REAL
  2135. C                                         T
  2136. C     M BY N RECTANGULAR MATRIX, FORMING U B RATHER THAN U.  HOUSEHOLDER
  2137. C     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED.
  2138. C
  2139. C     ON INPUT
  2140. C
  2141. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2142. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2143. C          DIMENSION STATEMENT.  NOTE THAT NM MUST BE AT LEAST
  2144. C          AS LARGE AS THE MAXIMUM OF M AND N.
  2145. C
  2146. C        M IS THE NUMBER OF ROWS OF A AND B.
  2147. C
  2148. C        N IS THE NUMBER OF COLUMNS OF A AND THE ORDER OF V.
  2149. C
  2150. C        A CONTAINS THE RECTANGULAR COEFFICIENT MATRIX OF THE SYSTEM.
  2151. C
  2152. C        IP IS THE NUMBER OF COLUMNS OF B.  IP CAN BE ZERO.
  2153. C
  2154. C        B CONTAINS THE CONSTANT COLUMN MATRIX OF THE SYSTEM
  2155. C          IF IP IS NOT ZERO.  OTHERWISE B IS NOT REFERENCED.
  2156. C
  2157. C     ON OUTPUT
  2158. C
  2159. C        A HAS BEEN OVERWRITTEN BY THE MATRIX V (ORTHOGONAL) OF THE
  2160. C          DECOMPOSITION IN ITS FIRST N ROWS AND COLUMNS.  IF AN
  2161. C          ERROR EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO
  2162. C          INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT.
  2163. C
  2164. C        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE
  2165. C          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN
  2166. C          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT
  2167. C          FOR INDICES IERR+1,IERR+2,...,N.
  2168. C
  2169. C                                   T
  2170. C        B HAS BEEN OVERWRITTEN BY U B.  IF AN ERROR EXIT IS MADE,
  2171. C                       T
  2172. C          THE ROWS OF U B CORRESPONDING TO INDICES OF CORRECT
  2173. C          SINGULAR VALUES SHOULD BE CORRECT.
  2174. C
  2175. C        IERR IS SET TO
  2176. C          ZERO       FOR NORMAL RETURN,
  2177. C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
  2178. C                     DETERMINED AFTER 30 ITERATIONS.
  2179. C
  2180. C        RV1 IS A TEMPORARY STORAGE ARRAY.
  2181. C
  2182. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  2183. C
  2184. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2185. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2186. C
  2187. C     THIS VERSION DATED AUGUST 1983.
  2188. C
  2189. C     ------------------------------------------------------------------
  2190. C
  2191.       SUBROUTINE ORTBAK(NM,LOW,IGH,A,ORT,M,Z)
  2192. C
  2193.       INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1
  2194.       REAL A(NM,IGH),ORT(IGH),Z(NM,M)
  2195.       REAL G
  2196. C
  2197. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTBAK,
  2198. C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
  2199. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
  2200. C
  2201. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL
  2202. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  2203. C     UPPER HESSENBERG MATRIX DETERMINED BY  ORTHES.
  2204. C
  2205. C     ON INPUT
  2206. C
  2207. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2208. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2209. C          DIMENSION STATEMENT.
  2210. C
  2211. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  2212. C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
  2213. C          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX.
  2214. C
  2215. C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
  2216. C          FORMATIONS USED IN THE REDUCTION BY  ORTHES
  2217. C          IN ITS STRICT LOWER TRIANGLE.
  2218. C
  2219. C        ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANS-
  2220. C          FORMATIONS USED IN THE REDUCTION BY  ORTHES.
  2221. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  2222. C
  2223. C        M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED.
  2224. C
  2225. C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN-
  2226. C          VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS.
  2227. C
  2228. C     ON OUTPUT
  2229. C
  2230. C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE
  2231. C          TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS.
  2232. C
  2233. C        ORT HAS BEEN ALTERED.
  2234. C
  2235. C     NOTE THAT ORTBAK PRESERVES VECTOR EUCLIDEAN NORMS.
  2236. C
  2237. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2238. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2239. C
  2240. C     THIS VERSION DATED AUGUST 1983.
  2241. C
  2242. C     ------------------------------------------------------------------
  2243. C
  2244.       SUBROUTINE ORTHES(NM,N,LOW,IGH,A,ORT)
  2245. C
  2246.       INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW
  2247.       REAL A(NM,N),ORT(IGH)
  2248.       REAL F,G,H,SCALE
  2249. C
  2250. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTHES,
  2251. C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
  2252. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
  2253. C
  2254. C     GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE
  2255. C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
  2256. C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
  2257. C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
  2258. C
  2259. C     ON INPUT
  2260. C
  2261. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2262. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2263. C          DIMENSION STATEMENT.
  2264. C
  2265. C        N IS THE ORDER OF THE MATRIX.
  2266. C
  2267. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  2268. C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
  2269. C          SET LOW=1, IGH=N.
  2270. C
  2271. C        A CONTAINS THE INPUT MATRIX.
  2272. C
  2273. C     ON OUTPUT
  2274. C
  2275. C        A CONTAINS THE HESSENBERG MATRIX.  INFORMATION ABOUT
  2276. C          THE ORTHOGONAL TRANSFORMATIONS USED IN THE REDUCTION
  2277. C          IS STORED IN THE REMAINING TRIANGLE UNDER THE
  2278. C          HESSENBERG MATRIX.
  2279. C
  2280. C        ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
  2281. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  2282. C
  2283. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2284. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2285. C
  2286. C     THIS VERSION DATED AUGUST 1983.
  2287. C
  2288. C     ------------------------------------------------------------------
  2289. C
  2290.       SUBROUTINE ORTRAN(NM,N,LOW,IGH,A,ORT,Z)
  2291. C
  2292.       INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1
  2293.       REAL A(NM,IGH),ORT(IGH),Z(NM,N)
  2294.       REAL G
  2295. C
  2296. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTRANS,
  2297. C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
  2298. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
  2299. C
  2300. C     THIS SUBROUTINE ACCUMULATES THE ORTHOGONAL SIMILARITY
  2301. C     TRANSFORMATIONS USED IN THE REDUCTION OF A REAL GENERAL
  2302. C     MATRIX TO UPPER HESSENBERG FORM BY  ORTHES.
  2303. C
  2304. C     ON INPUT
  2305. C
  2306. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2307. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2308. C          DIMENSION STATEMENT.
  2309. C
  2310. C        N IS THE ORDER OF THE MATRIX.
  2311. C
  2312. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
  2313. C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
  2314. C          SET LOW=1, IGH=N.
  2315. C
  2316. C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
  2317. C          FORMATIONS USED IN THE REDUCTION BY  ORTHES
  2318. C          IN ITS STRICT LOWER TRIANGLE.
  2319. C
  2320. C        ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANS-
  2321. C          FORMATIONS USED IN THE REDUCTION BY  ORTHES.
  2322. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  2323. C
  2324. C     ON OUTPUT
  2325. C
  2326. C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
  2327. C          REDUCTION BY  ORTHES.
  2328. C
  2329. C        ORT HAS BEEN ALTERED.
  2330. C
  2331. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2332. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2333. C
  2334. C     THIS VERSION DATED AUGUST 1983.
  2335. C
  2336. C     ------------------------------------------------------------------
  2337. C
  2338. C
  2339.       SUBROUTINE QZHES(NM,N,A,B,MATZ,Z)
  2340. C
  2341.       INTEGER I,J,K,L,N,LB,L1,NM,NK1,NM1,NM2
  2342.       REAL A(NM,N),B(NM,N),Z(NM,N)
  2343.       REAL R,S,T,U1,U2,V1,V2,RHO
  2344.       LOGICAL MATZ
  2345. C
  2346. C     THIS SUBROUTINE IS THE FIRST STEP OF THE QZ ALGORITHM
  2347. C     FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS,
  2348. C     SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART.
  2349. C
  2350. C     THIS SUBROUTINE ACCEPTS A PAIR OF REAL GENERAL MATRICES AND
  2351. C     REDUCES ONE OF THEM TO UPPER HESSENBERG FORM AND THE OTHER
  2352. C     TO UPPER TRIANGULAR FORM USING ORTHOGONAL TRANSFORMATIONS.
  2353. C     IT IS USUALLY FOLLOWED BY  QZIT,  QZVAL  AND, POSSIBLY,  QZVEC.
  2354. C
  2355. C     ON INPUT
  2356. C
  2357. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2358. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2359. C          DIMENSION STATEMENT.
  2360. C
  2361. C        N IS THE ORDER OF THE MATRICES.
  2362. C
  2363. C        A CONTAINS A REAL GENERAL MATRIX.
  2364. C
  2365. C        B CONTAINS A REAL GENERAL MATRIX.
  2366. C
  2367. C        MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS
  2368. C          ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING
  2369. C          EIGENVECTORS, AND TO .FALSE. OTHERWISE.
  2370. C
  2371. C     ON OUTPUT
  2372. C
  2373. C        A HAS BEEN REDUCED TO UPPER HESSENBERG FORM.  THE ELEMENTS
  2374. C          BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO.
  2375. C
  2376. C        B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM.  THE ELEMENTS
  2377. C          BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO.
  2378. C
  2379. C        Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS IF
  2380. C          MATZ HAS BEEN SET TO .TRUE.  OTHERWISE, Z IS NOT REFERENCED.
  2381. C
  2382. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2383. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2384. C
  2385. C     THIS VERSION DATED AUGUST 1983.
  2386. C
  2387. C     ------------------------------------------------------------------
  2388. C
  2389. C
  2390.       SUBROUTINE QZIT(NM,N,A,B,EPS1,MATZ,Z,IERR)
  2391. C
  2392.       INTEGER I,J,K,L,N,EN,K1,K2,LD,LL,L1,NA,NM,ISH,ITN,ITS,KM1,LM1,
  2393.      X        ENM2,IERR,LOR1,ENORN
  2394.       REAL A(NM,N),B(NM,N),Z(NM,N)
  2395.       REAL R,S,T,A1,A2,A3,EP,SH,U1,U2,U3,V1,V2,V3,ANI,A11,
  2396.      X       A12,A21,A22,A33,A34,A43,A44,BNI,B11,B12,B22,B33,B34,
  2397.      X       B44,EPSA,EPSB,EPS1,ANORM,BNORM,EPSLON
  2398.       LOGICAL MATZ,NOTLAS
  2399. C
  2400. C     THIS SUBROUTINE IS THE SECOND STEP OF THE QZ ALGORITHM
  2401. C     FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS,
  2402. C     SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART,
  2403. C     AS MODIFIED IN TECHNICAL NOTE NASA TN D-7305(1973) BY WARD.
  2404. C
  2405. C     THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM
  2406. C     IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER TRIANGULAR FORM.
  2407. C     IT REDUCES THE HESSENBERG MATRIX TO QUASI-TRIANGULAR FORM USING
  2408. C     ORTHOGONAL TRANSFORMATIONS WHILE MAINTAINING THE TRIANGULAR FORM
  2409. C     OF THE OTHER MATRIX.  IT IS USUALLY PRECEDED BY  QZHES  AND
  2410. C     FOLLOWED BY  QZVAL  AND, POSSIBLY,  QZVEC.
  2411. C
  2412. C     ON INPUT
  2413. C
  2414. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2415. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2416. C          DIMENSION STATEMENT.
  2417. C
  2418. C        N IS THE ORDER OF THE MATRICES.
  2419. C
  2420. C        A CONTAINS A REAL UPPER HESSENBERG MATRIX.
  2421. C
  2422. C        B CONTAINS A REAL UPPER TRIANGULAR MATRIX.
  2423. C
  2424. C        EPS1 IS A TOLERANCE USED TO DETERMINE NEGLIGIBLE ELEMENTS.
  2425. C          EPS1 = 0.0 (OR NEGATIVE) MAY BE INPUT, IN WHICH CASE AN
  2426. C          ELEMENT WILL BE NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF
  2427. C          ERROR TIMES THE NORM OF ITS MATRIX.  IF THE INPUT EPS1 IS
  2428. C          POSITIVE, THEN AN ELEMENT WILL BE CONSIDERED NEGLIGIBLE
  2429. C          IF IT IS LESS THAN EPS1 TIMES THE NORM OF ITS MATRIX.  A
  2430. C          POSITIVE VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION,
  2431. C          BUT LESS ACCURATE RESULTS.
  2432. C
  2433. C        MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS
  2434. C          ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING
  2435. C          EIGENVECTORS, AND TO .FALSE. OTHERWISE.
  2436. C
  2437. C        Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE
  2438. C          TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION
  2439. C          BY  QZHES, IF PERFORMED, OR ELSE THE IDENTITY MATRIX.
  2440. C          IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED.
  2441. C
  2442. C     ON OUTPUT
  2443. C
  2444. C        A HAS BEEN REDUCED TO QUASI-TRIANGULAR FORM.  THE ELEMENTS
  2445. C          BELOW THE FIRST SUBDIAGONAL ARE STILL ZERO AND NO TWO
  2446. C          CONSECUTIVE SUBDIAGONAL ELEMENTS ARE NONZERO.
  2447. C
  2448. C        B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS
  2449. C          HAVE BEEN ALTERED.  THE LOCATION B(N,1) IS USED TO STORE
  2450. C          EPS1 TIMES THE NORM OF B FOR LATER USE BY  QZVAL  AND  QZVEC.
  2451. C
  2452. C        Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS
  2453. C          (FOR BOTH STEPS) IF MATZ HAS BEEN SET TO .TRUE..
  2454. C
  2455. C        IERR IS SET TO
  2456. C          ZERO       FOR NORMAL RETURN,
  2457. C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
  2458. C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
  2459. C
  2460. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2461. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2462. C
  2463. C     THIS VERSION DATED AUGUST 1983.
  2464. C
  2465. C     ------------------------------------------------------------------
  2466. C
  2467.       SUBROUTINE QZVAL(NM,N,A,B,ALFR,ALFI,BETA,MATZ,Z)
  2468. C
  2469.       INTEGER I,J,N,EN,NA,NM,NN,ISW
  2470.       REAL A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N)
  2471.       REAL C,D,E,R,S,T,AN,A1,A2,BN,CQ,CZ,DI,DR,EI,TI,TR,U1,
  2472.      X       U2,V1,V2,A1I,A11,A12,A2I,A21,A22,B11,B12,B22,SQI,SQR,
  2473.      X       SSI,SSR,SZI,SZR,A11I,A11R,A12I,A12R,A22I,A22R,EPSB
  2474.       LOGICAL MATZ
  2475. C
  2476. C     THIS SUBROUTINE IS THE THIRD STEP OF THE QZ ALGORITHM
  2477. C     FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS,
  2478. C     SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART.
  2479. C
  2480. C     THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM
  2481. C     IN QUASI-TRIANGULAR FORM AND THE OTHER IN UPPER TRIANGULAR FORM.
  2482. C     IT REDUCES THE QUASI-TRIANGULAR MATRIX FURTHER, SO THAT ANY
  2483. C     REMAINING 2-BY-2 BLOCKS CORRESPOND TO PAIRS OF COMPLEX
  2484. C     EIGENVALUES, AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE
  2485. C     GENERALIZED EIGENVALUES.  IT IS USUALLY PRECEDED BY  QZHES
  2486. C     AND  QZIT  AND MAY BE FOLLOWED BY  QZVEC.
  2487. C
  2488. C     ON INPUT
  2489. C
  2490. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2491. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2492. C          DIMENSION STATEMENT.
  2493. C
  2494. C        N IS THE ORDER OF THE MATRICES.
  2495. C
  2496. C        A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX.
  2497. C
  2498. C        B CONTAINS A REAL UPPER TRIANGULAR MATRIX.  IN ADDITION,
  2499. C          LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB)
  2500. C          COMPUTED AND SAVED IN  QZIT.
  2501. C
  2502. C        MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS
  2503. C          ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING
  2504. C          EIGENVECTORS, AND TO .FALSE. OTHERWISE.
  2505. C
  2506. C        Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE
  2507. C          TRANSFORMATION MATRIX PRODUCED IN THE REDUCTIONS BY QZHES
  2508. C          AND QZIT, IF PERFORMED, OR ELSE THE IDENTITY MATRIX.
  2509. C          IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED.
  2510. C
  2511. C     ON OUTPUT
  2512. C
  2513. C        A HAS BEEN REDUCED FURTHER TO A QUASI-TRIANGULAR MATRIX
  2514. C          IN WHICH ALL NONZERO SUBDIAGONAL ELEMENTS CORRESPOND TO
  2515. C          PAIRS OF COMPLEX EIGENVALUES.
  2516. C
  2517. C        B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS
  2518. C          HAVE BEEN ALTERED.  B(N,1) IS UNALTERED.
  2519. C
  2520. C        ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE
  2521. C          DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX THAT WOULD BE
  2522. C          OBTAINED IF A WERE REDUCED COMPLETELY TO TRIANGULAR FORM
  2523. C          BY UNITARY TRANSFORMATIONS.  NON-ZERO VALUES OF ALFI OCCUR
  2524. C          IN PAIRS, THE FIRST MEMBER POSITIVE AND THE SECOND NEGATIVE.
  2525. C
  2526. C        BETA CONTAINS THE DIAGONAL ELEMENTS OF THE CORRESPONDING B,
  2527. C          NORMALIZED TO BE REAL AND NON-NEGATIVE.  THE GENERALIZED
  2528. C          EIGENVALUES ARE THEN THE RATIOS ((ALFR+I*ALFI)/BETA).
  2529. C
  2530. C        Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS
  2531. C          (FOR ALL THREE STEPS) IF MATZ HAS BEEN SET TO .TRUE.
  2532. C
  2533. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2534. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2535. C
  2536. C     THIS VERSION DATED AUGUST 1983.
  2537. C
  2538. C     ------------------------------------------------------------------
  2539. C
  2540.       SUBROUTINE QZVEC(NM,N,A,B,ALFR,ALFI,BETA,Z)
  2541. C
  2542.       INTEGER I,J,K,M,N,EN,II,JJ,NA,NM,NN,ISW,ENM2
  2543.       REAL A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N)
  2544.       REAL D,Q,R,S,T,W,X,Y,DI,DR,RA,RR,SA,TI,TR,T1,T2,W1,X1,
  2545.      X       ZZ,Z1,ALFM,ALMI,ALMR,BETM,EPSB
  2546. C
  2547. C     THIS SUBROUTINE IS THE OPTIONAL FOURTH STEP OF THE QZ ALGORITHM
  2548. C     FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS,
  2549. C     SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART.
  2550. C
  2551. C     THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM IN
  2552. C     QUASI-TRIANGULAR FORM (IN WHICH EACH 2-BY-2 BLOCK CORRESPONDS TO
  2553. C     A PAIR OF COMPLEX EIGENVALUES) AND THE OTHER IN UPPER TRIANGULAR
  2554. C     FORM.  IT COMPUTES THE EIGENVECTORS OF THE TRIANGULAR PROBLEM AND
  2555. C     TRANSFORMS THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM.
  2556. C     IT IS USUALLY PRECEDED BY  QZHES,  QZIT, AND  QZVAL.
  2557. C
  2558. C     ON INPUT
  2559. C
  2560. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2561. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2562. C          DIMENSION STATEMENT.
  2563. C
  2564. C        N IS THE ORDER OF THE MATRICES.
  2565. C
  2566. C        A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX.
  2567. C
  2568. C        B CONTAINS A REAL UPPER TRIANGULAR MATRIX.  IN ADDITION,
  2569. C          LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB)
  2570. C          COMPUTED AND SAVED IN  QZIT.
  2571. C
  2572. C        ALFR, ALFI, AND BETA  ARE VECTORS WITH COMPONENTS WHOSE
  2573. C          RATIOS ((ALFR+I*ALFI)/BETA) ARE THE GENERALIZED
  2574. C          EIGENVALUES.  THEY ARE USUALLY OBTAINED FROM  QZVAL.
  2575. C
  2576. C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
  2577. C          REDUCTIONS BY  QZHES,  QZIT, AND  QZVAL, IF PERFORMED.
  2578. C          IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ARE
  2579. C          DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX.
  2580. C
  2581. C     ON OUTPUT
  2582. C
  2583. C        A IS UNALTERED.  ITS SUBDIAGONAL ELEMENTS PROVIDE INFORMATION
  2584. C           ABOUT THE STORAGE OF THE COMPLEX EIGENVECTORS.
  2585. C
  2586. C        B HAS BEEN DESTROYED.
  2587. C
  2588. C        ALFR, ALFI, AND BETA ARE UNALTERED.
  2589. C
  2590. C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS.
  2591. C          IF ALFI(I) .EQ. 0.0, THE I-TH EIGENVALUE IS REAL AND
  2592. C            THE I-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR.
  2593. C          IF ALFI(I) .NE. 0.0, THE I-TH EIGENVALUE IS COMPLEX.
  2594. C            IF ALFI(I) .GT. 0.0, THE EIGENVALUE IS THE FIRST OF
  2595. C              A COMPLEX PAIR AND THE I-TH AND (I+1)-TH COLUMNS
  2596. C              OF Z CONTAIN ITS EIGENVECTOR.
  2597. C            IF ALFI(I) .LT. 0.0, THE EIGENVALUE IS THE SECOND OF
  2598. C              A COMPLEX PAIR AND THE (I-1)-TH AND I-TH COLUMNS
  2599. C              OF Z CONTAIN THE CONJUGATE OF ITS EIGENVECTOR.
  2600. C          EACH EIGENVECTOR IS NORMALIZED SO THAT THE MODULUS
  2601. C          OF ITS LARGEST COMPONENT IS 1.0 .
  2602. C
  2603. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2604. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2605. C
  2606. C     THIS VERSION DATED AUGUST 1983.
  2607. C
  2608. C     ------------------------------------------------------------------
  2609. C
  2610.       SUBROUTINE RATQR(N,EPS1,D,E,E2,M,W,IND,BD,TYPE,IDEF,IERR)
  2611. C
  2612.       INTEGER I,J,K,M,N,II,JJ,K1,IDEF,IERR,JDEF
  2613.       REAL D(N),E(N),E2(N),W(N),BD(N)
  2614.       REAL F,P,Q,R,S,EP,QP,ERR,TOT,EPS1,DELTA,EPSLON
  2615.       INTEGER IND(N)
  2616.       LOGICAL TYPE
  2617. C
  2618. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE RATQR,
  2619. C     NUM. MATH. 11, 264-272(1968) BY REINSCH AND BAUER.
  2620. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 257-265(1971).
  2621. C
  2622. C     THIS SUBROUTINE FINDS THE ALGEBRAICALLY SMALLEST OR LARGEST
  2623. C     EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE
  2624. C     RATIONAL QR METHOD WITH NEWTON CORRECTIONS.
  2625. C
  2626. C     ON INPUT
  2627. C
  2628. C        N IS THE ORDER OF THE MATRIX.
  2629. C
  2630. C        EPS1 IS A THEORETICAL ABSOLUTE ERROR TOLERANCE FOR THE
  2631. C          COMPUTED EIGENVALUES.  IF THE INPUT EPS1 IS NON-POSITIVE,
  2632. C          OR INDEED SMALLER THAN ITS DEFAULT VALUE, IT IS RESET
  2633. C          AT EACH ITERATION TO THE RESPECTIVE DEFAULT VALUE,
  2634. C          NAMELY, THE PRODUCT OF THE RELATIVE MACHINE PRECISION
  2635. C          AND THE MAGNITUDE OF THE CURRENT EIGENVALUE ITERATE.
  2636. C          THE THEORETICAL ABSOLUTE ERROR IN THE K-TH EIGENVALUE
  2637. C          IS USUALLY NOT GREATER THAN K TIMES EPS1.
  2638. C
  2639. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  2640. C
  2641. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  2642. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  2643. C
  2644. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  2645. C          E2(1) IS ARBITRARY.
  2646. C
  2647. C        M IS THE NUMBER OF EIGENVALUES TO BE FOUND.
  2648. C
  2649. C        IDEF SHOULD BE SET TO 1 IF THE INPUT MATRIX IS KNOWN TO BE
  2650. C          POSITIVE DEFINITE, TO -1 IF THE INPUT MATRIX IS KNOWN TO
  2651. C          BE NEGATIVE DEFINITE, AND TO 0 OTHERWISE.
  2652. C
  2653. C        TYPE SHOULD BE SET TO .TRUE. IF THE SMALLEST EIGENVALUES
  2654. C          ARE TO BE FOUND, AND TO .FALSE. IF THE LARGEST EIGENVALUES
  2655. C          ARE TO BE FOUND.
  2656. C
  2657. C     ON OUTPUT
  2658. C
  2659. C        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS
  2660. C          (LAST) DEFAULT VALUE.
  2661. C
  2662. C        D AND E ARE UNALTERED (UNLESS W OVERWRITES D).
  2663. C
  2664. C        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
  2665. C          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
  2666. C          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
  2667. C          E2(1) IS SET TO 0.0E0 IF THE SMALLEST EIGENVALUES HAVE BEEN
  2668. C          FOUND, AND TO 2.0E0 IF THE LARGEST EIGENVALUES HAVE BEEN
  2669. C          FOUND.  E2 IS OTHERWISE UNALTERED (UNLESS OVERWRITTEN BY BD).
  2670. C
  2671. C        W CONTAINS THE M ALGEBRAICALLY SMALLEST EIGENVALUES IN
  2672. C          ASCENDING ORDER, OR THE M LARGEST EIGENVALUES IN
  2673. C          DESCENDING ORDER.  IF AN ERROR EXIT IS MADE BECAUSE OF
  2674. C          AN INCORRECT SPECIFICATION OF IDEF, NO EIGENVALUES
  2675. C          ARE FOUND.  IF THE NEWTON ITERATES FOR A PARTICULAR
  2676. C          EIGENVALUE ARE NOT MONOTONE, THE BEST ESTIMATE OBTAINED
  2677. C          IS RETURNED AND IERR IS SET.  W MAY COINCIDE WITH D.
  2678. C
  2679. C        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
  2680. C          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
  2681. C          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
  2682. C          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC..
  2683. C
  2684. C        BD CONTAINS REFINED BOUNDS FOR THE THEORETICAL ERRORS OF THE
  2685. C          CORRESPONDING EIGENVALUES IN W.  THESE BOUNDS ARE USUALLY
  2686. C          WITHIN THE TOLERANCE SPECIFIED BY EPS1.  BD MAY COINCIDE
  2687. C          WITH E2.
  2688. C
  2689. C        IERR IS SET TO
  2690. C          ZERO       FOR NORMAL RETURN,
  2691. C          6*N+1      IF  IDEF  IS SET TO 1 AND  TYPE  TO .TRUE.
  2692. C                     WHEN THE MATRIX IS NOT POSITIVE DEFINITE, OR
  2693. C                     IF  IDEF  IS SET TO -1 AND  TYPE  TO .FALSE.
  2694. C                     WHEN THE MATRIX IS NOT NEGATIVE DEFINITE,
  2695. C          5*N+K      IF SUCCESSIVE ITERATES TO THE K-TH EIGENVALUE
  2696. C                     ARE NOT MONOTONE INCREASING, WHERE K REFERS
  2697. C                     TO THE LAST SUCH OCCURRENCE.
  2698. C
  2699. C     NOTE THAT SUBROUTINE TRIDIB IS GENERALLY FASTER AND MORE
  2700. C     ACCURATE THAN RATQR IF THE EIGENVALUES ARE CLUSTERED.
  2701. C
  2702. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2703. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2704. C
  2705. C     THIS VERSION DATED AUGUST 1983.
  2706. C
  2707. C     ------------------------------------------------------------------
  2708. C
  2709.       SUBROUTINE REBAK(NM,N,B,DL,M,Z)
  2710. C
  2711.       INTEGER I,J,K,M,N,I1,II,NM
  2712.       REAL B(NM,N),DL(N),Z(NM,M)
  2713.       REAL X
  2714. C
  2715. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE REBAKA,
  2716. C     NUM. MATH. 11, 99-110(1968) BY MARTIN AND WILKINSON.
  2717. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971).
  2718. C
  2719. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A GENERALIZED
  2720. C     SYMMETRIC EIGENSYSTEM BY BACK TRANSFORMING THOSE OF THE
  2721. C     DERIVED SYMMETRIC MATRIX DETERMINED BY  REDUC.
  2722. C
  2723. C     ON INPUT
  2724. C
  2725. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2726. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2727. C          DIMENSION STATEMENT.
  2728. C
  2729. C        N IS THE ORDER OF THE MATRIX SYSTEM.
  2730. C
  2731. C        B CONTAINS INFORMATION ABOUT THE SIMILARITY TRANSFORMATION
  2732. C          (CHOLESKY DECOMPOSITION) USED IN THE REDUCTION BY  REDUC
  2733. C          IN ITS STRICT LOWER TRIANGLE.
  2734. C
  2735. C        DL CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATION.
  2736. C
  2737. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
  2738. C
  2739. C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
  2740. C          IN ITS FIRST M COLUMNS.
  2741. C
  2742. C     ON OUTPUT
  2743. C
  2744. C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
  2745. C          IN ITS FIRST M COLUMNS.
  2746. C
  2747. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2748. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2749. C
  2750. C     THIS VERSION DATED AUGUST 1983.
  2751. C
  2752. C     ------------------------------------------------------------------
  2753. C
  2754.       SUBROUTINE REBAKB(NM,N,B,DL,M,Z)
  2755. C
  2756.       INTEGER I,J,K,M,N,I1,II,NM
  2757.       REAL B(NM,N),DL(N),Z(NM,M)
  2758.       REAL X
  2759. C
  2760. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE REBAKB,
  2761. C     NUM. MATH. 11, 99-110(1968) BY MARTIN AND WILKINSON.
  2762. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971).
  2763. C
  2764. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A GENERALIZED
  2765. C     SYMMETRIC EIGENSYSTEM BY BACK TRANSFORMING THOSE OF THE
  2766. C     DERIVED SYMMETRIC MATRIX DETERMINED BY  REDUC2.
  2767. C
  2768. C     ON INPUT
  2769. C
  2770. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2771. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2772. C          DIMENSION STATEMENT.
  2773. C
  2774. C        N IS THE ORDER OF THE MATRIX SYSTEM.
  2775. C
  2776. C        B CONTAINS INFORMATION ABOUT THE SIMILARITY TRANSFORMATION
  2777. C          (CHOLESKY DECOMPOSITION) USED IN THE REDUCTION BY  REDUC2
  2778. C          IN ITS STRICT LOWER TRIANGLE.
  2779. C
  2780. C        DL CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATION.
  2781. C
  2782. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
  2783. C
  2784. C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
  2785. C          IN ITS FIRST M COLUMNS.
  2786. C
  2787. C     ON OUTPUT
  2788. C
  2789. C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
  2790. C          IN ITS FIRST M COLUMNS.
  2791. C
  2792. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2793. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2794. C
  2795. C     THIS VERSION DATED AUGUST 1983.
  2796. C
  2797. C     ------------------------------------------------------------------
  2798. C
  2799.       SUBROUTINE REDUC(NM,N,A,B,DL,IERR)
  2800. C
  2801.       INTEGER I,J,K,N,I1,J1,NM,NN,IERR
  2802.       REAL A(NM,N),B(NM,N),DL(N)
  2803.       REAL X,Y
  2804. C
  2805. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE REDUC1,
  2806. C     NUM. MATH. 11, 99-110(1968) BY MARTIN AND WILKINSON.
  2807. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971).
  2808. C
  2809. C     THIS SUBROUTINE REDUCES THE GENERALIZED SYMMETRIC EIGENPROBLEM
  2810. C     AX=(LAMBDA)BX, WHERE B IS POSITIVE DEFINITE, TO THE STANDARD
  2811. C     SYMMETRIC EIGENPROBLEM USING THE CHOLESKY FACTORIZATION OF B.
  2812. C
  2813. C     ON INPUT
  2814. C
  2815. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2816. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2817. C          DIMENSION STATEMENT.
  2818. C
  2819. C        N IS THE ORDER OF THE MATRICES A AND B.  IF THE CHOLESKY
  2820. C          FACTOR L OF B IS ALREADY AVAILABLE, N SHOULD BE PREFIXED
  2821. C          WITH A MINUS SIGN.
  2822. C
  2823. C        A AND B CONTAIN THE REAL SYMMETRIC INPUT MATRICES.  ONLY THE
  2824. C          FULL UPPER TRIANGLES OF THE MATRICES NEED BE SUPPLIED.  IF
  2825. C          N IS NEGATIVE, THE STRICT LOWER TRIANGLE OF B CONTAINS,
  2826. C          INSTEAD, THE STRICT LOWER TRIANGLE OF ITS CHOLESKY FACTOR L.
  2827. C
  2828. C        DL CONTAINS, IF N IS NEGATIVE, THE DIAGONAL ELEMENTS OF L.
  2829. C
  2830. C     ON OUTPUT
  2831. C
  2832. C        A CONTAINS IN ITS FULL LOWER TRIANGLE THE FULL LOWER TRIANGLE
  2833. C          OF THE SYMMETRIC MATRIX DERIVED FROM THE REDUCTION TO THE
  2834. C          STANDARD FORM.  THE STRICT UPPER TRIANGLE OF A IS UNALTERED.
  2835. C
  2836. C        B CONTAINS IN ITS STRICT LOWER TRIANGLE THE STRICT LOWER
  2837. C          TRIANGLE OF ITS CHOLESKY FACTOR L.  THE FULL UPPER
  2838. C          TRIANGLE OF B IS UNALTERED.
  2839. C
  2840. C        DL CONTAINS THE DIAGONAL ELEMENTS OF L.
  2841. C
  2842. C        IERR IS SET TO
  2843. C          ZERO       FOR NORMAL RETURN,
  2844. C          7*N+1      IF B IS NOT POSITIVE DEFINITE.
  2845. C
  2846. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2847. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2848. C
  2849. C     THIS VERSION DATED AUGUST 1983.
  2850. C
  2851. C     ------------------------------------------------------------------
  2852. C
  2853.       SUBROUTINE REDUC2(NM,N,A,B,DL,IERR)
  2854. C
  2855.       INTEGER I,J,K,N,I1,J1,NM,NN,IERR
  2856.       REAL A(NM,N),B(NM,N),DL(N)
  2857.       REAL X,Y
  2858. C
  2859. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE REDUC2,
  2860. C     NUM. MATH. 11, 99-110(1968) BY MARTIN AND WILKINSON.
  2861. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971).
  2862. C
  2863. C     THIS SUBROUTINE REDUCES THE GENERALIZED SYMMETRIC EIGENPROBLEMS
  2864. C     ABX=(LAMBDA)X OR BAY=(LAMBDA)Y, WHERE B IS POSITIVE DEFINITE,
  2865. C     TO THE STANDARD SYMMETRIC EIGENPROBLEM USING THE CHOLESKY
  2866. C     FACTORIZATION OF B.
  2867. C
  2868. C     ON INPUT
  2869. C
  2870. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2871. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2872. C          DIMENSION STATEMENT.
  2873. C
  2874. C        N IS THE ORDER OF THE MATRICES A AND B.  IF THE CHOLESKY
  2875. C          FACTOR L OF B IS ALREADY AVAILABLE, N SHOULD BE PREFIXED
  2876. C          WITH A MINUS SIGN.
  2877. C
  2878. C        A AND B CONTAIN THE REAL SYMMETRIC INPUT MATRICES.  ONLY THE
  2879. C          FULL UPPER TRIANGLES OF THE MATRICES NEED BE SUPPLIED.  IF
  2880. C          N IS NEGATIVE, THE STRICT LOWER TRIANGLE OF B CONTAINS,
  2881. C          INSTEAD, THE STRICT LOWER TRIANGLE OF ITS CHOLESKY FACTOR L.
  2882. C
  2883. C        DL CONTAINS, IF N IS NEGATIVE, THE DIAGONAL ELEMENTS OF L.
  2884. C
  2885. C     ON OUTPUT
  2886. C
  2887. C        A CONTAINS IN ITS FULL LOWER TRIANGLE THE FULL LOWER TRIANGLE
  2888. C          OF THE SYMMETRIC MATRIX DERIVED FROM THE REDUCTION TO THE
  2889. C          STANDARD FORM.  THE STRICT UPPER TRIANGLE OF A IS UNALTERED.
  2890. C
  2891. C        B CONTAINS IN ITS STRICT LOWER TRIANGLE THE STRICT LOWER
  2892. C          TRIANGLE OF ITS CHOLESKY FACTOR L.  THE FULL UPPER
  2893. C          TRIANGLE OF B IS UNALTERED.
  2894. C
  2895. C        DL CONTAINS THE DIAGONAL ELEMENTS OF L.
  2896. C
  2897. C        IERR IS SET TO
  2898. C          ZERO       FOR NORMAL RETURN,
  2899. C          7*N+1      IF B IS NOT POSITIVE DEFINITE.
  2900. C
  2901. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2902. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2903. C
  2904. C     THIS VERSION DATED AUGUST 1983.
  2905. C
  2906. C     ------------------------------------------------------------------
  2907. C
  2908.       SUBROUTINE RG(NM,N,A,WR,WI,MATZ,Z,IV1,FV1,IERR)
  2909. C
  2910.       INTEGER N,NM,IS1,IS2,IERR,MATZ
  2911.       REAL A(NM,N),WR(N),WI(N),Z(NM,N),FV1(N)
  2912.       INTEGER IV1(N)
  2913. C
  2914. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  2915. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  2916. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  2917. C     OF A REAL GENERAL MATRIX.
  2918. C
  2919. C     ON INPUT
  2920. C
  2921. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  2922. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2923. C        DIMENSION STATEMENT.
  2924. C
  2925. C        N  IS THE ORDER OF THE MATRIX  A.
  2926. C
  2927. C        A  CONTAINS THE REAL GENERAL MATRIX.
  2928. C
  2929. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  2930. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  2931. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  2932. C
  2933. C     ON OUTPUT
  2934. C
  2935. C        WR  AND  WI  CONTAIN THE REAL AND IMAGINARY PARTS,
  2936. C        RESPECTIVELY, OF THE EIGENVALUES.  COMPLEX CONJUGATE
  2937. C        PAIRS OF EIGENVALUES APPEAR CONSECUTIVELY WITH THE
  2938. C        EIGENVALUE HAVING THE POSITIVE IMAGINARY PART FIRST.
  2939. C
  2940. C        Z  CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS
  2941. C        IF MATZ IS NOT ZERO.  IF THE J-TH EIGENVALUE IS REAL, THE
  2942. C        J-TH COLUMN OF  Z  CONTAINS ITS EIGENVECTOR.  IF THE J-TH
  2943. C        EIGENVALUE IS COMPLEX WITH POSITIVE IMAGINARY PART, THE
  2944. C        J-TH AND (J+1)-TH COLUMNS OF  Z  CONTAIN THE REAL AND
  2945. C        IMAGINARY PARTS OF ITS EIGENVECTOR.  THE CONJUGATE OF THIS
  2946. C        VECTOR IS THE EIGENVECTOR FOR THE CONJUGATE EIGENVALUE.
  2947. C
  2948. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  2949. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR HQR
  2950. C           AND HQR2.  THE NORMAL COMPLETION CODE IS ZERO.
  2951. C
  2952. C        IV1  AND  FV1  ARE TEMPORARY STORAGE ARRAYS.
  2953. C
  2954. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  2955. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  2956. C
  2957. C     THIS VERSION DATED AUGUST 1983.
  2958. C
  2959. C     ------------------------------------------------------------------
  2960. C
  2961.       SUBROUTINE RGG(NM,N,A,B,ALFR,ALFI,BETA,MATZ,Z,IERR)
  2962. C
  2963.       INTEGER N,NM,IERR,MATZ
  2964.       REAL A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N)
  2965.       LOGICAL TF
  2966. C
  2967. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  2968. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  2969. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  2970. C     FOR THE REAL GENERAL GENERALIZED EIGENPROBLEM  AX = (LAMBDA)BX.
  2971. C
  2972. C     ON INPUT
  2973. C
  2974. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  2975. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2976. C        DIMENSION STATEMENT.
  2977. C
  2978. C        N  IS THE ORDER OF THE MATRICES  A  AND  B.
  2979. C
  2980. C        A  CONTAINS A REAL GENERAL MATRIX.
  2981. C
  2982. C        B  CONTAINS A REAL GENERAL MATRIX.
  2983. C
  2984. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  2985. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  2986. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  2987. C
  2988. C     ON OUTPUT
  2989. C
  2990. C        ALFR  AND  ALFI  CONTAIN THE REAL AND IMAGINARY PARTS,
  2991. C        RESPECTIVELY, OF THE NUMERATORS OF THE EIGENVALUES.
  2992. C
  2993. C        BETA  CONTAINS THE DENOMINATORS OF THE EIGENVALUES,
  2994. C        WHICH ARE THUS GIVEN BY THE RATIOS  (ALFR+I*ALFI)/BETA.
  2995. C        COMPLEX CONJUGATE PAIRS OF EIGENVALUES APPEAR CONSECUTIVELY
  2996. C        WITH THE EIGENVALUE HAVING THE POSITIVE IMAGINARY PART FIRST.
  2997. C
  2998. C        Z  CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS
  2999. C        IF MATZ IS NOT ZERO.  IF THE J-TH EIGENVALUE IS REAL, THE
  3000. C        J-TH COLUMN OF  Z  CONTAINS ITS EIGENVECTOR.  IF THE J-TH
  3001. C        EIGENVALUE IS COMPLEX WITH POSITIVE IMAGINARY PART, THE
  3002. C        J-TH AND (J+1)-TH COLUMNS OF  Z  CONTAIN THE REAL AND
  3003. C        IMAGINARY PARTS OF ITS EIGENVECTOR.  THE CONJUGATE OF THIS
  3004. C        VECTOR IS THE EIGENVECTOR FOR THE CONJUGATE EIGENVALUE.
  3005. C
  3006. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3007. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR QZIT.
  3008. C           THE NORMAL COMPLETION CODE IS ZERO.
  3009. C
  3010. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3011. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3012. C
  3013. C     THIS VERSION DATED AUGUST 1983.
  3014. C
  3015. C     ------------------------------------------------------------------
  3016. C
  3017.       SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR)
  3018. C
  3019.       INTEGER N,NM,IERR,MATZ
  3020.       REAL A(NM,N),W(N),Z(NM,N),FV1(N),FV2(N)
  3021. C
  3022. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  3023. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  3024. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  3025. C     OF A REAL SYMMETRIC MATRIX.
  3026. C
  3027. C     ON INPUT
  3028. C
  3029. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  3030. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3031. C        DIMENSION STATEMENT.
  3032. C
  3033. C        N  IS THE ORDER OF THE MATRIX  A.
  3034. C
  3035. C        A  CONTAINS THE REAL SYMMETRIC MATRIX.
  3036. C
  3037. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  3038. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  3039. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  3040. C
  3041. C     ON OUTPUT
  3042. C
  3043. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
  3044. C
  3045. C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO.
  3046. C
  3047. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3048. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT
  3049. C           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO.
  3050. C
  3051. C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
  3052. C
  3053. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3054. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3055. C
  3056. C     THIS VERSION DATED AUGUST 1983.
  3057. C
  3058. C     ------------------------------------------------------------------
  3059. C
  3060.       SUBROUTINE RSB(NM,N,MB,A,W,MATZ,Z,FV1,FV2,IERR)
  3061. C
  3062.       INTEGER N,MB,NM,IERR,MATZ
  3063.       REAL A(NM,MB),W(N),Z(NM,N),FV1(N),FV2(N)
  3064.       LOGICAL TF
  3065. C
  3066. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  3067. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  3068. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  3069. C     OF A REAL SYMMETRIC BAND MATRIX.
  3070. C
  3071. C     ON INPUT
  3072. C
  3073. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  3074. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3075. C        DIMENSION STATEMENT.
  3076. C
  3077. C        N  IS THE ORDER OF THE MATRIX  A.
  3078. C
  3079. C        MB  IS THE HALF BAND WIDTH OF THE MATRIX, DEFINED AS THE
  3080. C        NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL
  3081. C        DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE
  3082. C        LOWER TRIANGLE OF THE MATRIX.
  3083. C
  3084. C        A  CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
  3085. C        BAND MATRIX.  ITS LOWEST SUBDIAGONAL IS STORED IN THE
  3086. C        LAST  N+1-MB  POSITIONS OF THE FIRST COLUMN, ITS NEXT
  3087. C        SUBDIAGONAL IN THE LAST  N+2-MB  POSITIONS OF THE
  3088. C        SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND
  3089. C        FINALLY ITS PRINCIPAL DIAGONAL IN THE  N  POSITIONS
  3090. C        OF THE LAST COLUMN.  CONTENTS OF STORAGES NOT PART
  3091. C        OF THE MATRIX ARE ARBITRARY.
  3092. C
  3093. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  3094. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  3095. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  3096. C
  3097. C     ON OUTPUT
  3098. C
  3099. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
  3100. C
  3101. C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO.
  3102. C
  3103. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3104. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT
  3105. C           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO.
  3106. C
  3107. C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
  3108. C
  3109. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3110. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3111. C
  3112. C     THIS VERSION DATED AUGUST 1983.
  3113. C
  3114. C     ------------------------------------------------------------------
  3115. C
  3116.       SUBROUTINE RSG(NM,N,A,B,W,MATZ,Z,FV1,FV2,IERR)
  3117. C
  3118.       INTEGER N,NM,IERR,MATZ
  3119.       REAL A(NM,N),B(NM,N),W(N),Z(NM,N),FV1(N),FV2(N)
  3120. C
  3121. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  3122. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  3123. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  3124. C     FOR THE REAL SYMMETRIC GENERALIZED EIGENPROBLEM  AX = (LAMBDA)BX.
  3125. C
  3126. C     ON INPUT
  3127. C
  3128. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  3129. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3130. C        DIMENSION STATEMENT.
  3131. C
  3132. C        N  IS THE ORDER OF THE MATRICES  A  AND  B.
  3133. C
  3134. C        A  CONTAINS A REAL SYMMETRIC MATRIX.
  3135. C
  3136. C        B  CONTAINS A POSITIVE DEFINITE REAL SYMMETRIC MATRIX.
  3137. C
  3138. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  3139. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  3140. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  3141. C
  3142. C     ON OUTPUT
  3143. C
  3144. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
  3145. C
  3146. C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO.
  3147. C
  3148. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3149. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT
  3150. C           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO.
  3151. C
  3152. C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
  3153. C
  3154. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3155. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3156. C
  3157. C     THIS VERSION DATED AUGUST 1983.
  3158. C
  3159. C     ------------------------------------------------------------------
  3160. C
  3161.       SUBROUTINE RSGAB(NM,N,A,B,W,MATZ,Z,FV1,FV2,IERR)
  3162. C
  3163.       INTEGER N,NM,IERR,MATZ
  3164.       REAL A(NM,N),B(NM,N),W(N),Z(NM,N),FV1(N),FV2(N)
  3165. C
  3166. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  3167. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  3168. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  3169. C     FOR THE REAL SYMMETRIC GENERALIZED EIGENPROBLEM  ABX = (LAMBDA)X.
  3170. C
  3171. C     ON INPUT
  3172. C
  3173. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  3174. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3175. C        DIMENSION STATEMENT.
  3176. C
  3177. C        N  IS THE ORDER OF THE MATRICES  A  AND  B.
  3178. C
  3179. C        A  CONTAINS A REAL SYMMETRIC MATRIX.
  3180. C
  3181. C        B  CONTAINS A POSITIVE DEFINITE REAL SYMMETRIC MATRIX.
  3182. C
  3183. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  3184. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  3185. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  3186. C
  3187. C     ON OUTPUT
  3188. C
  3189. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
  3190. C
  3191. C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO.
  3192. C
  3193. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3194. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT
  3195. C           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO.
  3196. C
  3197. C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
  3198. C
  3199. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3200. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3201. C
  3202. C     THIS VERSION DATED AUGUST 1983.
  3203. C
  3204. C     ------------------------------------------------------------------
  3205. C
  3206.       SUBROUTINE RSGBA(NM,N,A,B,W,MATZ,Z,FV1,FV2,IERR)
  3207. C
  3208.       INTEGER N,NM,IERR,MATZ
  3209.       REAL A(NM,N),B(NM,N),W(N),Z(NM,N),FV1(N),FV2(N)
  3210. C
  3211. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  3212. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  3213. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  3214. C     FOR THE REAL SYMMETRIC GENERALIZED EIGENPROBLEM  BAX = (LAMBDA)X.
  3215. C
  3216. C     ON INPUT
  3217. C
  3218. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  3219. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3220. C        DIMENSION STATEMENT.
  3221. C
  3222. C        N  IS THE ORDER OF THE MATRICES  A  AND  B.
  3223. C
  3224. C        A  CONTAINS A REAL SYMMETRIC MATRIX.
  3225. C
  3226. C        B  CONTAINS A POSITIVE DEFINITE REAL SYMMETRIC MATRIX.
  3227. C
  3228. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  3229. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  3230. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  3231. C
  3232. C     ON OUTPUT
  3233. C
  3234. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
  3235. C
  3236. C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO.
  3237. C
  3238. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3239. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT
  3240. C           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO.
  3241. C
  3242. C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
  3243. C
  3244. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3245. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3246. C
  3247. C     THIS VERSION DATED AUGUST 1983.
  3248. C
  3249. C     ------------------------------------------------------------------
  3250. C
  3251.       SUBROUTINE RSM(NM,N,A,W,M,Z,FWORK,IWORK,IERR)
  3252. C
  3253.       INTEGER N,NM,M,IWORK(N),IERR
  3254.       REAL A(NM,N),W(N),Z(NM,M),FWORK(1)
  3255. C
  3256. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  3257. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  3258. C     TO FIND ALL OF THE EIGENVALUES AND SOME OF THE EIGENVECTORS
  3259. C     OF A REAL SYMMETRIC MATRIX.
  3260. C
  3261. C     ON INPUT
  3262. C
  3263. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  3264. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3265. C        DIMENSION STATEMENT.
  3266. C
  3267. C        N  IS THE ORDER OF THE MATRIX  A.
  3268. C
  3269. C        A  CONTAINS THE REAL SYMMETRIC MATRIX.
  3270. C
  3271. C        M  THE EIGENVECTORS CORRESPONDING TO THE FIRST M EIGENVALUES
  3272. C           ARE TO BE COMPUTED.
  3273. C           IF M = 0 THEN NO EIGENVECTORS ARE COMPUTED.
  3274. C           IF M = N THEN ALL OF THE EIGENVECTORS ARE COMPUTED.
  3275. C
  3276. C     ON OUTPUT
  3277. C
  3278. C        W  CONTAINS ALL N EIGENVALUES IN ASCENDING ORDER.
  3279. C
  3280. C        Z  CONTAINS THE ORTHONORMAL EIGENVECTORS ASSOCIATED WITH
  3281. C           THE FIRST M EIGENVALUES.
  3282. C
  3283. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3284. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT,
  3285. C           IMTQLV AND TINVIT.  THE NORMAL COMPLETION CODE IS ZERO.
  3286. C
  3287. C        FWORK  IS A TEMPORARY STORAGE ARRAY OF DIMENSION 8*N.
  3288. C
  3289. C        IWORK  IS AN INTEGER TEMPORARY STORAGE ARRAY OF DIMENSION N.
  3290. C
  3291. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3292. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3293. C
  3294. C     THIS VERSION DATED AUGUST 1983.
  3295. C
  3296. C     ------------------------------------------------------------------
  3297. C
  3298.       SUBROUTINE RSP(NM,N,NV,A,W,MATZ,Z,FV1,FV2,IERR)
  3299. C
  3300.       INTEGER I,J,N,NM,NV,IERR,MATZ
  3301.       REAL A(NV),W(N),Z(NM,N),FV1(N),FV2(N)
  3302. C
  3303. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  3304. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  3305. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  3306. C     OF A REAL SYMMETRIC PACKED MATRIX.
  3307. C
  3308. C     ON INPUT
  3309. C
  3310. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  3311. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3312. C        DIMENSION STATEMENT.
  3313. C
  3314. C        N  IS THE ORDER OF THE MATRIX  A.
  3315. C
  3316. C        NV  IS AN INTEGER VARIABLE SET EQUAL TO THE
  3317. C        DIMENSION OF THE ARRAY  A  AS SPECIFIED FOR
  3318. C        A  IN THE CALLING PROGRAM.  NV  MUST NOT BE
  3319. C        LESS THAN  N*(N+1)/2.
  3320. C
  3321. C        A  CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
  3322. C        PACKED MATRIX STORED ROW-WISE.
  3323. C
  3324. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  3325. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  3326. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  3327. C
  3328. C     ON OUTPUT
  3329. C
  3330. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
  3331. C
  3332. C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO.
  3333. C
  3334. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3335. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT
  3336. C           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO.
  3337. C
  3338. C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
  3339. C
  3340. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3341. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3342. C
  3343. C     THIS VERSION DATED AUGUST 1983.
  3344. C
  3345. C     ------------------------------------------------------------------
  3346. C
  3347.       SUBROUTINE RST(NM,N,W,E,MATZ,Z,IERR)
  3348. C
  3349.       INTEGER I,J,N,NM,IERR,MATZ
  3350.       REAL W(N),E(N),Z(NM,N)
  3351. C
  3352. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  3353. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  3354. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  3355. C     OF A REAL SYMMETRIC TRIDIAGONAL MATRIX.
  3356. C
  3357. C     ON INPUT
  3358. C
  3359. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  3360. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3361. C        DIMENSION STATEMENT.
  3362. C
  3363. C        N  IS THE ORDER OF THE MATRIX.
  3364. C
  3365. C        W  CONTAINS THE DIAGONAL ELEMENTS OF THE REAL
  3366. C        SYMMETRIC TRIDIAGONAL MATRIX.
  3367. C
  3368. C        E  CONTAINS THE SUBDIAGONAL ELEMENTS OF THE MATRIX IN
  3369. C        ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  3370. C
  3371. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  3372. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  3373. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  3374. C
  3375. C     ON OUTPUT
  3376. C
  3377. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
  3378. C
  3379. C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO.
  3380. C
  3381. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3382. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR IMTQL1
  3383. C           AND IMTQL2.  THE NORMAL COMPLETION CODE IS ZERO.
  3384. C
  3385. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3386. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3387. C
  3388. C     THIS VERSION DATED AUGUST 1983.
  3389. C
  3390. C     ------------------------------------------------------------------
  3391. C
  3392.       SUBROUTINE RT(NM,N,A,W,MATZ,Z,FV1,IERR)
  3393. C
  3394.       INTEGER N,NM,IERR,MATZ
  3395.       REAL A(NM,3),W(N),Z(NM,N),FV1(N)
  3396. C
  3397. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  3398. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  3399. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  3400. C     OF A SPECIAL REAL TRIDIAGONAL MATRIX.
  3401. C
  3402. C     ON INPUT
  3403. C
  3404. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  3405. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3406. C        DIMENSION STATEMENT.
  3407. C
  3408. C        N  IS THE ORDER OF THE MATRIX  A.
  3409. C
  3410. C        A  CONTAINS THE SPECIAL REAL TRIDIAGONAL MATRIX IN ITS
  3411. C        FIRST THREE COLUMNS.  THE SUBDIAGONAL ELEMENTS ARE STORED
  3412. C        IN THE LAST  N-1  POSITIONS OF THE FIRST COLUMN, THE
  3413. C        DIAGONAL ELEMENTS IN THE SECOND COLUMN, AND THE SUPERDIAGONAL
  3414. C        ELEMENTS IN THE FIRST  N-1  POSITIONS OF THE THIRD COLUMN.
  3415. C        ELEMENTS  A(1,1)  AND  A(N,3)  ARE ARBITRARY.
  3416. C
  3417. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  3418. C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
  3419. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  3420. C
  3421. C     ON OUTPUT
  3422. C
  3423. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
  3424. C
  3425. C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO.
  3426. C
  3427. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
  3428. C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR IMTQL1
  3429. C           AND IMTQL2.  THE NORMAL COMPLETION CODE IS ZERO.
  3430. C
  3431. C        FV1  IS A TEMPORARY STORAGE ARRAY.
  3432. C
  3433. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3434. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3435. C
  3436. C     THIS VERSION DATED AUGUST 1983.
  3437. C
  3438. C     ------------------------------------------------------------------
  3439. C
  3440.       SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1)
  3441. C
  3442.       INTEGER I,J,K,L,M,N,II,I1,KK,K1,LL,L1,MN,NM,ITS,IERR
  3443.       REAL A(NM,N),W(N),U(NM,N),V(NM,N),RV1(N)
  3444.       REAL C,F,G,H,S,X,Y,Z,TST1,TST2,SCALE,PYTHAG
  3445.       LOGICAL MATU,MATV
  3446. C
  3447. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE SVD,
  3448. C     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH.
  3449. C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
  3450. C
  3451. C     THIS SUBROUTINE DETERMINES THE SINGULAR VALUE DECOMPOSITION
  3452. C          T
  3453. C     A=USV  OF A REAL M BY N RECTANGULAR MATRIX.  HOUSEHOLDER
  3454. C     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED.
  3455. C
  3456. C     ON INPUT
  3457. C
  3458. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  3459. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3460. C          DIMENSION STATEMENT.  NOTE THAT NM MUST BE AT LEAST
  3461. C          AS LARGE AS THE MAXIMUM OF M AND N.
  3462. C
  3463. C        M IS THE NUMBER OF ROWS OF A (AND U).
  3464. C
  3465. C        N IS THE NUMBER OF COLUMNS OF A (AND U) AND THE ORDER OF V.
  3466. C
  3467. C        A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED.
  3468. C
  3469. C        MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE
  3470. C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
  3471. C
  3472. C        MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE
  3473. C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
  3474. C
  3475. C     ON OUTPUT
  3476. C
  3477. C        A IS UNALTERED (UNLESS OVERWRITTEN BY U OR V).
  3478. C
  3479. C        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE
  3480. C          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN
  3481. C          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT
  3482. C          FOR INDICES IERR+1,IERR+2,...,N.
  3483. C
  3484. C        U CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE
  3485. C          DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE.  OTHERWISE
  3486. C          U IS USED AS A TEMPORARY ARRAY.  U MAY COINCIDE WITH A.
  3487. C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING
  3488. C          TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT.
  3489. C
  3490. C        V CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF
  3491. C          MATV HAS BEEN SET TO .TRUE.  OTHERWISE V IS NOT REFERENCED.
  3492. C          V MAY ALSO COINCIDE WITH A IF U IS NOT NEEDED.  IF AN ERROR
  3493. C          EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO INDICES OF
  3494. C          CORRECT SINGULAR VALUES SHOULD BE CORRECT.
  3495. C
  3496. C        IERR IS SET TO
  3497. C          ZERO       FOR NORMAL RETURN,
  3498. C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
  3499. C                     DETERMINED AFTER 30 ITERATIONS.
  3500. C
  3501. C        RV1 IS A TEMPORARY STORAGE ARRAY.
  3502. C
  3503. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  3504. C
  3505. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3506. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3507. C
  3508. C     THIS VERSION DATED AUGUST 1983.
  3509. C
  3510. C     ------------------------------------------------------------------
  3511. C
  3512.       SUBROUTINE TINVIT(NM,N,D,E,E2,M,W,IND,Z,
  3513.      X                  IERR,RV1,RV2,RV3,RV4,RV6)
  3514. C
  3515.       INTEGER I,J,M,N,P,Q,R,S,II,IP,JJ,NM,ITS,TAG,IERR,GROUP
  3516.       REAL D(N),E(N),E2(N),W(M),Z(NM,M),
  3517.      X       RV1(N),RV2(N),RV3(N),RV4(N),RV6(N)
  3518.       REAL U,V,UK,XU,X0,X1,EPS2,EPS3,EPS4,NORM,ORDER,EPSLON,
  3519.      X       PYTHAG
  3520.       INTEGER IND(M)
  3521. C
  3522. C     THIS SUBROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH-
  3523. C     NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
  3524. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
  3525. C
  3526. C     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
  3527. C     SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
  3528. C     USING INVERSE ITERATION.
  3529. C
  3530. C     ON INPUT
  3531. C
  3532. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  3533. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3534. C          DIMENSION STATEMENT.
  3535. C
  3536. C        N IS THE ORDER OF THE MATRIX.
  3537. C
  3538. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  3539. C
  3540. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  3541. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  3542. C
  3543. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E,
  3544. C          WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
  3545. C          E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
  3546. C          THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM
  3547. C          OF THE MAGNITUDES OF D(I) AND D(I-1).  E2(1) MUST CONTAIN
  3548. C          0.0E0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0E0
  3549. C          IF THE EIGENVALUES ARE IN DESCENDING ORDER.  IF  BISECT,
  3550. C          TRIDIB, OR  IMTQLV  HAS BEEN USED TO FIND THE EIGENVALUES,
  3551. C          THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE.
  3552. C
  3553. C        M IS THE NUMBER OF SPECIFIED EIGENVALUES.
  3554. C
  3555. C        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
  3556. C
  3557. C        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
  3558. C          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
  3559. C          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
  3560. C          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.
  3561. C
  3562. C     ON OUTPUT
  3563. C
  3564. C        ALL INPUT ARRAYS ARE UNALTERED.
  3565. C
  3566. C        Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
  3567. C          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO.
  3568. C
  3569. C        IERR IS SET TO
  3570. C          ZERO       FOR NORMAL RETURN,
  3571. C          -R         IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
  3572. C                     EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS.
  3573. C
  3574. C        RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
  3575. C
  3576. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  3577. C
  3578. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3579. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3580. C
  3581. C     THIS VERSION DATED AUGUST 1983.
  3582. C
  3583. C     ------------------------------------------------------------------
  3584. C
  3585.       SUBROUTINE TQL1(N,D,E,IERR)
  3586. C
  3587.       INTEGER I,J,L,M,N,II,L1,L2,MML,IERR
  3588.       REAL D(N),E(N)
  3589.       REAL C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2,TST1,TST2,PYTHAG
  3590. C
  3591. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL1,
  3592. C     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
  3593. C     WILKINSON.
  3594. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
  3595. C
  3596. C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
  3597. C     TRIDIAGONAL MATRIX BY THE QL METHOD.
  3598. C
  3599. C     ON INPUT
  3600. C
  3601. C        N IS THE ORDER OF THE MATRIX.
  3602. C
  3603. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  3604. C
  3605. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  3606. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  3607. C
  3608. C      ON OUTPUT
  3609. C
  3610. C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
  3611. C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
  3612. C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
  3613. C          THE SMALLEST EIGENVALUES.
  3614. C
  3615. C        E HAS BEEN DESTROYED.
  3616. C
  3617. C        IERR IS SET TO
  3618. C          ZERO       FOR NORMAL RETURN,
  3619. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
  3620. C                     DETERMINED AFTER 30 ITERATIONS.
  3621. C
  3622. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  3623. C
  3624. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3625. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3626. C
  3627. C     THIS VERSION DATED AUGUST 1983.
  3628. C
  3629. C     ------------------------------------------------------------------
  3630. C
  3631.       SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
  3632. C
  3633.       INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR
  3634.       REAL D(N),E(N),Z(NM,N)
  3635.       REAL C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2,TST1,TST2,PYTHAG
  3636. C
  3637. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
  3638. C     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
  3639. C     WILKINSON.
  3640. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
  3641. C
  3642. C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
  3643. C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
  3644. C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
  3645. C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
  3646. C     FULL MATRIX TO TRIDIAGONAL FORM.
  3647. C
  3648. C     ON INPUT
  3649. C
  3650. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  3651. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3652. C          DIMENSION STATEMENT.
  3653. C
  3654. C        N IS THE ORDER OF THE MATRIX.
  3655. C
  3656. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  3657. C
  3658. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  3659. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  3660. C
  3661. C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
  3662. C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
  3663. C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
  3664. C          THE IDENTITY MATRIX.
  3665. C
  3666. C      ON OUTPUT
  3667. C
  3668. C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
  3669. C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
  3670. C          UNORDERED FOR INDICES 1,2,...,IERR-1.
  3671. C
  3672. C        E HAS BEEN DESTROYED.
  3673. C
  3674. C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
  3675. C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
  3676. C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
  3677. C          EIGENVALUES.
  3678. C
  3679. C        IERR IS SET TO
  3680. C          ZERO       FOR NORMAL RETURN,
  3681. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
  3682. C                     DETERMINED AFTER 30 ITERATIONS.
  3683. C
  3684. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  3685. C
  3686. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3687. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3688. C
  3689. C     THIS VERSION DATED AUGUST 1983.
  3690. C
  3691. C     ------------------------------------------------------------------
  3692. C
  3693.       SUBROUTINE TQLRAT(N,D,E2,IERR)
  3694. C
  3695.       INTEGER I,J,L,M,N,II,L1,MML,IERR
  3696.       REAL D(N),E2(N)
  3697.       REAL B,C,F,G,H,P,R,S,T,EPSLON,PYTHAG
  3698. C
  3699. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
  3700. C     ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
  3701. C
  3702. C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
  3703. C     TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD.
  3704. C
  3705. C     ON INPUT
  3706. C
  3707. C        N IS THE ORDER OF THE MATRIX.
  3708. C
  3709. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  3710. C
  3711. C        E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE
  3712. C          INPUT MATRIX IN ITS LAST N-1 POSITIONS.  E2(1) IS ARBITRARY.
  3713. C
  3714. C      ON OUTPUT
  3715. C
  3716. C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
  3717. C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
  3718. C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
  3719. C          THE SMALLEST EIGENVALUES.
  3720. C
  3721. C        E2 HAS BEEN DESTROYED.
  3722. C
  3723. C        IERR IS SET TO
  3724. C          ZERO       FOR NORMAL RETURN,
  3725. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
  3726. C                     DETERMINED AFTER 30 ITERATIONS.
  3727. C
  3728. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  3729. C
  3730. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3731. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3732. C
  3733. C     THIS VERSION DATED AUGUST 1983.
  3734. C
  3735. C     ------------------------------------------------------------------
  3736. C
  3737.       SUBROUTINE TRBAK1(NM,N,A,E,M,Z)
  3738. C
  3739.       INTEGER I,J,K,L,M,N,NM
  3740.       REAL A(NM,N),E(N),Z(NM,M)
  3741.       REAL S
  3742. C
  3743. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK1,
  3744. C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
  3745. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  3746. C
  3747. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
  3748. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  3749. C     SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  TRED1.
  3750. C
  3751. C     ON INPUT
  3752. C
  3753. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  3754. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3755. C          DIMENSION STATEMENT.
  3756. C
  3757. C        N IS THE ORDER OF THE MATRIX.
  3758. C
  3759. C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
  3760. C          FORMATIONS USED IN THE REDUCTION BY  TRED1
  3761. C          IN ITS STRICT LOWER TRIANGLE.
  3762. C
  3763. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  3764. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  3765. C
  3766. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
  3767. C
  3768. C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
  3769. C          IN ITS FIRST M COLUMNS.
  3770. C
  3771. C     ON OUTPUT
  3772. C
  3773. C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
  3774. C          IN ITS FIRST M COLUMNS.
  3775. C
  3776. C     NOTE THAT TRBAK1 PRESERVES VECTOR EUCLIDEAN NORMS.
  3777. C
  3778. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3779. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3780. C
  3781. C     THIS VERSION DATED AUGUST 1983.
  3782. C
  3783. C     ------------------------------------------------------------------
  3784. C
  3785.       SUBROUTINE TRBAK3(NM,N,NV,A,M,Z)
  3786. C
  3787.       INTEGER I,J,K,L,M,N,IK,IZ,NM,NV
  3788.       REAL A(NV),Z(NM,M)
  3789.       REAL H,S
  3790. C
  3791. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
  3792. C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
  3793. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  3794. C
  3795. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
  3796. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
  3797. C     SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  TRED3.
  3798. C
  3799. C     ON INPUT
  3800. C
  3801. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  3802. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3803. C          DIMENSION STATEMENT.
  3804. C
  3805. C        N IS THE ORDER OF THE MATRIX.
  3806. C
  3807. C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
  3808. C          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  3809. C
  3810. C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS
  3811. C          USED IN THE REDUCTION BY  TRED3  IN ITS FIRST
  3812. C          N*(N+1)/2 POSITIONS.
  3813. C
  3814. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
  3815. C
  3816. C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
  3817. C          IN ITS FIRST M COLUMNS.
  3818. C
  3819. C     ON OUTPUT
  3820. C
  3821. C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
  3822. C          IN ITS FIRST M COLUMNS.
  3823. C
  3824. C     NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS.
  3825. C
  3826. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3827. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3828. C
  3829. C     THIS VERSION DATED AUGUST 1983.
  3830. C
  3831. C     ------------------------------------------------------------------
  3832. C
  3833.       SUBROUTINE TRED1(NM,N,A,D,E,E2)
  3834. C
  3835.       INTEGER I,J,K,L,N,II,NM,JP1
  3836.       REAL A(NM,N),D(N),E(N),E2(N)
  3837.       REAL F,G,H,SCALE
  3838. C
  3839. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1,
  3840. C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
  3841. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  3842. C
  3843. C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX
  3844. C     TO A SYMMETRIC TRIDIAGONAL MATRIX USING
  3845. C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
  3846. C
  3847. C     ON INPUT
  3848. C
  3849. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  3850. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3851. C          DIMENSION STATEMENT.
  3852. C
  3853. C        N IS THE ORDER OF THE MATRIX.
  3854. C
  3855. C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
  3856. C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
  3857. C
  3858. C     ON OUTPUT
  3859. C
  3860. C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
  3861. C          FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER
  3862. C          TRIANGLE.  THE FULL UPPER TRIANGLE OF A IS UNALTERED.
  3863. C
  3864. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
  3865. C
  3866. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  3867. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
  3868. C
  3869. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  3870. C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
  3871. C
  3872. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3873. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3874. C
  3875. C     THIS VERSION DATED AUGUST 1983.
  3876. C
  3877. C     ------------------------------------------------------------------
  3878. C
  3879.       SUBROUTINE TRED2(NM,N,A,D,E,Z)
  3880. C
  3881.       INTEGER I,J,K,L,N,II,NM,JP1
  3882.       REAL A(NM,N),D(N),E(N),Z(NM,N)
  3883.       REAL F,G,H,HH,SCALE
  3884. C
  3885. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
  3886. C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
  3887. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  3888. C
  3889. C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A
  3890. C     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
  3891. C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
  3892. C
  3893. C     ON INPUT
  3894. C
  3895. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  3896. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  3897. C          DIMENSION STATEMENT.
  3898. C
  3899. C        N IS THE ORDER OF THE MATRIX.
  3900. C
  3901. C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
  3902. C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
  3903. C
  3904. C     ON OUTPUT
  3905. C
  3906. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
  3907. C
  3908. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  3909. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
  3910. C
  3911. C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
  3912. C          PRODUCED IN THE REDUCTION.
  3913. C
  3914. C        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED.
  3915. C
  3916. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3917. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3918. C
  3919. C     THIS VERSION DATED AUGUST 1983.
  3920. C
  3921. C     ------------------------------------------------------------------
  3922. C
  3923.       SUBROUTINE TRED3(N,NV,A,D,E,E2)
  3924. C
  3925.       INTEGER I,J,K,L,N,II,IZ,JK,NV,JM1
  3926.       REAL A(NV),D(N),E(N),E2(N)
  3927.       REAL F,G,H,HH,SCALE
  3928. C
  3929. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
  3930. C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
  3931. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  3932. C
  3933. C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS
  3934. C     A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
  3935. C     USING ORTHOGONAL SIMILARITY TRANSFORMATIONS.
  3936. C
  3937. C     ON INPUT
  3938. C
  3939. C        N IS THE ORDER OF THE MATRIX.
  3940. C
  3941. C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
  3942. C          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  3943. C
  3944. C        A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
  3945. C          INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
  3946. C          ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
  3947. C
  3948. C     ON OUTPUT
  3949. C
  3950. C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL
  3951. C          TRANSFORMATIONS USED IN THE REDUCTION.
  3952. C
  3953. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
  3954. C
  3955. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  3956. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
  3957. C
  3958. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  3959. C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
  3960. C
  3961. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  3962. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  3963. C
  3964. C     THIS VERSION DATED AUGUST 1983.
  3965. C
  3966. C     ------------------------------------------------------------------
  3967. C
  3968. C
  3969.       SUBROUTINE TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,M,W,IND,IERR,RV4,RV5)
  3970. C
  3971.       INTEGER I,J,K,L,M,N,P,Q,R,S,II,M1,M2,M11,M22,TAG,IERR,ISTURM
  3972.       REAL D(N),E(N),E2(N),W(M),RV4(N),RV5(N)
  3973.       REAL U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,TST1,TST2,EPSLON
  3974.       INTEGER IND(M)
  3975. C
  3976. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BISECT,
  3977. C     NUM. MATH. 9, 386-393(1967) BY BARTH, MARTIN, AND WILKINSON.
  3978. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 249-256(1971).
  3979. C
  3980. C     THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL
  3981. C     SYMMETRIC MATRIX BETWEEN SPECIFIED BOUNDARY INDICES,
  3982. C     USING BISECTION.
  3983. C
  3984. C     ON INPUT
  3985. C
  3986. C        N IS THE ORDER OF THE MATRIX.
  3987. C
  3988. C        EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED
  3989. C          EIGENVALUES.  IF THE INPUT EPS1 IS NON-POSITIVE,
  3990. C          IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE,
  3991. C          NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE
  3992. C          PRECISION AND THE 1-NORM OF THE SUBMATRIX.
  3993. C
  3994. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  3995. C
  3996. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  3997. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  3998. C
  3999. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  4000. C          E2(1) IS ARBITRARY.
  4001. C
  4002. C        M11 SPECIFIES THE LOWER BOUNDARY INDEX FOR THE DESIRED
  4003. C          EIGENVALUES.
  4004. C
  4005. C        M SPECIFIES THE NUMBER OF EIGENVALUES DESIRED.  THE UPPER
  4006. C          BOUNDARY INDEX M22 IS THEN OBTAINED AS M22=M11+M-1.
  4007. C
  4008. C     ON OUTPUT
  4009. C
  4010. C        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS
  4011. C          (LAST) DEFAULT VALUE.
  4012. C
  4013. C        D AND E ARE UNALTERED.
  4014. C
  4015. C        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
  4016. C          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
  4017. C          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
  4018. C          E2(1) IS ALSO SET TO ZERO.
  4019. C
  4020. C        LB AND UB DEFINE AN INTERVAL CONTAINING EXACTLY THE DESIRED
  4021. C          EIGENVALUES.
  4022. C
  4023. C        W CONTAINS, IN ITS FIRST M POSITIONS, THE EIGENVALUES
  4024. C          BETWEEN INDICES M11 AND M22 IN ASCENDING ORDER.
  4025. C
  4026. C        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
  4027. C          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
  4028. C          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
  4029. C          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC..
  4030. C
  4031. C        IERR IS SET TO
  4032. C          ZERO       FOR NORMAL RETURN,
  4033. C          3*N+1      IF MULTIPLE EIGENVALUES AT INDEX M11 MAKE
  4034. C                     UNIQUE SELECTION IMPOSSIBLE,
  4035. C          3*N+2      IF MULTIPLE EIGENVALUES AT INDEX M22 MAKE
  4036. C                     UNIQUE SELECTION IMPOSSIBLE.
  4037. C
  4038. C        RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS.
  4039. C
  4040. C     NOTE THAT SUBROUTINE TQL1, IMTQL1, OR TQLRAT IS GENERALLY FASTER
  4041. C     THAN TRIDIB, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND.
  4042. C
  4043. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  4044. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  4045. C
  4046. C     THIS VERSION DATED AUGUST 1983.
  4047. C
  4048. C     ------------------------------------------------------------------
  4049. C
  4050.       SUBROUTINE TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,Z,
  4051.      X                  IERR,RV1,RV2,RV3,RV4,RV5,RV6)
  4052. C
  4053.       INTEGER I,J,K,M,N,P,Q,R,S,II,IP,JJ,MM,M1,M2,NM,ITS,
  4054.      X        IERR,GROUP,ISTURM
  4055.       REAL D(N),E(N),E2(N),W(MM),Z(NM,MM),
  4056.      X       RV1(N),RV2(N),RV3(N),RV4(N),RV5(N),RV6(N)
  4057.       REAL U,V,LB,T1,T2,UB,UK,XU,X0,X1,EPS1,EPS2,EPS3,EPS4,
  4058.      X       NORM,TST1,TST2,EPSLON,PYTHAG
  4059. C
  4060. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRISTURM
  4061. C     BY PETERS AND WILKINSON.
  4062. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
  4063. C
  4064. C     THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL
  4065. C     SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL AND THEIR
  4066. C     ASSOCIATED EIGENVECTORS, USING BISECTION AND INVERSE ITERATION.
  4067. C
  4068. C     ON INPUT
  4069. C
  4070. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  4071. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  4072. C          DIMENSION STATEMENT.
  4073. C
  4074. C        N IS THE ORDER OF THE MATRIX.
  4075. C
  4076. C        EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED
  4077. C          EIGENVALUES.  IT SHOULD BE CHOSEN COMMENSURATE WITH
  4078. C          RELATIVE PERTURBATIONS IN THE MATRIX ELEMENTS OF THE
  4079. C          ORDER OF THE RELATIVE MACHINE PRECISION.  IF THE
  4080. C          INPUT EPS1 IS NON-POSITIVE, IT IS RESET FOR EACH
  4081. C          SUBMATRIX TO A DEFAULT VALUE, NAMELY, MINUS THE
  4082. C          PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE
  4083. C          1-NORM OF THE SUBMATRIX.
  4084. C
  4085. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  4086. C
  4087. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  4088. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  4089. C
  4090. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  4091. C          E2(1) IS ARBITRARY.
  4092. C
  4093. C        LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES.
  4094. C          IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND.
  4095. C
  4096. C        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF
  4097. C          EIGENVALUES IN THE INTERVAL.  WARNING. IF MORE THAN
  4098. C          MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL,
  4099. C          AN ERROR RETURN IS MADE WITH NO VALUES OR VECTORS FOUND.
  4100. C
  4101. C     ON OUTPUT
  4102. C
  4103. C        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS
  4104. C          (LAST) DEFAULT VALUE.
  4105. C
  4106. C        D AND E ARE UNALTERED.
  4107. C
  4108. C        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
  4109. C          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
  4110. C          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
  4111. C          E2(1) IS ALSO SET TO ZERO.
  4112. C
  4113. C        M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB).
  4114. C
  4115. C        W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER IF THE MATRIX
  4116. C          DOES NOT SPLIT.  IF THE MATRIX SPLITS, THE EIGENVALUES ARE
  4117. C          IN ASCENDING ORDER FOR EACH SUBMATRIX.  IF A VECTOR ERROR
  4118. C          EXIT IS MADE, W CONTAINS THOSE VALUES ALREADY FOUND.
  4119. C
  4120. C        Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
  4121. C          IF AN ERROR EXIT IS MADE, Z CONTAINS THOSE VECTORS
  4122. C          ALREADY FOUND.
  4123. C
  4124. C        IERR IS SET TO
  4125. C          ZERO       FOR NORMAL RETURN,
  4126. C          3*N+1      IF M EXCEEDS MM.
  4127. C          4*N+R      IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
  4128. C                     EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS.
  4129. C
  4130. C        RV1, RV2, RV3, RV4, RV5, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
  4131. C
  4132. C     THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM
  4133. C     APPEARS IN TSTURM IN-LINE.
  4134. C
  4135. C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
  4136. C
  4137. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
  4138. C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
  4139. C
  4140. C     THIS VERSION DATED AUGUST 1983.
  4141. C
  4142. C     ------------------------------------------------------------------
  4143. C
  4144.