1CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2C     this subroutine is moved to the top of the
3C     source file in order to resolve the mis-
4C     matching parameter types later in code.
5C     2005-10-11 Tommi Hassinen
6CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
7      SUBROUTINE EC08C(A,B,VALUE,VEC,N,IV,W)
8C
9C TO FIND THE EIGENVALUES AND VECTORS OF A TRI-DIAGONAL
10C  HERMITIAN MATRIX.
11      REAL VALUE(*),W(*),PCK(2),ONE,ZERO,VEC(*)
12      COMPLEX A(*),B(*),DN,UPCK
13      EQUIVALENCE (PCK(1),UPCK)
14C  WE TREAT VEC AS IF IT IS DEFINED AS COMPLEX VEC(IV,N)
15C  IN THE CALLING PROGRAM.
16      DATA ONE, ZERO/1.0,0.0/
17      IV2=IV+IV
18      N2=N+N
19      IL=IV2*(N-1)+1
20      W(1)=A(1)
21C
22C  THE HERMITIAN FORM IS TRANSFORMED INTO A REAL FORM
23      IF(N-1)80,80,10
24   10 DO 20 I=2,N
25         W(I)=A(I)
26   20 W(I+N)=CABS(B(I))
27C
28C  FIND THE EIGENVALUES AND VECTORS OF THE REAL FORM
29   30 CALL EA08C(W,W(N+1),VALUE,VEC,N,IV2,W(N2+1))
30C
31C THE VECTORS IN VEC AT THIS POINT ARE REAL,WE NOW EXPAND THEM
32C INTO VEC AS THOUGH THEY WERE COMPLEX.
33      DO 50 I=1,IL,IV2
34         DO 40 J=1,N
35            K=N-J
36            L=K+K
37   40    VEC(I+L)=VEC(I+K)
38   50 VEC(I+1)=ZERO
39      IF(N.LE.1)GO TO 80
40      DN=ONE
41      K=1
42C
43C TRANSFORM VECTORS OF REAL FORM TO THOSE OF COMPLEX FORM.
44      DO 70 I=4,N2,2
45         K=K+1
46         UPCK=ONE
47         IF(W(K+N).NE.ZERO)UPCK=DN*CONJG(B(K))/W(K+N)
48         I1=IL-1+I
49         DO 60 J=I,I1,IV2
50            VEC(J)=VEC(J-1)*PCK(2)
51   60    VEC(J-1)=VEC(J-1)*PCK(1)
52   70 DN=UPCK
53   80 RETURN
54      END
55CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
56C     this subroutine is moved to the top of the
57C     source file in order to resolve the mis-
58C     matching parameter types later in code.
59C     2005-10-11 Tommi Hassinen
60CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
61      SUBROUTINE ME08B (A,Q,B,N,IA)
62      REAL A(IA,*),Q(2,*),B(IA,*)
63      DO 10 I=1,N
64         A(1,I)=A(1,I) -Q(1,1)*B(1,I)+Q(2,1)*B(2,I)
65     1 -Q(1,I)*B(1,1)+Q(2,I)*B(2,1)
66   10 A(2,I)=A(2,I)-Q(2,1)*B(1,I)-Q(1,1)*B(2,I)
67     1 +Q(2,I)*B(1,1)+Q(1,I)*B(2,1)
68      RETURN
69      END
70      SUBROUTINE CDIAG(A,VALUE,VEC,N, NEED)
71C
72C TO FIND THE EIGENVALUES AND EIGENVECTORS OF A HERMITIAN MATRIX.
73      REAL VALUE(*),H
74      COMPLEX W(6000)
75      COMPLEX A(N,*),VEC(N,*)
76      COMPLEX  FM06AS
77      IA=N
78      IV=N
79C
80C REDUCE MATRIX TO A TRI-DIAGONAL HERMITIAN MATRIX.
81      CALL ME08A(A,W,W(N+1),N,IA,W(2*N+1))
82C
83C FIND THE EIGENVALUES AND EIGENVECTORS OF THE TRI-DIAGONAL MATRIX
84      CALL EC08C(W,W(N+1),VALUE,VEC,N,IV,W(2*N+1))
85      IF(NEED.EQ.0) GOTO 50
86      IF(N.LT.2)RETURN
87C
88C THE EIGENVECTORS OF THE ORIGINAL MATRIX ARE NOW FOUND BY
89C BACK TRANSFORMATION USING INFORMATION STORE IN THE UPPER
90C TRIANGLE OF MATRIX A (BY ME08)
91      DO 40 II=3,N
92         I=N-II+1
93         H=W(N+I+1)*CONJG(A(I,I+1))
94         IF(H)10,40,10
95   10    DO30L=1,N
96            I1=I+1
97            S=FM06AS(N-I,A(I,I+1),IA,VEC(I+1,L),1)
98            S=S/H
99            DO20K=I1,N
100   20       VEC(K,L)=VEC(K,L)+CONJG(A(I,K))*S
101   30    CONTINUE
102   40 CONTINUE
103   50 CALL SORT(VALUE,VEC,N)
104      RETURN
105      END
106      SUBROUTINE EA08C(A,B,VALUE,VEC,M,IV,W)
107C  STANDARD FORTRAN 66 (A VERIFIED PFORT SUBROUTINE)
108      DIMENSION A(*),B(*),VALUE(*),VEC(*),W(*)
109      DATA EPS/1.E-6/,A34/0.0/
110C     THIS USES QR ITERATION TO FIND THE EIGENVALUES AND EIGENVECTORS
111C  OF THE SYMMETRIC TRIDIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS ARE
112C  A(I),I=1,M AND OFF-DIAGONAL ELEMENTS ARE B(I),I=2,M.  THE ARRAY
113C  W IS USED FOR WORKSPACE AND MUST HAVE DIMENSION AT LEAST 2*M.
114C  WE TREAT VEC AS IF IT HAD DIMENSIONS (IV,M).
115      SML=EPS*FLOAT(M)
116      CALL EA09C(A,B,W(M+1),M,W)
117C     SET VEC TO THE IDENTITY MATRIX.
118      DO 20 I=1,M
119         VALUE(I)=A(I)
120         W(I)=B(I)
121         K=(I-1)*IV+1
122         L=K+M-1
123         DO 10 J=K,L
124   10    VEC(J)=0.
125         KI=K+I-1
126   20 VEC(KI)=1.
127      K=0
128      IF(M.EQ.1)RETURN
129      DO 100 N3=2,M
130         N2=M+2-N3
131C     EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2
132         MN2=M+N2
133         ROOT=W(MN2)
134         DO 80 ITER=1,20
135            BB=(VALUE(N2)-VALUE(N2-1))*0.5
136            CC=W(N2)*W(N2)
137            A22=VALUE(N2)
138            IF(CC.NE.0.0)A22=A22+CC/(BB+SIGN(1.0,BB)*SQRT(BB*BB+CC))
139            DO 30 I=1,N2
140               MI=M+I
141               IF(ABS(ROOT-A22).LE.ABS(W(MI)-A22))GO TO 30
142               ROOT=W(MI)
143               MN=M+N2
144               W(MI)=W(MN)
145               W(MN)=ROOT
146   30       CONTINUE
147            DO 40 II=2,N2
148               N1=2+N2-II
149               IF(ABS(W(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*SML)GO
150     1 TO 50
151   40       CONTINUE
152            N1=1
153   50       IF(N2.EQ.N1) GO TO 100
154            N2M1=N2-1
155            IF(ITER.GE.3)ROOT=A22
156            K=K+1
157            A22=VALUE(N1)
158            A12=A22-ROOT
159            A23=W(N1+1)
160            A13=A23
161            DO70 I=N1,N2M1
162               A33=VALUE(I+1)
163               IF(I.NE.N2M1)A34=W  (I+2)
164               S=SIGN(SQRT(A12*A12+A13*A13),A12)
165               SI=A13/S
166               CO=A12/S
167               JK=I*IV+1
168               J1=JK-IV
169               J2=J1+MIN0(M,I+K)-1
170               DO 60 JI=J1,J2
171                  V1=VEC(JI)
172                  V2=VEC(JK)
173                  VEC(JI)=V1*CO+V2*SI
174                  VEC(JK)=V2*CO-V1*SI
175   60          JK=JK+1
176               IF(I.NE.N1)  W(I)=S
177               A11=CO*A22+SI*A23
178               A12=CO*A23+SI*A33
179               A13=SI*A34
180               A21=CO*A23-SI*A22
181               A22=CO*A33-SI*A23
182               A23=CO*A34
183               VALUE(I)=A11*CO+A12*SI
184               A12=-A11*SI+A12*CO
185               W(I+1)=A12
186   70       A22=A22*CO-A21*SI
187   80    VALUE(N2)=A22
188         WRITE(6,90)
189   90    FORMAT(48H1CYCLE DETECTED IN SUBROUTINE EA08 -STOPPING NOW)
190         STOP
191  100 CONTINUE
192C     RAYLEIGH QUOTIENT
193      DO 120 J=1,M
194         K=(J-1)*IV
195         XX=VEC(K+1)**2
196         XAX=XX*A(1)
197         DO 110 I=2,M
198            KI=K+I
199            XX=XX+VEC(KI)**2
200  110    XAX=XAX+VEC(KI)*(2.*B(I)*VEC(KI-1)+A(I)*VEC(KI))
201  120 VALUE(J)=XAX/XX
202      RETURN
203      END
204      SUBROUTINE EA09C(A,B,VALUE,M,OFF)
205C  STANDARD FORTRAN 66 (A VERFIED PFORT SUBROUTINE)
206      DIMENSION A(*),B(*),VALUE(*),OFF(*)
207      DATA A34/0.0/,EPS/1.0E-6/
208      SML=EPS*FLOAT(M)
209      VALUE(1)=A(1)
210      IF(M.EQ.1)RETURN
211      DO 10 I=2,M
212         VALUE(I)=A(I)
213   10 OFF(I)=B(I)
214C     EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2
215      MAXIT=10*M
216      DO 80 ITER=1,MAXIT
217         DO 40 N3=2,M
218            N2=M+2-N3
219            DO 20 II=2,N2
220               N1=2+N2-II
221               IF(ABS(OFF(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*SML)
222     1GO TO 30
223   20       CONTINUE
224            N1=1
225   30       IF(N2.NE.N1) GO TO 50
226   40    CONTINUE
227         RETURN
228C     ROOT  IS THE EIGENVALUE OF THE BOTTOM 2*2 MATRIX THAT IS NEAREST
229C     TO THE LAST MATRIX ELEMENT AND IS USED TO ACCELERATE THE
230C     CONVERGENCE
231   50    BB=(VALUE(N2)-VALUE(N2-1))*0.5
232         CC=OFF(N2)*OFF(N2)
233         SBB=1.
234         IF(BB.LT.0.)SBB=-1.
235         ROOT=VALUE(N2)+CC/(BB+SBB*SQRT(BB*BB+CC))
236         N2M1=N2-1
237   60    A22=VALUE(N1)
238         A12=A22-ROOT
239         A23=OFF(N1+1)
240         A13=A23
241         DO 70 I=N1,N2M1
242            A33=VALUE(I+1)
243            IF(I.NE.N2M1)A34=OFF(I+2)
244            S=SQRT(A12*A12+A13*A13)
245            SI=A13/S
246            CO=A12/S
247            IF(I.NE.N1)OFF(I)=S
248            A11=CO*A22+SI*A23
249            A12=CO*A23+SI*A33
250            A13=SI*A34
251            A21=CO*A23-SI*A22
252            A22=CO*A33-SI*A23
253            A23=CO*A34
254            VALUE(I)=A11*CO+A12*SI
255            A12=-A11*SI+A12*CO
256            OFF(I+1)=A12
257   70    A22=A22*CO-A21*SI
258   80 VALUE(N2)=A22
259      WRITE(6,90)
260   90 FORMAT(39H1LOOPING DETECTED IN EA09-STOPPING NOW )
261      STOP
262      END
263      COMPLEX FUNCTION FM06AS(N,A,IA,B,IB)
264      IMPLICIT COMPLEX (A-H,O-Z)
265      COMPLEX A(*), B(*)
266*******************************************************
267*
268*    FM06AS - A FUNCTION ROUTINE TO COMPUTE THE VALUE OF THE
269*      INNER PRODUCT, OR DOT PRODUCT, OF TWO SINGLE PRECISION
270*      COMPLEX VECTORS, ACCUMULATING THE INTERMEDIATE PRODUCTS
271*      DOUBLE PRECISION.  THE ELEMENTS OF EACH VECTOR CAN BE
272*      STORED IN ANY FIXED DISPLACEMENT FROM NEIGHBOURING
273*      ELEMENTS.
274*
275*    COMPUTES: SUM(J=1,N) A((J-1)*IA+1)*B((J-1)*IB+1)
276*
277*          W = FM06AS(N,A,IA,B,IB)
278*
279*    N   INTEGER SCALAR; (USER:*); LENGTH OF THE VECTORS A AND B.
280*        IF N <= 0 THE INNER PRODUCT VALUE IS DEFINED TO BE ZERO.
281*    A   COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY
282*        CONTAINING THE 1ST VECTOR.  THE FORTRAN CONVENTION OF STORING
283*        REAL AND IMAGINARY PARTS IN ADJACENT WORDS IS ASSUMED.
284*    IA  INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF
285*        AN ELEMENT IN THE ARRAY A TO ITS NEIGHBOUR, I.E. THE VECTOR
286*        ELEMENTS ARE IN A(1), A(IA+1), A(2*IA+1),...
287*        IF IA < 0 THE ELEMENTS ARE ASSUMED TO BE STORED IN
288*        A(1-(N-1)*IA), A(1-(N-2)*IA),..., A(1-IA), A(1).
289*    B   COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY
290*        CONTAINING THE SECOND VECTOR.  TREAT LIKE A.
291*    IB  INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF
292*        AN ELEMENT IN B TO ITS NEIGHBOUR. TREAT LIKE IA.
293*    FM06AS  COMPLEX FUNCTION; (*:FM06AS); THE INNER PRODUCT VALUE.
294*        IT IS RETURNED DOUBLE PRECISION, THE REAL PART IN FLT PNT
295*        REG 0 AND THE IMAGINARY PART IN FLT PNT REG 2.
296*
297*    THIS ROUTINE IS WRITTEN IN FORTRAN.
298*
299*--------------------------------------------------------*
300      SUM=(0.0,0.0)
301      DO 10 I=1,N
302   10 SUM=SUM+A((I-1)*IA+1)*B((I-1)*IB+1)
303      FM06AS=SUM
304      RETURN
305      END
306      COMPLEX FUNCTION FM06BS(N,A,IA,B,IB)
307      IMPLICIT COMPLEX (A-H,O-Z)
308      COMPLEX A(*), B(*)
309*******************************************************
310*
311*    FM06BS - A FUNCTION ROUTINE TO COMPUTE THE VALUE OF THE
312*      INNER PRODUCT, OR DOT PRODUCT, OF A SIGLE    PRECISION
313*      COMPLEX VECTORS, ACCUMULATING THE INTERMEDIATE PRODUCTS
314*      DOUBLE PRECISION.  THE ELEMENTS OF EACH VECTOR CAN BE
315*      STORED IN ANY FIXED DISPLACEMENT FROM NEIGHBOURING
316*      ELEMENTS.
317*
318*    COMPUTES: SUM(J=1,N) A((J-1)*IA+1)*B((J-1)*IB+1)
319*
320*          W = FM06BS(N,A,IA,B,IB)
321*
322*    N   INTEGER SCALAR; (USER:*); LENGTH OF THE VECTORS A AND B.
323*        IF N <= 0 THE INNER PRODUCT VALUE IS DEFINED TO BE ZERO.
324*    A   COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY
325*        CONTAINING THE 1ST VECTOR.  THE FORTRAN CONVENTION OF STORING
326*        REAL AND IMAGINARY PARTS IN ADJACENT WORDS IS ASSUMED.
327*    IA  INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF
328*        AN ELEMENT IN THE ARRAY A TO ITS NEIGHBOUR, I.E. THE VECTOR
329*        ELEMENTS ARE IN A(1), A(IA+1), A(2*IA+1),...
330*        IF IA < 0 THE ELEMENTS ARE ASSUMED TO BE STORED IN
331*        A(1-(N-1)*IA), A(1-(N-2)*IA),..., A(1-IA), A(1).
332*    B   COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY
333*        CONTAINING THE SECOND VECTOR.  TREAT LIKE A.
334*    IB  INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF
335*        AN ELEMENT IN B TO ITS NEIGHBOUR. TREAT LIKE IA.
336*    FM06BS  COMPLEX FUNCTION; (*:FM06BS); THE INNER PRODUCT VALUE.
337*        IT IS RETURNED DOUBLE PRECISION, THE REAL PART IN FLT PNT
338*        REG 0 AND THE IMAGINARY PART IN FLT PNT REG 2.
339*
340*    THIS ROUTINE IS WRITTEN IN FORTRAN.
341*
342*--------------------------------------------------------*
343      SUM=(0.0,0.0)
344      DO 10 I=1,N
345   10 SUM=SUM+A((I-1)*IA+1)*CONJG(B((I-1)*IB+1))
346      FM06BS=SUM
347      RETURN
348      END
349      SUBROUTINE ME08A(A,ALPHA,BETA,N,IA,Q)
350      COMPLEX  A(IA,*),ALPHA(*),BETA(*),Q(*),CW,QJ
351      COMPLEX  FM06AS,FM06BS
352      REAL PP,ZERO,P5,H
353      REAL S1,PP1
354      DATA ZERO/0.0/, P5 /0.50/
355C**************************************************************
356C THE REDUCTION OF FULL HERMITIAN MATRIX INTO TRI-DIAGONAL HERMITIAN
357C FORM IS DONE IN N-2 STEPS.AT THE I TH STEP ZEROS ARE INTRODUCED IN
358C THE I TH ROW AND COLUMNS WITHOUT DESTROYING PREVIOUSLY PRODUCED ZEROS
359C
360C                                                    H
361C AT THE I TH STEP WE HAVE  A =P A   P  WITH P =I-U U  /K
362C                            I  I I-1 I       I    I I   I
363C
364C WHERE U =(0,0,...,A     (1+S /T ),A     ,...,A   )
365C        I           I,I+1    I  I   I,I+2      I,N
366C        2   N          2        2                2
367C       S = SUM  ] A   ]     K =S +T   AND  T = SQRT(]A     ] S )
368C        I   J=I+1  I,J       I  I  I        I         I,I+1   I
369C
370C COMPUTATIONAL DETAILS AT THE I TH STAGE ARE (1) FORM S  ,K   THEN
371C                                                       I   I
372C
373C                                     H                      H     H
374C (2) Q =A   U  /K    (3) Q =Q -.5U (U Q /K ) (4) A =A   -U Q  -Q U
375C      I  I-1 I   I        I  I    I  I I  I       I  I-1  I I   I I
376C
377C  THE VECTORS U BEING APPROPIATELY IN A.
378      IF(N.LE.0)GO TO 90
379      DO 10 J=1,N
380         ALPHA(J)=A(J,J)
381         DO 10 I=1,J
382   10 A(I,J)=CONJG(A(J,I))
383      IF(N-2)90,80,20
384   20 N2=N-2
385      DO 60 I=1,N2
386         I1=I+1
387C                       (1)
388         CW=FM06BS(N-I,A(I,I+1),IA,A(I,I+1),IA)
389         PP=CW
390         PP1=SQRT(PP)
391         BETA(I+1)=CMPLX(-PP1,ZERO)
392         S1=CABS(A(I,I+1))
393         IF(S1.GT.ZERO)BETA(I+1)=BETA(I+1)*A(I,I+1)/S1
394         IF(PP.LE.1.D-15)GO TO 60
395         H=PP+PP1*S1
396         A(I,I+1)=A(I,I+1)-BETA(I+1)
397C                       (2)
398         DO 30 K=I1,N
399            QJ=FM06AS(-(I-K),A(I+1,K),1,A(I,I+1),IA)
400            QJ=FM06BS(N-K,A(K,K+1),IA,A(I,K+1),IA)+CONJG(QJ)
401   30    Q(K)=QJ/H
402C                       (3)
403         CW=FM06AS(N-I,A(I,I+1),IA,Q(I+1),1)
404         PP=CW*P5/H
405         DO 40 K=I1,N
406   40    Q(K)=Q(K)-PP*CONJG(A(I,K))
407C                       (4)
408         DO 50 K=I1,N
409   50    CALL ME08B (A(K,K),Q(K),A(I,K),N-K+1,IA*2)
410   60 CONTINUE
411      DO 70 I=2,N
412         QJ=ALPHA(I)
413         ALPHA(I)=A(I,I)
414   70 A(I,I)=QJ
415   80 BETA(N)=A(N-1,N)
416   90 RETURN
417      END
418      SUBROUTINE SORT(VAL,VEC,N)
419      COMPLEX VEC(N,*), SUM
420      REAL    VAL(*)
421      DO 30 I=1,N
422         X=1.E9
423         DO 10 J=I,N
424            IF(VAL(J).LT.X) THEN
425               K=J
426               X=VAL(J)
427            ENDIF
428   10    CONTINUE
429         DO 20 J=1,N
430            SUM=VEC(J,K)
431            VEC(J,K)=VEC(J,I)
432   20    VEC(J,I)=SUM
433         VAL(K)=VAL(I)
434         VAL(I)=X
435   30 CONTINUE
436      RETURN
437      END
438