1!
2C ******************************************************************
3C *  HERMDP.FOR  SOLVES EIGENVALUE PROBLEM                          *
4C *                                       H X = E X                 *
5C *  WHERE H IS A HERMITIAN MATRIX, X ARE THE EIGENVECTORS AND E    *
6C *  THE EIGENVALUES.                                               *
7C *                             OTTO F. SANKEY                      *
8C *                             DEPT. OF PHYSICS                    *
9C *                             ARIZONA STATE UNIVERSITY            *
10C *                             602 965-4334                        *
11C *                                                                 *
12C *             DOUBLE PRECISION   5/15/85                          *
13c *             updated            2/10/89 ofs.                     *
14C *******************************************************************
15C
16C TO USE:
17C               SUBROUTINE HERMDP(H,ZR,ZI,E,IDIM,M,W1,W2,icall)
18C H     = INPUT HAMILTONIAN.
19c         it contains real and imaginary parts of h as follows:
20c         h(i,j) = real (h(i,j))  , [ i.ge.j, put Real(lower triangle)
21c                                             in lower triangle].
22c         h(j,i) = imag (h(i,j))  , [ i.ge.j, put imag(lower triangle)
23c                                             in lower triangle].
24c
25C  IDIM = DIMENSION OF H EXACTLY AS IN CALLING ROUTINE.
26C     M = DIMENSION OF MATRIX YOU WANT DIAGONALIZED (M.LE.IDIM)
27C ZR,ZI = OUTPUT EIGENVECTORS (R MEANS REAL PART AND I IMAG. PART).
28c         zr(i,j) is the real part of the i'th component of the j'th
29c         eigenvector -- zr(component, eigenvector).
30C     E = OUTPUT EIGENVALUES
31C  icall= 3 for eigenvectors AND eigenvalues, icall.ne.3 for eigenvalues
32c                                             only.
33C W1, W2 ARE WORK VECTORS, W1(IDIM),W2(2,IDIM)
34c aid is a work vector dimensioned 501. for bigger matrices, please
35c                                       redimension.
36c
37c program computes trace & compares with sum of eigenvalues as a check.
38C
39C CALL SUBROUTINES: DPHTRI,DPIMTQ,DPHTBK
40c ==============================================================
41        SUBROUTINE hERMDP(H,ZR,ZI,E,IDIM,M,W1,W2,icall)
42c ==============================================================
43c        DOUBLE PRECISION H(IDIM,M),E(M),AID(501)
44        integer :: idim, m, icall
45        DOUBLE PRECISION H(IDIM,M),E(M),AID(5501),w3(5501)
46        DOUBLE PRECISION ZR(IDIM,M),ZI(IDIM,M),W1(M),W2(2,M)
47        DOUBLE PRECISION EIGSUM,TRACE
48c ==============================================================
49        COMMON /AIDIAG/AID
50c     ==============================================================
51        integer :: i, j, ier
52
53c        IF(M.GT.501)WRITE(*,*)' ERROR IN HERMDP ******* M>501'
54c        IF(M.GT.501)GO TO 999
55c      write(*,*)' welcome to hermdp.f ......'
56        IF(M.GT.5501)WRITE(*,*)' ERROR IN HERMDP ******* M>5501'
57        IF(M.GT.5501)GO TO 999
58C FIND TRACE OF MATRIX, to test sum of eigenvalues.
59       TRACE=0.D0
60       DO 10 I=1,M
6110     TRACE=TRACE+H(I,I)
62C INITIALIZE ZR,ZI. must set zr to identity.
63        DO 1 I=1,M
64        DO 2 J=1,M
652       ZR(I,J)=0.D0
661       ZR(I,I)=1.D0
67c ==============================================================
68c first call.
69        CALL DPHTRI(IDIM,M,H,E,W1,w3,W2)
70       IF(M.LT.300)GO TO 3
71       WRITE(*,*)' FIRST CALL FINISHED'
723      continue
73c ==============================================================
74c second call.
75       CALL DPIMTQ(IDIM,M,E,W1,ZR,IER)
76       IF(M.LT.300)GO TO 4
77       WRITE(*,*)' SECOND CALL FINISHED'
784      continue
79c ==============================================================
80c third call: for eigenvectors.
81       if(icall.eq.3)then
82           CALL DPHTBK(IDIM,M,H,W2,M,ZR,ZI)
83           IF(M.GT.300)WRITE(*,*)' THIRD CALL FINISHED'
84       else
85           write(*,*)' **caution** in hermdp, eigenvcs not desired'
86       end if
87       IF(M.GT.300.and.icall.eq.3)WRITE(*,*)' THIRD CALL FINISHED'
88c ==============================================================
89C SUM EIGENVALUES AND COMPARE WITH TRACE.
90       EIGSUM=0.D0
91       DO 11 I=1,M
9211     EIGSUM=EIGSUM+E(I)
93       IF(100.D0*DABS(EIGSUM-TRACE).LT.DABS(TRACE)*1.D-9)GOTO 12
94       WRITE(*,*)' ERROR: hERMDP SUM OF EIGENVALUES NOT EQUAL TO TRACE'
95       WRITE(*,*)' EIGSUM=',EIGSUM,' TRACE=',TRACE
96       IF(DABS(EIGSUM-TRACE).GT.DABS(TRACE)*1.D-08)
97     X  WRITE(*,*)' BAD ERROR IN hERMDP ******* !!!!!! ******'
9812      IF(IER.NE.0) WRITE(*,*)' ERROR IN HERMDP IER IN HERMDP=',IER
99c ==============================================================
100        RETURN
101999     STOP
102        END
103c ==============================================================
104C
105      SUBROUTINE DPHTRI(NM,N,AR,D,E,E2,TAU)
106C
107      INTEGER I,J,K,L,N,II,NM,JP1
108      DOUBLE PRECISION AR(NM,N),D(N),E(N),E2(N),TAU(2,N)
109c      DOUBLE PRECISION F,FI,G,GI,H,HH,SI,SCALE,AID(501)
110      DOUBLE PRECISION F,FI,G,GI,H,HH,SI,SCALE,AID(5501)
111      DOUBLE PRECISION DSQRT,CDABS,DABS
112      double precision ajunk
113       COMMON /AIDIAG/AID
114C      COMPLEX*16 DCMPLX
115C
116C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
117C     THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968)
118C     BY MARTIN, REINSCH, AND WILKINSON.
119C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
120C
121C     THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX
122C     TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING
123C     UNITARY SIMILARITY TRANSFORMATIONS.
124C
125C     ON INPUT-
126C
127C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
128C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
129C          DIMENSION STATEMENT,
130C
131C        N IS THE ORDER OF THE MATRIX,
132C
133C        AR AND AI CONTAIN THE REAL AND AIMAGINARY PARTS,
134C          RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX.
135C          ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
136C
137C     ON OUTPUT-
138C
139C       AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
140C          FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER
141C          TRIANGLES.  THEIR STRICT UPPER TRIANGLES AND THE
142C          DIAGONAL OF AR ARE UNALTERED,
143C
144C        D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX,
145C
146C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
147C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO,
148C
149C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
150C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED,
151C
152C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
153C
154C     ARITHMETIC IS REAL EXCEPT FOR THE USE OF THE SUBROUTINES
155C     CABS AND CMPLX IN COMPUTING COMPLEX ABSOLUTE VALUES.
156C
157C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
158C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
159C
160C     ------------------------------------------------------------------
161C
162      TAU(1,N) = 1.0D0
163      TAU(2,N) = 0.0D0
164c       DO 324 I=1,501
165       DO 324 I=1,5501
166324      AID(I)=0.D0
167C
168C
169      DO 100 I = 1, N
170  100 D(I) = AR(I,I)
171C
172C     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
173      DO  300 II = 1, N
174         I = N + 1 - II
175         L = I - 1
176         H = 0.0D0
177         SCALE = 0.0D0
178         IF (L .LT. 1) GO TO 130
179C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
180         DO 120 K = 1, L
181  120    SCALE = SCALE + DABS(AR(I,K)) + DABS(AR(K,I))
182C
183         IF (SCALE .NE. 0.0) GO TO 140
184         TAU(1,L) = 1.0D0
185         TAU(2,L) = 0.0D0
186  130    E(I) = 0.0D0
187         E2(I) = 0.0D0
188         GO TO 290
189C
190  140    DO 150 K = 1, L
191            AR(I,K) = AR(I,K) / SCALE
192            AR(K,I) = AR(K,I) / SCALE
193            H = H + AR(I,K) * AR(I,K) + AR(K,I) * AR(K,I)
194  150    CONTINUE
195C
196         E2(I) = SCALE * SCALE * H
197         G = DSQRT(H)
198         E(I) = SCALE * G
199         F=DSQRT(AR(I,L)**2+AR(L,I)**2)
200C         F = CDABS(DCMPLX(AR(I,L),AI(I,L)))
201C     ********** FORM NEXT DIAGONAL ELEMENT OF MATRIX T **********
202         IF (F .EQ. 0.0D0) GO TO 160
203         TAU(1,L) = (AR(L,I) * TAU(2,I) - AR(I,L) * TAU(1,I)) / F
204         SI = (AR(I,L) * TAU(2,I) + AR(L,I) * TAU(1,I)) / F
205         H = H + F * G
206         G = 1.0D0 + G / F
207         AR(I,L) = G * AR(I,L)
208         AR(L,I) = G * AR(L,I)
209         IF (L .EQ. 1) GO TO 270
210         GO TO 170
211  160    TAU(1,L) = -TAU(1,I)
212         SI = TAU(2,I)
213         AR(I,L) = G
214  170    F = 0.0D0
215C
216         DO 240 J = 1, L
217            G = 0.0D0
218            GI = 0.0D0
219C     ********** FORM ELEMENT OF A*U **********
220            DO 180 K = 1, J
221       AJUNK=AR(K,J)
222       IF(K.EQ.J)AJUNK=AID(K)
223               G = G + AR(J,K) * AR(I,K) + AJUNK * AR(K,I)
224               GI = GI - AR(J,K) * AR(K,I) + AJUNK * AR(I,K)
225  180       CONTINUE
226C
227            JP1 = J + 1
228            IF (L .LT. JP1) GO TO 220
229C
230            DO 200 K = JP1, L
231               G = G + AR(K,J) * AR(I,K) - AR(J,K) * AR(K,I)
232               GI = GI - AR(K,J) * AR(K,I)- AR(J,K) * AR(I,K)
233  200       CONTINUE
234C     ********** FORM ELEMENT OF P **********
235  220       E(J) = G / H
236            TAU(2,J) = GI / H
237            F = F + E(J) * AR(I,J) - TAU(2,J) * AR(J,I)
238  240    CONTINUE
239C
240         HH = F / (H + H)
241C     ********** FORM REDUCED A **********
242         DO 260 J = 1, L
243            F = AR(I,J)
244            G = E(J) - HH * F
245            E(J) = G
246            FI = -AR(J,I)
247            GI = TAU(2,J) - HH * FI
248            TAU(2,J) = -GI
249C
250            DO 260 K = 1, J
251       AJUNK=AR(K,J)
252       IF(J.EQ.K)AJUNK=AID(K)
253               AR(J,K) = AR(J,K) - F * E(K) - G * AR(I,K)
254     X                           + FI * TAU(2,K) + GI * AR(K,I)
255       IF(K.EQ.J)GO TO 17
256               AR(K,J) = AR(K,J) - F * TAU(2,K) - G * AR(K,I)
257     X                           - FI * E(K) - GI * AR(I,K)
258       GO TO 260
25917               AID(K)=AID(K) - F * TAU(2,K) - G * AR(K,I)
260     X                           - FI * E(K) - GI * AR(I,K)
261  260    CONTINUE
262C
263  270    DO 280 K = 1, L
264            AR(I,K) = SCALE * AR(I,K)
265            AR(K,I) = SCALE * AR(K,I)
266  280    CONTINUE
267C
268         TAU(2,L) = -SI
269  290    HH = D(I)
270         D(I) = AR(I,I)
271         AR(I,I) = HH
272         AID(I) = SCALE*DSQRT(H)
273  300 CONTINUE
274C
275      RETURN
276C     ********** LAST CARD OF DPHTRI **********
277      END
278c ==============================================================
279      SUBROUTINE DPIMTQ(NM,N,D,E,Z,IERR)
280C
281      INTEGER I,J,K,L,M,N,II,NM,MML,IERR
282      DOUBLE PRECISION D(N),E(N),Z(NM,N)
283      DOUBLE PRECISION B,C,F,G,P,R,S,MACHEP
284      DOUBLE PRECISION DSQRT,DABS,DSIGN
285C
286C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE DPIMTQ,
287C     NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON,
288C     AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
289C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
290C
291C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
292C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD.
293C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
294C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
295C     FULL MATRIX TO TRIDIAGONAL FORM.
296C
297C     ON INPUT-
298C
299C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
300C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
301C          DIMENSION STATEMENT,
302C
303C        N IS THE ORDER OF THE MATRIX,
304C
305C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
306C
307C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
308C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
309C
310C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
311C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
312C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
313C          THE IDENTITY MATRIX.
314C
315C      ON OUTPUT-
316C
317C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
318C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
319C          UNORDERED FOR INDICES 1,2,...,IERR-1,
320C
321C        E HAS BEEN DESTROYED,
322C
323C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
324C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
325C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
326C          EIGENVALUES,
327C
328C        IERR IS SET TO
329C          ZERO       FOR NORMAL RETURN,
330C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
331C                     DETERMINED AFTER 30 ITERATIONS.
332C
333C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
334C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
335C
336C     ------------------------------------------------------------------
337C
338C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
339C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
340C     MACHEP=2**(-26) SP.    2**(-56) DP. (26 AND 56 ARE MAXIMUMS!)
341C     WE USE 54 FOR SAFETY. (1.D0+MACHEP.GT.1.D0) DEFINES MINIMUM MACHEP.
342C                **********
343      MACHEP=2.D0**(-54)
344C
345      IERR = 0
346      IF (N .EQ. 1) GO TO 1001
347C
348      DO 100 I = 2, N
349  100 E(I-1) = E(I)
350C
351      E(N) = 0.0D0
352C
353      DO 240 L = 1, N
354         J = 0
355C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
356  105    DO 110 M = L, N
357            IF (M .EQ. N) GO TO 120
358            IF (DABS(E(M)) .LE. MACHEP * (DABS(D(M)) + DABS(D(M+1))))
359     X         GO TO 120
360  110    CONTINUE
361C
362  120    P = D(L)
363         IF (M .EQ. L) GO TO 240
364         IF (J .EQ. 30) GO TO 1000
365         J = J + 1
366C     ********** FORM SHIFT **********
367         G = (D(L+1) - P) / (2.0D0 * E(L))
368         R = DSQRT(G*G+1.0D0)
369         G = D(M) - P + E(L) / (G + DSIGN(R,G))
370         S = 1.0D0
371         C = 1.0D0
372         P = 0.0D0
373         MML = M - L
374C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
375         DO 200 II = 1, MML
376            I = M - II
377            F = S * E(I)
378            B = C * E(I)
379            IF (DABS(F)  .LT. DABS(G))  GO TO 150
380            C = G / F
381            R =DSQRT(C*C+1.0D0)
382            E(I+1) = F * R
383            S = 1.0D0 / R
384            C = C * S
385            GO TO 160
386  150       S = F / G
387            R = DSQRT(S*S+1.0D0)
388            E(I+1) = G * R
389            C = 1.0D0 / R
390            S = S * C
391  160       G = D(I+1) - P
392            R = (D(I) - G) * S + 2.0D0 * C * B
393            P = S * R
394            D(I+1) = G + P
395            G = C * R - B
396C     ********** FORM VECTOR **********
397            DO 180 K = 1, N
398               F = Z(K,I+1)
399               Z(K,I+1) = S * Z(K,I) + C * F
400               Z(K,I) = C * Z(K,I) - S * F
401  180       CONTINUE
402C
403  200    CONTINUE
404C
405         D(L) = D(L) - P
406         E(L) = G
407         E(M) = 0.0D0
408         GO TO 105
409  240 CONTINUE
410C     ********** ORDER EIGENVALUES AND EIGENVECTORS **********
411      DO 300 II = 2, N
412         I = II - 1
413         K = I
414         P = D(I)
415C
416         DO 260 J = II, N
417            IF (D(J) .GE. P) GO TO 260
418            K = J
419            P = D(J)
420  260    CONTINUE
421C
422         IF (K .EQ. I) GO TO 300
423         D(K) = D(I)
424         D(I) = P
425C
426         DO 280 J = 1, N
427            P = Z(J,I)
428            Z(J,I) = Z(J,K)
429            Z(J,K) = P
430  280    CONTINUE
431C
432  300 CONTINUE
433C
434      GO TO 1001
435C     ********** SET ERROR -- NO CONVERGENCE TO AN
436C               EIGENVALUE AFTER 30 ITERATIONS **********
437 1000 IERR = L
438 1001 RETURN
439C     ********** LAST CARD OF DPIMTQ **********
440      END
441c ==============================================================
442      SUBROUTINE DPHTBK(NM,N,AR,TAU,M,ZR,ZI)
443C
444      INTEGER I,J,K,L,M,N,NM
445      DOUBLE PRECISION AR(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M)
446c      DOUBLE PRECISION H,S,SI,AID(501)
447      DOUBLE PRECISION H,S,SI,AID(5501)
448       COMMON /AIDIAG/AID
449C
450C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
451C     THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968)
452C     BY MARTIN, REINSCH, AND WILKINSON.
453C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
454C
455C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN
456C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
457C     REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  DPHTRI.
458C
459C     ON INPUT-
460C
461C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
462C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
463C          DIMENSION STATEMENT,
464C
465C        N IS THE ORDER OF THE MATRIX,
466C
467C        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
468C          FORMATIONS USED IN THE REDUCTION BY  DPHTRI  IN THEIR
469C          FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR.
470C
471C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS,
472C
473C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
474C
475C        ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
476C          IN ITS FIRST M COLUMNS.
477C
478C     ON OUTPUT-
479C
480C        ZR AND ZI CONTAIN THE REAL AND AIMAGINARY PARTS,
481C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
482C          IN THEIR FIRST M COLUMNS.
483C
484C     NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR
485C     IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED.
486C
487C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
488C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
489C
490C     ------------------------------------------------------------------
491        IF(M .EQ. 0 ) GOTO 200
492C     ********** TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC
493C                TRIDIAGONAL MATRIX TO THOSE OF THE HERMITIAN
494C                TRIDIAGONAL MATRIX. **********
495      DO 50 K = 1, N
496C
497         DO 50 J = 1, M
498            ZI(K,J) = - ZR(K,J) * TAU(2,K)
499            ZR(K,J) = ZR(K,J) * TAU(1,K)
500   50 CONTINUE
501C
502      IF (N .EQ. 1) GO TO 200
503C     ********** RECOVER AND APPLY THE HOUSEHOLDER MATRICES **********
504      DO 140 I = 2, N
505         L = I - 1
506         H = AID(I)
507         IF (H .EQ. 0.0D0) GO TO 140
508C
509         DO 130 J = 1, M
510            S = 0.0D0
511            SI = 0.0D0
512C
513            DO 110 K = 1, L
514               S = S + AR(I,K) * ZR(K,J) - AR(K,I) * ZI(K,J)
515               SI = SI + AR(I,K) * ZI(K,J) + AR(K,I) * ZR(K,J)
516  110       CONTINUE
517C
518        S=(S/H)/H
519        SI=(SI/H)/H
520C
521            DO 120 K = 1, L
522               ZR(K,J) = ZR(K,J) - S * AR(I,K) - SI * AR(K,I)
523               ZI(K,J) = ZI(K,J) - SI * AR(I,K) + S * AR(K,I)
524  120       CONTINUE
525C
526  130    CONTINUE
527C
528  140 CONTINUE
529C
530  200 RETURN
531C     ********** LAST CARD OF DPHTBK **********
532      END
533