1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
|
C
C***********************************************************************
C
SUBROUTINE DLABCM(N, NBAND, NL, NR, A, EIGVAL,
1 LDE, EIGVEC, ATOL, ARTOL, BOUND, ATEMP, D, VTEMP)
C
C THIS SUBROUTINE ORGANIZES THE CALCULATION OF THE EIGENVALUES
C FOR THE BNDEIG PACKAGE. EIGENVALUES ARE COMPUTED BY
C A MODIFIED RAYLEIGH QUOTIENT ITERATION. THE EIGENVALUE COUNT
C OBTAINED BY EACH FACTORIZATION IS USED TO OCCASIONALLY OVERRIDE
C THE COMPUTED RAYLEIGH QUOTIENT WITH A DIFFERENT SHIFT TO
C INSURE CONVERGENCE TO THE DESIRED EIGENVALUES.
C
C FORMAL PARAMETERS.
C
INTEGER N, NBAND, NL, NR, LDE
DOUBLE PRECISION A(NBAND,1), EIGVAL(1),
1 EIGVEC(LDE,1), ATOL, ARTOL, BOUND(2,1), ATEMP(1),
2 D(1), VTEMP(1)
C
C
C LOCAL VARIABLES
C
LOGICAL FLAG
INTEGER I, J, L, M, NUML, NUMVEC, NVAL
DOUBLE PRECISION ERRB, GAP, RESID, RQ, SIGMA, VNORM
C
C
C FUNCTIONS CALLED
C
INTEGER MIN0
DOUBLE PRECISION DMAX1, DMIN1, DDOT, DNRM2
C
C SUBROUTINES CALLED
C
C DLABAX, DLABFC, DLARAN, DAXPY, DCOPY, DSCAL
C
C REPLACE ZERO VECTORS BY RANDOM
C
NVAL = NR - NL + 1
FLAG = .FALSE.
DO 5 I = 1, NVAL
IF(DDOT(N, EIGVEC(1,I), 1, EIGVEC(1,I), 1) .EQ. 0.0D0)
1 CALL DLARAN(N,EIGVEC(1,I))
5 CONTINUE
C
C LOOP OVER EIGENVALUES
C
SIGMA = BOUND(2,NVAL+1)
DO 400 J = 1, NVAL
NUML = J
C
C PREPARE TO COMPUTE FIRST RAYLEIGH QUOTIENT
C
10 CALL DLABAX(N, NBAND, A, EIGVEC(1,J), VTEMP)
VNORM = DNRM2(N, VTEMP, 1)
IF(VNORM .EQ. 0.0D0) GO TO 20
CALL DSCAL(N, 1.0D0/VNORM, VTEMP, 1)
CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1)
CALL DAXPY(N, -SIGMA, EIGVEC(1,J), 1, VTEMP, 1)
C
C LOOP OVER SHIFTS
C
C COMPUTE RAYLEIGH QUOTIENT, RESIDUAL NORM, AND CURRENT TOLERANCE
C
20 VNORM = DNRM2(N, EIGVEC(1,J), 1)
IF(VNORM .NE. 0.0D0) GO TO 30
CALL DLARAN(N, EIGVEC(1,J))
GO TO 10
C
30 RQ = SIGMA + DDOT(N, EIGVEC(1,J), 1, VTEMP, 1)
1 /VNORM/VNORM
CALL DAXPY(N, SIGMA-RQ, EIGVEC(1,J), 1, VTEMP, 1)
RESID = DMAX1(ATOL, DNRM2(N, VTEMP, 1)/VNORM)
CALL DSCAL(N, 1.0/VNORM, EIGVEC(1,J), 1)
C
C ACCEPT EIGENVALUE IF THE INTERVAL IS SMALL ENOUGH
C
IF(BOUND(2,J+1) - BOUND(1,J+1) .LT. 3.0D0*ATOL) GO TO 300
C
C COMPUTE MINIMAL ERROR BOUND
C
ERRB = RESID
GAP = DMIN1(BOUND(1,J+2) - RQ, RQ - BOUND(2,J))
IF(GAP .GT. RESID) ERRB = DMAX1(ATOL, RESID*RESID/GAP)
C
C TENTATIVE NEW SHIFT
C
SIGMA = 0.5D0*(BOUND(1,J+1) + BOUND(2,J+1))
C
C CHECK FOR TERMINALTION
C
IF(RESID .GT. 2.0D0*ATOL) GO TO 40
IF(RQ - ERRB .GT. BOUND(2,J) .AND.
1 RQ + ERRB .LT. BOUND(1,J+2)) GO TO 310
C
C RQ IS TO THE LEFT OF THE INTERVAL
C
40 IF(RQ .GE. BOUND(1,J+1)) GO TO 50
IF(RQ - ERRB .GT. BOUND(2,J)) GO TO 100
IF(RQ + ERRB .LT. BOUND(1,J+1)) CALL DLARAN(N,EIGVEC(1,J))
GO TO 200
C
C RQ IS TO THE RIGHT OF THE INTERVAL
C
50 IF(RQ .LE. BOUND(2,J+1)) GO TO 100
IF(RQ + ERRB .LT. BOUND(1,J+2)) GO TO 100
C
C SAVE THE REJECTED VECTOR IF INDICATED
C
IF(RQ - ERRB .LE. BOUND(2,J+1)) GO TO 200
DO 60 I = J, NVAL
IF(BOUND(2,I+1) .GT. RQ) GO TO 70
60 CONTINUE
GO TO 80
C
70 CALL DCOPY(N, EIGVEC(1,J), 1, EIGVEC(1,I), 1)
C
80 CALL DLARAN(N, EIGVEC(1,J))
GO TO 200
C
C PERTURB RQ TOWARD THE MIDDLE
C
100 IF(SIGMA .LT. RQ) SIGMA = DMAX1(SIGMA, RQ-ERRB)
IF(SIGMA .GT. RQ) SIGMA = DMIN1(SIGMA, RQ+ERRB)
C
C FACTOR AND SOLVE
C
200 DO 210 I = J, NVAL
IF(SIGMA .LT. BOUND(1,I+1)) GO TO 220
210 CONTINUE
I = NVAL + 1
220 NUMVEC = I - J
NUMVEC = MIN0(NUMVEC, NBAND + 2)
IF(RESID .LT. ARTOL) NUMVEC = MIN0(1,NUMVEC)
CALL DCOPY(N, EIGVEC(1,J), 1, VTEMP, 1)
CALL DLABFC(N, NBAND, A, SIGMA, NUMVEC, LDE,
1 EIGVEC(1,J), NUML, 2*NBAND-1, ATEMP, D, ATOL)
C
C PARTIALLY SCALE EXTRA VECTORS TO PREVENT UNDERFLOW OR OVERFLOW
C
IF(NUMVEC .EQ. 1) GO TO 227
L = NUMVEC - 1
DO 225 I = 1,L
M = J + I
CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,M), 1)
225 CONTINUE
C
C UPDATE INTERVALS
C
227 NUML = NUML - NL + 1
IF(NUML .GE. 0) BOUND(2,1) = DMIN1(BOUND(2,1), SIGMA)
DO 230 I = J, NVAL
IF(SIGMA .LT. BOUND(1,I+1)) GO TO 20
IF(NUML .LT. I) BOUND(1,I+1) = SIGMA
IF(NUML .GE. I) BOUND(2,I+1) = SIGMA
230 CONTINUE
IF(NUML .LT. NVAL + 1) BOUND(1,NVAL+2) = DMAX1(SIGMA,
1 BOUND(1,NVAL+2))
GO TO 20
C
C ACCEPT AN EIGENPAIR
C
300 CALL DLARAN(N, EIGVEC(1,J))
FLAG = .TRUE.
GO TO 310
C
305 FLAG = .FALSE.
RQ = 0.5D0*(BOUND(1,J+1) + BOUND(2,J+1))
CALL DLABFC(N, NBAND, A, RQ, NUMVEC, LDE,
1 EIGVEC(1,J), NUML, 2*NBAND-1, ATEMP, D, ATOL)
VNORM = DNRM2(N, EIGVEC(1,J), 1)
IF(VNORM .NE. 0.0) CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1)
C
C ORTHOGONALIZE THE NEW EIGENVECTOR AGAINST THE OLD ONES
C
310 EIGVAL(J) = RQ
IF(J .EQ. 1) GO TO 330
M = J - 1
DO 320 I = 1, M
CALL DAXPY(N, -DDOT(N,EIGVEC(1,I),1,EIGVEC(1,J),1),
1 EIGVEC(1,I), 1, EIGVEC(1,J), 1)
320 CONTINUE
330 VNORM = DNRM2(N, EIGVEC(1,J), 1)
IF(VNORM .EQ. 0.0D0) GO TO 305
CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1)
C
C ORTHOGONALIZE LATER VECTORS AGAINST THE CONVERGED ONE
C
IF(FLAG) GO TO 305
IF(J .EQ. NVAL) RETURN
M = J + 1
DO 340 I = M, NVAL
CALL DAXPY(N, -DDOT(N,EIGVEC(1,J),1,EIGVEC(1,I),1),
1 EIGVEC(1,J), 1, EIGVEC(1,I), 1)
340 CONTINUE
400 CONTINUE
RETURN
C
500 CONTINUE
END
|
