← Back to branch summary

~nickpapior/siesta/trunk-kpoint-dos

« back to all changes in this revision

Viewing changes to Src/cdiag.f

  • Committer: Alberto Garcia
  • Date: 2004-11-25 18:49:43 UTC
  • Revision ID: Arch-1:siesta@uam.es--2004%siesta-devel--reference--0.11--patch-1
Siesta 0.11 -- imported from CVS
Import from cvs using date instead of siesta-0-11-release tag, since
the Pseudo structure was not properly integrated at that time and
did not get the tag.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/cdiag.f
1
 
      SUBROUTINE CDIAG (H,NH,S,NS,N,E,Z,NZ,AUX)
2
 
 
3
 
*  SOLVES THE COMPLEX GENERAL EIGENVALUE PROBLEM   H*Z = E*S*Z
4
 
*  WHERE S AND H ARE COMPLEX HERMITIAN MATRICES
5
 
 
6
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
7
 
      DIMENSION H(2*NH*N),S(2*NS*N),Z(2*NZ*N),E(N),AUX(N,5)
8
 
 
9
 
*  START TIME COUNT
10
 
      CALL TIMER('CDIAG',1)
11
 
 
12
 
*  REDUCE GENERAL EIGENVALUE PROBLEM TO A NORMAL ONE
13
 
      CALL REDUCC (NH,NS,N,H,S,AUX)
14
 
 
15
 
*  NEXT LINE FOR IMSL
16
 
C     CALL EIGCH (H,NH,11,E,Z,NZ,AUX(1,2),IER)
17
 
 
18
 
*  NEXT LINE FOR SLATEC. AUX 4*N REALS
19
 
C     CALL CHIEV (H,NH,N,E,Z,NZ,AUX(1,2),1,INFO)
20
 
 
21
 
*  NEXT THREE LINES FOR IBM-ESSL. AUX 4*N REALS
22
 
C     CALL PACK (NH,N,H,H,+1)
23
 
C     CALL ZHPEV (1,H,E,Z,NZ,N,AUX(1,2),4*N)
24
 
C     CALL PACK (NH,N,H,H,-1)
25
 
 
26
 
*  NEXT LINES FOR EISPACK. AUX 2*N REALS.
27
 
      IF (NH.NE.NZ) THEN
28
 
        WRITE (6,*) 'CDIAG: NH AND NZ MUST BE EQUAL. NH,NZ =',NH,NZ
29
 
        STOP 2
30
 
      ENDIF
31
 
      CALL TRANS (2,NH*N,H,H,Z)
32
 
      CALL EISCH1 (NH,N,H,H(NH*N+1),E,Z,Z(NZ*N+1),IER,AUX(1,2))
33
 
      IF (IER.NE.0) THEN
34
 
        WRITE (6,*) 'CDIAG: AN ERROR OCCURRED IN EISCH1. IER =',IER
35
 
        STOP 3
36
 
      ENDIF
37
 
      CALL TRANS (NZ*N,2,Z,Z,H)
38
 
 
39
 
*  UNDO THE REDUCTION OF GENERAL EIGENVALUE PROBLEM BY REDUCC.
40
 
      CALL REBAKK (NS,NZ,N,S,AUX,Z)
41
 
 
42
 
*  STOP TIME COUNT
43
 
      CALL TIMER('CDIAG',2)
44
 
      END
45
 
 
46
 
 
47
 
      SUBROUTINE TRANS (N1,N2,A,B,AUX)
48
 
 
49
 
*  TRANSPOSES MATRIX A(N1,N2) TO B(N2,N1). A AND B MAY BE
50
 
*  THE SAME PHYSICAL ARRAY. WRITTEN BY J.SOLER
51
 
 
52
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
53
 
      DIMENSION A(N1,N2),B(N2,N1),AUX(N1*N2)
54
 
      K=0
55
 
      DO 20 I=1,N2
56
 
        DO 10 J=1,N1
57
 
          K=K+1
58
 
          AUX(K)=A(J,I)
59
 
   10   CONTINUE
60
 
   20 CONTINUE
61
 
      K=0
62
 
      DO 40 I=1,N2
63
 
        DO 30 J=1,N1
64
 
          K=K+1
65
 
          B(I,J)=AUX(K)
66
 
   30   CONTINUE
67
 
   40 CONTINUE
68
 
      END
69
 
 
70
 
 
71
 
      SUBROUTINE PACK (NA,N,A,AP,ISN)
72
 
 
73
 
*  PACKS (ISN=1) OR UNPACKS (ISN=-1) A COMPLEX HERMITIAN
74
 
*  MATRIX TO OR FROM PACKED LOWER MODE.  WRITTEN BY J.SOLER
75
 
 
76
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
77
 
      DIMENSION A(2,NA,N),AP(2,*)
78
 
      PARAMETER (ZERO=0.D0)
79
 
      IF (ISN.GT.0) THEN
80
 
         K=N+1
81
 
         DO 20 I=2,N
82
 
            DO 10 J=I,N
83
 
               AP(1,K)=A(1,J,I)
84
 
               AP(2,K)=A(2,J,I)
85
 
               K=K+1
86
 
  10        CONTINUE
87
 
  20     CONTINUE
88
 
      ELSE
89
 
         K=(N*(N+1))/2
90
 
         DO 40 I=N,2,-1
91
 
            DO 30 J=N,I,-1
92
 
               A(1,J,I)=AP(1,K)
93
 
               A(2,J,I)=AP(2,K)
94
 
               K=K-1
95
 
  30        CONTINUE
96
 
  40     CONTINUE
97
 
         DO 60 I=1,N
98
 
            DO 50 J=1,I-1
99
 
               A(1,J,I)= A(1,I,J)
100
 
               A(2,J,I)=-A(2,I,J)
101
 
  50        CONTINUE
102
 
            A(2,I,I)=ZERO
103
 
  60     CONTINUE
104
 
      ENDIF
105
 
      END
106
 
 
107
 
 
108
 
      SUBROUTINE REDUCC (NH,NS,N,H,S,DL)
109
 
 
110
 
*  GENERALIZATION OF REDUC1 (COMPLEX)
111
 
*  CHOLESKY FACTORIZATION OF S ( S=L*TRANSPOSE(L) )
112
 
*    USING SCHMIDT' ORTHONORMALIZATION METHOD
113
 
*    I)  REPRESENTS A NON-ORTHOGONAL BASIS FUNCTION
114
 
*    I>  REPRESENTS A ORTHONORMAL WAVE FUNCTION
115
 
*  OBTAINED FROM M.METHFESSEL. SLIGHTLY REWRITTEN BY J.SOLER.
116
 
 
117
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
118
 
      DIMENSION H(2,NH,N),S(2,NS,N),DL(N)
119
 
      PARAMETER (ZERO=0.D0,ONE=1.D0)
120
 
 
121
 
*  OBTAIN TRANSFER MATRIX <IJ) FOR I<J USING CHOLESKY FACTORIZATION
122
 
      DO 30 I=1,N
123
 
         DO 20 J=I,N
124
 
            XR=S(1,I,J)
125
 
            XI=S(2,I,J)
126
 
            DO 10 K=1,I-1
127
 
               XR=XR-S(1,J,K)*S(1,I,K)-S(2,J,K)*S(2,I,K)
128
 
               XI=XI+S(1,J,K)*S(2,I,K)-S(2,J,K)*S(1,I,K)
129
 
  10        CONTINUE
130
 
            IF (J.EQ.I) THEN
131
 
               IF (XR.LE.ZERO) THEN
132
 
                  WRITE (6,*) 'REDUCC: MATRIX S NOT POS. DEFINITE. I=',I
133
 
                  STOP
134
 
               END IF
135
 
               DL(I)=DSQRT(XR)
136
 
               Y=ONE/DL(I)
137
 
            ELSE
138
 
               S(1,J,I)=XR*Y
139
 
               S(2,J,I)=XI*Y
140
 
            END IF
141
 
  20     CONTINUE
142
 
  30  CONTINUE
143
 
 
144
 
*  OBTAIN MATRIX <IHJ) FOR I.LE.J
145
 
      DO 60 I=1,N
146
 
         Y=ONE/DL(I)
147
 
         DO 50 J=I,N
148
 
            XR=H(1,J,I)
149
 
            XI=H(2,J,I)
150
 
            DO 40 K=1,I-1
151
 
               XR=XR-H(1,J,K)*S(1,I,K)+H(2,J,K)*S(2,I,K)
152
 
               XI=XI-H(1,J,K)*S(2,I,K)-H(2,J,K)*S(1,I,K)
153
 
  40        CONTINUE
154
 
            H(1,J,I)=XR*Y
155
 
            H(2,J,I)=XI*Y
156
 
  50     CONTINUE
157
 
  60  CONTINUE
158
 
 
159
 
*  OBTAIN MATRIX <IHJ> FOR I.LE.J
160
 
      DO 100 I=1,N
161
 
         DO 90 J=I,N
162
 
            XR=H(1,J,I)
163
 
            XI=H(2,J,I)
164
 
            DO 70 K=1,I-1
165
 
               XR=XR-S(1,J,K)*H(1,I,K)+S(2,J,K)*H(2,I,K)
166
 
               XI=XI+S(1,J,K)*H(2,I,K)+S(2,J,K)*H(1,I,K)
167
 
  70        CONTINUE
168
 
            DO 80 K=I,J-1
169
 
               XR=XR-S(1,J,K)*H(1,K,I)-S(2,J,K)*H(2,K,I)
170
 
               XI=XI-S(1,J,K)*H(2,K,I)+S(2,J,K)*H(1,K,I)
171
 
  80        CONTINUE
172
 
            H(1,J,I)=XR/DL(J)
173
 
            H(2,J,I)=XI/DL(J)
174
 
  90     CONTINUE
175
 
 100  CONTINUE
176
 
 
177
 
*  OBTAIN MATRIX <IHJ> FOR I.GE.J
178
 
      DO 120 I=1,N
179
 
         H(2,I,I)=ZERO
180
 
         DO 110 J=I+1,N
181
 
            H(1,I,J)= H(1,J,I)
182
 
            H(2,I,J)=-H(2,J,I)
183
 
 110     CONTINUE
184
 
 120  CONTINUE
185
 
      END
186
 
 
187
 
 
188
 
      SUBROUTINE REBAKK (NS,NZ,N,S,DL,Z)
189
 
 
190
 
*  BACK TRANSFORMATION OF THE EIGENVECTORS OF THE
191
 
*  GENERAL EIGENVALUE PROBLEM:  A*X = LAMBDA*B*X
192
 
*  THE FORWARD TRANSFORMATION WAS PERFORMED BY ROUTINE 'REDUCC'
193
 
 
194
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
195
 
      DIMENSION S(2,NS,N),Z(2,NZ,N),DL(N)
196
 
      DO 30 I=1,N
197
 
         DO 20 J=N,1,-1
198
 
            XR=Z(1,J,I)
199
 
            XI=Z(2,J,I)
200
 
            DO 10 K=J+1,N
201
 
               XR=XR-S(1,K,J)*Z(1,K,I)+S(2,K,J)*Z(2,K,I)
202
 
               XI=XI-S(1,K,J)*Z(2,K,I)-S(2,K,J)*Z(1,K,I)
203
 
  10        CONTINUE
204
 
            Z(1,J,I)=XR/DL(J)
205
 
            Z(2,J,I)=XI/DL(J)
206
 
  20     CONTINUE
207
 
  30  CONTINUE
208
 
      END
209
 
 
210
 
 
211
 
      SUBROUTINE EISCH1(NM,N,AR,AI,WR,ZR,ZI,IERR,WORK)
212
 
C
213
 
C     ALL EIGENVALUES AND CORRESPONDING EIGENVECTORS OF A COMPLEX
214
 
C     HERMITIAN MATRIX
215
 
C
216
 
*        ADAPTED TO DOUBLE PRECISION BY J.SOLER 30/10/89
217
 
*        FOR SINGLE PRECISION SIMPLY COMMENT OUT NEXT LINE
218
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
219
 
*        PARAMETERS INTRODUCED BY J.SOLER
220
 
      PARAMETER (ZERO=0.D0,ONE=1.D0)
221
 
      DIMENSION AR(NM,NM),AI(NM,NM),WR(N),ZR(NM,NM),ZI(NM,NM),WORK(*)
222
 
      CALL HTRIDI(NM,N,AR,AI,WR,ZI,ZI,WORK)
223
 
      DO 100 I=1,N
224
 
      DO 50 J=1,N
225
 
   50 ZR(I,J)=ZERO
226
 
  100 ZR(I,I)=ONE
227
 
      CALL TQL2C(NM,N,WR,ZI,ZR,IERR)
228
 
      IF(IERR.NE.0) RETURN
229
 
      CALL HTRIBK(NM,N,AR,AI,WORK,N,ZR,ZI)
230
 
      RETURN
231
 
      END
232
 
 
233
 
      SUBROUTINE HTRIDI(NM,N,AR,AI,D,E,E2,TAU)
234
 
*        ADAPTED TO DOUBLE PRECISION BY J.SOLER 30/10/89
235
 
*        FOR SINGLE PRECISION SIMPLY COMMENT OUT NEXT LINE
236
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
237
 
*        PARAMETERS INTRODUCED BY J.SOLER
238
 
      PARAMETER (ZERO=0.D0,ONE=1.D0)
239
 
      DIMENSION AR(NM,N),AI(NM,N),D(N),E(N),E2(N),TAU(2,N)
240
 
*        INTEGER I,J,K,L,N,II,NM,JP1
241
 
*        REAL F,FI,G,GI,H,HH,SI,SCALE
242
 
      TAU(1,N) = ONE
243
 
      TAU(2,N) = ZERO
244
 
      DO 100 I = 1, N
245
 
  100 D(I) = AR(I,I)
246
 
      DO  300 II = 1, N
247
 
         I = N + 1 - II
248
 
         L = I - 1
249
 
         H = ZERO
250
 
         SCALE = ZERO
251
 
         IF (L .LT. 1) GO TO 130
252
 
         DO 120 K = 1, L
253
 
  120    SCALE = SCALE + ABS(AR(I,K)) + ABS(AI(I,K))
254
 
         IF (SCALE .NE. ZERO) GO TO 140
255
 
         TAU(1,L) = ONE
256
 
         TAU(2,L) = ZERO
257
 
  130    E(I) = ZERO
258
 
         E2(I) = ZERO
259
 
         GO TO 290
260
 
  140    DO 150 K = 1, L
261
 
            AR(I,K) = AR(I,K) / SCALE
262
 
            AI(I,K) = AI(I,K) / SCALE
263
 
            H = H + AR(I,K) * AR(I,K) + AI(I,K) * AI(I,K)
264
 
  150    CONTINUE
265
 
         E2(I) = SCALE * SCALE * H
266
 
         G = SQRT(H)
267
 
         E(I) = SCALE * G
268
 
*           NEXT LINE CHANGED BY J.SOLER.
269
 
*        F = CABS(CMPLX(AR(I,L),AI(I,L)))
270
 
         F = SQRT(AR(I,L)**2+AI(I,L)**2)
271
 
         IF (F .EQ. ZERO) GO TO 160
272
 
         TAU(1,L) = (AI(I,L) * TAU(2,I) - AR(I,L) * TAU(1,I)) / F
273
 
         SI = (AR(I,L) * TAU(2,I) + AI(I,L) * TAU(1,I)) / F
274
 
         H = H + F * G
275
 
         G = ONE + G / F
276
 
         AR(I,L) = G * AR(I,L)
277
 
         AI(I,L) = G * AI(I,L)
278
 
         IF (L .EQ. 1) GO TO 270
279
 
         GO TO 170
280
 
  160    TAU(1,L) = -TAU(1,I)
281
 
         SI = TAU(2,I)
282
 
         AR(I,L) = G
283
 
  170    F = ZERO
284
 
         DO 240 J = 1, L
285
 
            G = ZERO
286
 
            GI = ZERO
287
 
            DO 180 K = 1, J
288
 
               G = G + AR(J,K) * AR(I,K) + AI(J,K) * AI(I,K)
289
 
               GI = GI - AR(J,K) * AI(I,K) + AI(J,K) * AR(I,K)
290
 
  180       CONTINUE
291
 
            JP1 = J + 1
292
 
            IF (L .LT. JP1) GO TO 220
293
 
            DO 200 K = JP1, L
294
 
               G = G + AR(K,J) * AR(I,K) - AI(K,J) * AI(I,K)
295
 
               GI = GI - AR(K,J) * AI(I,K) - AI(K,J) * AR(I,K)
296
 
  200       CONTINUE
297
 
  220       E(J) = G / H
298
 
            TAU(2,J) = GI / H
299
 
            F = F + E(J) * AR(I,J) - TAU(2,J) * AI(I,J)
300
 
  240    CONTINUE
301
 
         HH = F / (H + H)
302
 
         DO 260 J = 1, L
303
 
            F = AR(I,J)
304
 
            G = E(J) - HH * F
305
 
            E(J) = G
306
 
            FI = -AI(I,J)
307
 
            GI = TAU(2,J) - HH * FI
308
 
            TAU(2,J) = -GI
309
 
            DO 260 K = 1, J
310
 
               AR(J,K) = AR(J,K) - F * E(K) - G * AR(I,K)
311
 
     X                           + FI * TAU(2,K) + GI * AI(I,K)
312
 
               AI(J,K) = AI(J,K) - F * TAU(2,K) - G * AI(I,K)
313
 
     X                           - FI * E(K) - GI * AR(I,K)
314
 
  260    CONTINUE
315
 
  270    DO 280 K = 1, L
316
 
            AR(I,K) = SCALE * AR(I,K)
317
 
            AI(I,K) = SCALE * AI(I,K)
318
 
  280    CONTINUE
319
 
         TAU(2,L) = -SI
320
 
  290    HH = D(I)
321
 
         D(I) = AR(I,I)
322
 
         AR(I,I) = HH
323
 
         AI(I,I) = SCALE * SCALE * H
324
 
  300 CONTINUE
325
 
      RETURN
326
 
      END
327
 
 
328
 
      SUBROUTINE TQL2C(NM,N,D,E,Z,IERR)
329
 
*        ADAPTED TO DOUBLE PRECISION BY J.SOLER 30/10/89
330
 
*        FOR SINGLE PRECISION SIMPLY COMMENT OUT NEXT LINE
331
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
332
 
*        PARAMETERS INTRODUCED BY J.SOLER
333
 
      PARAMETER (ZERO=0.D0,ONE=1.D0,TWO=2.D0,EPMACH=2.D0**(-23))
334
 
      DIMENSION D(N),E(N),Z(NM,N)
335
 
*        INTEGER I,J,K,L,M,N,II,NM,MML,IERR
336
 
*        VARIABLE MACHEP CHANGED TO PARAMETER EPMACH. J.SOLER.
337
 
*        REAL B,C,F,G,H,P,R,S,MACHEP
338
 
*        MACHEP=TWO**(-23)
339
 
      IERR = 0
340
 
      IF (N .EQ. 1) GO TO 1001
341
 
      DO 100 I = 2, N
342
 
  100 E(I-1) = E(I)
343
 
      F = ZERO
344
 
      B = ZERO
345
 
      E(N) = ZERO
346
 
      DO 240 L = 1, N
347
 
         J = 0
348
 
*        H = MACHEP * (ABS(D(L)) + ABS(E(L)))
349
 
         H = EPMACH * (ABS(D(L)) + ABS(E(L)))
350
 
         IF (B .LT. H) B = H
351
 
         DO 110 M = L, N
352
 
            IF (ABS(E(M)) .LE. B) GO TO 120
353
 
  110    CONTINUE
354
 
  120    IF (M .EQ. L) GO TO 220
355
 
  130    IF (J .EQ. 30) GO TO 1000
356
 
         J = J + 1
357
 
         P = (D(L+1) - D(L)) / (TWO * E(L))
358
 
         R = SQRT(P*P+ONE)
359
 
         H = D(L) - E(L) / (P + SIGN(R,P))
360
 
         DO 140 I = L, N
361
 
  140    D(I) = D(I) - H
362
 
         F = F + H
363
 
         P = D(M)
364
 
         C = ONE
365
 
         S = ZERO
366
 
         MML = M - L
367
 
         DO 200 II = 1, MML
368
 
            I = M - II
369
 
            G = C * E(I)
370
 
            H = C * P
371
 
            IF (ABS(P) .LT. ABS(E(I))) GO TO 150
372
 
            C = E(I) / P
373
 
            R = SQRT(C*C+ONE)
374
 
            E(I+1) = S * P * R
375
 
            S = C / R
376
 
            C = ONE / R
377
 
            GO TO 160
378
 
  150       C = P / E(I)
379
 
            R = SQRT(C*C+ONE)
380
 
            E(I+1) = S * E(I) * R
381
 
            S = ONE / R
382
 
            C = C * S
383
 
  160       P = C * D(I) - S * G
384
 
            D(I+1) = H + S * (C * G + S * D(I))
385
 
            DO 180 K = 1, N
386
 
               H = Z(K,I+1)
387
 
               Z(K,I+1) = S * Z(K,I) + C * H
388
 
               Z(K,I) = C * Z(K,I) - S * H
389
 
  180       CONTINUE
390
 
  200    CONTINUE
391
 
         E(L) = S * P
392
 
         D(L) = C * P
393
 
         IF (ABS(E(L)) .GT. B) GO TO 130
394
 
  220    D(L) = D(L) + F
395
 
  240 CONTINUE
396
 
      DO 300 II = 2, N
397
 
         I = II - 1
398
 
         K = I
399
 
         P = D(I)
400
 
         DO 260 J = II, N
401
 
            IF (D(J) .GE. P) GO TO 260
402
 
            K = J
403
 
            P = D(J)
404
 
  260    CONTINUE
405
 
         IF (K .EQ. I) GO TO 300
406
 
         D(K) = D(I)
407
 
         D(I) = P
408
 
         DO 280 J = 1, N
409
 
            P = Z(J,I)
410
 
            Z(J,I) = Z(J,K)
411
 
            Z(J,K) = P
412
 
  280    CONTINUE
413
 
  300 CONTINUE
414
 
      GO TO 1001
415
 
 1000 IERR = L
416
 
 1001 RETURN
417
 
      END
418
 
 
419
 
      SUBROUTINE HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)
420
 
*        ADAPTED TO DOUBLE PRECISION BY J.SOLER 30/10/89
421
 
*        FOR SINGLE PRECISION SIMPLY COMMENT OUT NEXT LINE
422
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
423
 
*        PARAMETER INTRODUCED BY J.SOLER
424
 
      PARAMETER (ZERO=0.D0)
425
 
      DIMENSION AR(NM,N),AI(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M)
426
 
*        REAL H,S,SI
427
 
*        INTEGER I,J,K,L,M,N,NM
428
 
      DO 50 K = 1, N
429
 
         DO 50 J = 1, M
430
 
            ZI(K,J) = (- ZR(K,J)) * TAU(2,K)
431
 
            ZR(K,J) = ZR(K,J) * TAU(1,K)
432
 
   50 CONTINUE
433
 
      IF (N .EQ. 1) GO TO 200
434
 
      DO 140 I = 2, N
435
 
         L = I - 1
436
 
         H = AI(I,I)
437
 
         IF (H .EQ. ZERO) GO TO 140
438
 
         DO 130 J = 1, M
439
 
            S = ZERO
440
 
            SI = ZERO
441
 
            DO 110 K = 1, L
442
 
               S = S + AR(I,K) * ZR(K,J) - AI(I,K) * ZI(K,J)
443
 
               SI = SI + AR(I,K) * ZI(K,J) + AI(I,K) * ZR(K,J)
444
 
  110       CONTINUE
445
 
            S = S / H
446
 
            SI = SI / H
447
 
            DO 120 K = 1, L
448
 
               ZR(K,J) = ZR(K,J) - S * AR(I,K) - SI * AI(I,K)
449
 
               ZI(K,J) = ZI(K,J) - SI * AR(I,K) + S * AI(I,K)
450
 
  120       CONTINUE
451
 
  130    CONTINUE
452
 
  140 CONTINUE
453
 
  200 RETURN
454
 
      END