Mercurial > octave-nkf
comparison libcruft/arpack/src/dsaup2.f @ 12274:9f5d2ef078e8 release-3-4-x
import ARPACK sources to libcruft from Debian package libarpack2 2.1+parpack96.dfsg-3+b1
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 28 Jan 2011 14:04:33 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
12273:83133b5bf392 | 12274:9f5d2ef078e8 |
---|---|
1 c----------------------------------------------------------------------- | |
2 c\BeginDoc | |
3 c | |
4 c\Name: dsaup2 | |
5 c | |
6 c\Description: | |
7 c Intermediate level interface called by dsaupd. | |
8 c | |
9 c\Usage: | |
10 c call dsaup2 | |
11 c ( IDO, BMAT, N, WHICH, NEV, NP, TOL, RESID, MODE, IUPD, | |
12 c ISHIFT, MXITER, V, LDV, H, LDH, RITZ, BOUNDS, Q, LDQ, WORKL, | |
13 c IPNTR, WORKD, INFO ) | |
14 c | |
15 c\Arguments | |
16 c | |
17 c IDO, BMAT, N, WHICH, NEV, TOL, RESID: same as defined in dsaupd. | |
18 c MODE, ISHIFT, MXITER: see the definition of IPARAM in dsaupd. | |
19 c | |
20 c NP Integer. (INPUT/OUTPUT) | |
21 c Contains the number of implicit shifts to apply during | |
22 c each Arnoldi/Lanczos iteration. | |
23 c If ISHIFT=1, NP is adjusted dynamically at each iteration | |
24 c to accelerate convergence and prevent stagnation. | |
25 c This is also roughly equal to the number of matrix-vector | |
26 c products (involving the operator OP) per Arnoldi iteration. | |
27 c The logic for adjusting is contained within the current | |
28 c subroutine. | |
29 c If ISHIFT=0, NP is the number of shifts the user needs | |
30 c to provide via reverse comunication. 0 < NP < NCV-NEV. | |
31 c NP may be less than NCV-NEV since a leading block of the current | |
32 c upper Tridiagonal matrix has split off and contains "unwanted" | |
33 c Ritz values. | |
34 c Upon termination of the IRA iteration, NP contains the number | |
35 c of "converged" wanted Ritz values. | |
36 c | |
37 c IUPD Integer. (INPUT) | |
38 c IUPD .EQ. 0: use explicit restart instead implicit update. | |
39 c IUPD .NE. 0: use implicit update. | |
40 c | |
41 c V Double precision N by (NEV+NP) array. (INPUT/OUTPUT) | |
42 c The Lanczos basis vectors. | |
43 c | |
44 c LDV Integer. (INPUT) | |
45 c Leading dimension of V exactly as declared in the calling | |
46 c program. | |
47 c | |
48 c H Double precision (NEV+NP) by 2 array. (OUTPUT) | |
49 c H is used to store the generated symmetric tridiagonal matrix | |
50 c The subdiagonal is stored in the first column of H starting | |
51 c at H(2,1). The main diagonal is stored in the arscnd column | |
52 c of H starting at H(1,2). If dsaup2 converges store the | |
53 c B-norm of the final residual vector in H(1,1). | |
54 c | |
55 c LDH Integer. (INPUT) | |
56 c Leading dimension of H exactly as declared in the calling | |
57 c program. | |
58 c | |
59 c RITZ Double precision array of length NEV+NP. (OUTPUT) | |
60 c RITZ(1:NEV) contains the computed Ritz values of OP. | |
61 c | |
62 c BOUNDS Double precision array of length NEV+NP. (OUTPUT) | |
63 c BOUNDS(1:NEV) contain the error bounds corresponding to RITZ. | |
64 c | |
65 c Q Double precision (NEV+NP) by (NEV+NP) array. (WORKSPACE) | |
66 c Private (replicated) work array used to accumulate the | |
67 c rotation in the shift application step. | |
68 c | |
69 c LDQ Integer. (INPUT) | |
70 c Leading dimension of Q exactly as declared in the calling | |
71 c program. | |
72 c | |
73 c WORKL Double precision array of length at least 3*(NEV+NP). (INPUT/WORKSPACE) | |
74 c Private (replicated) array on each PE or array allocated on | |
75 c the front end. It is used in the computation of the | |
76 c tridiagonal eigenvalue problem, the calculation and | |
77 c application of the shifts and convergence checking. | |
78 c If ISHIFT .EQ. O and IDO .EQ. 3, the first NP locations | |
79 c of WORKL are used in reverse communication to hold the user | |
80 c supplied shifts. | |
81 c | |
82 c IPNTR Integer array of length 3. (OUTPUT) | |
83 c Pointer to mark the starting locations in the WORKD for | |
84 c vectors used by the Lanczos iteration. | |
85 c ------------------------------------------------------------- | |
86 c IPNTR(1): pointer to the current operand vector X. | |
87 c IPNTR(2): pointer to the current result vector Y. | |
88 c IPNTR(3): pointer to the vector B * X when used in one of | |
89 c the spectral transformation modes. X is the current | |
90 c operand. | |
91 c ------------------------------------------------------------- | |
92 c | |
93 c WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION) | |
94 c Distributed array to be used in the basic Lanczos iteration | |
95 c for reverse communication. The user should not use WORKD | |
96 c as temporary workspace during the iteration !!!!!!!!!! | |
97 c See Data Distribution Note in dsaupd. | |
98 c | |
99 c INFO Integer. (INPUT/OUTPUT) | |
100 c If INFO .EQ. 0, a randomly initial residual vector is used. | |
101 c If INFO .NE. 0, RESID contains the initial residual vector, | |
102 c possibly from a previous run. | |
103 c Error flag on output. | |
104 c = 0: Normal return. | |
105 c = 1: All possible eigenvalues of OP has been found. | |
106 c NP returns the size of the invariant subspace | |
107 c spanning the operator OP. | |
108 c = 2: No shifts could be applied. | |
109 c = -8: Error return from trid. eigenvalue calculation; | |
110 c This should never happen. | |
111 c = -9: Starting vector is zero. | |
112 c = -9999: Could not build an Lanczos factorization. | |
113 c Size that was built in returned in NP. | |
114 c | |
115 c\EndDoc | |
116 c | |
117 c----------------------------------------------------------------------- | |
118 c | |
119 c\BeginLib | |
120 c | |
121 c\References: | |
122 c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in | |
123 c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), | |
124 c pp 357-385. | |
125 c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly | |
126 c Restarted Arnoldi Iteration", Rice University Technical Report | |
127 c TR95-13, Department of Computational and Applied Mathematics. | |
128 c 3. B.N. Parlett, "The Symmetric Eigenvalue Problem". Prentice-Hall, | |
129 c 1980. | |
130 c 4. B.N. Parlett, B. Nour-Omid, "Towards a Black Box Lanczos Program", | |
131 c Computer Physics Communications, 53 (1989), pp 169-179. | |
132 c 5. B. Nour-Omid, B.N. Parlett, T. Ericson, P.S. Jensen, "How to | |
133 c Implement the Spectral Transformation", Math. Comp., 48 (1987), | |
134 c pp 663-673. | |
135 c 6. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos | |
136 c Algorithm for Solving Sparse Symmetric Generalized Eigenproblems", | |
137 c SIAM J. Matr. Anal. Apps., January (1993). | |
138 c 7. L. Reichel, W.B. Gragg, "Algorithm 686: FORTRAN Subroutines | |
139 c for Updating the QR decomposition", ACM TOMS, December 1990, | |
140 c Volume 16 Number 4, pp 369-377. | |
141 c | |
142 c\Routines called: | |
143 c dgetv0 ARPACK initial vector generation routine. | |
144 c dsaitr ARPACK Lanczos factorization routine. | |
145 c dsapps ARPACK application of implicit shifts routine. | |
146 c dsconv ARPACK convergence of Ritz values routine. | |
147 c dseigt ARPACK compute Ritz values and error bounds routine. | |
148 c dsgets ARPACK reorder Ritz values and error bounds routine. | |
149 c dsortr ARPACK sorting routine. | |
150 c ivout ARPACK utility routine that prints integers. | |
151 c arscnd ARPACK utility routine for timing. | |
152 c dvout ARPACK utility routine that prints vectors. | |
153 c dlamch LAPACK routine that determines machine constants. | |
154 c dcopy Level 1 BLAS that copies one vector to another. | |
155 c ddot Level 1 BLAS that computes the scalar product of two vectors. | |
156 c dnrm2 Level 1 BLAS that computes the norm of a vector. | |
157 c dscal Level 1 BLAS that scales a vector. | |
158 c dswap Level 1 BLAS that swaps two vectors. | |
159 c | |
160 c\Author | |
161 c Danny Sorensen Phuong Vu | |
162 c Richard Lehoucq CRPC / Rice University | |
163 c Dept. of Computational & Houston, Texas | |
164 c Applied Mathematics | |
165 c Rice University | |
166 c Houston, Texas | |
167 c | |
168 c\Revision history: | |
169 c 12/15/93: Version ' 2.4' | |
170 c xx/xx/95: Version ' 2.4'. (R.B. Lehoucq) | |
171 c | |
172 c\SCCS Information: @(#) | |
173 c FILE: saup2.F SID: 2.7 DATE OF SID: 5/19/98 RELEASE: 2 | |
174 c | |
175 c\EndLib | |
176 c | |
177 c----------------------------------------------------------------------- | |
178 c | |
179 subroutine dsaup2 | |
180 & ( ido, bmat, n, which, nev, np, tol, resid, mode, iupd, | |
181 & ishift, mxiter, v, ldv, h, ldh, ritz, bounds, | |
182 & q, ldq, workl, ipntr, workd, info ) | |
183 c | |
184 c %----------------------------------------------------% | |
185 c | Include files for debugging and timing information | | |
186 c %----------------------------------------------------% | |
187 c | |
188 include 'debug.h' | |
189 include 'stat.h' | |
190 c | |
191 c %------------------% | |
192 c | Scalar Arguments | | |
193 c %------------------% | |
194 c | |
195 character bmat*1, which*2 | |
196 integer ido, info, ishift, iupd, ldh, ldq, ldv, mxiter, | |
197 & n, mode, nev, np | |
198 Double precision | |
199 & tol | |
200 c | |
201 c %-----------------% | |
202 c | Array Arguments | | |
203 c %-----------------% | |
204 c | |
205 integer ipntr(3) | |
206 Double precision | |
207 & bounds(nev+np), h(ldh,2), q(ldq,nev+np), resid(n), | |
208 & ritz(nev+np), v(ldv,nev+np), workd(3*n), | |
209 & workl(3*(nev+np)) | |
210 c | |
211 c %------------% | |
212 c | Parameters | | |
213 c %------------% | |
214 c | |
215 Double precision | |
216 & one, zero | |
217 parameter (one = 1.0D+0, zero = 0.0D+0) | |
218 c | |
219 c %---------------% | |
220 c | Local Scalars | | |
221 c %---------------% | |
222 c | |
223 character wprime*2 | |
224 logical cnorm, getv0, initv, update, ushift | |
225 integer ierr, iter, j, kplusp, msglvl, nconv, nevbef, nev0, | |
226 & np0, nptemp, nevd2, nevm2, kp(3) | |
227 Double precision | |
228 & rnorm, temp, eps23 | |
229 save cnorm, getv0, initv, update, ushift, | |
230 & iter, kplusp, msglvl, nconv, nev0, np0, | |
231 & rnorm, eps23 | |
232 c | |
233 c %----------------------% | |
234 c | External Subroutines | | |
235 c %----------------------% | |
236 c | |
237 external dcopy, dgetv0, dsaitr, dscal, dsconv, dseigt, dsgets, | |
238 & dsapps, dsortr, dvout, ivout, arscnd, dswap | |
239 c | |
240 c %--------------------% | |
241 c | External Functions | | |
242 c %--------------------% | |
243 c | |
244 Double precision | |
245 & ddot, dnrm2, dlamch | |
246 external ddot, dnrm2, dlamch | |
247 c | |
248 c %---------------------% | |
249 c | Intrinsic Functions | | |
250 c %---------------------% | |
251 c | |
252 intrinsic min | |
253 c | |
254 c %-----------------------% | |
255 c | Executable Statements | | |
256 c %-----------------------% | |
257 c | |
258 if (ido .eq. 0) then | |
259 c | |
260 c %-------------------------------% | |
261 c | Initialize timing statistics | | |
262 c | & message level for debugging | | |
263 c %-------------------------------% | |
264 c | |
265 call arscnd (t0) | |
266 msglvl = msaup2 | |
267 c | |
268 c %---------------------------------% | |
269 c | Set machine dependent constant. | | |
270 c %---------------------------------% | |
271 c | |
272 eps23 = dlamch('Epsilon-Machine') | |
273 eps23 = eps23**(2.0D+0/3.0D+0) | |
274 c | |
275 c %-------------------------------------% | |
276 c | nev0 and np0 are integer variables | | |
277 c | hold the initial values of NEV & NP | | |
278 c %-------------------------------------% | |
279 c | |
280 nev0 = nev | |
281 np0 = np | |
282 c | |
283 c %-------------------------------------% | |
284 c | kplusp is the bound on the largest | | |
285 c | Lanczos factorization built. | | |
286 c | nconv is the current number of | | |
287 c | "converged" eigenvlues. | | |
288 c | iter is the counter on the current | | |
289 c | iteration step. | | |
290 c %-------------------------------------% | |
291 c | |
292 kplusp = nev0 + np0 | |
293 nconv = 0 | |
294 iter = 0 | |
295 c | |
296 c %--------------------------------------------% | |
297 c | Set flags for computing the first NEV steps | | |
298 c | of the Lanczos factorization. | | |
299 c %--------------------------------------------% | |
300 c | |
301 getv0 = .true. | |
302 update = .false. | |
303 ushift = .false. | |
304 cnorm = .false. | |
305 c | |
306 if (info .ne. 0) then | |
307 c | |
308 c %--------------------------------------------% | |
309 c | User provides the initial residual vector. | | |
310 c %--------------------------------------------% | |
311 c | |
312 initv = .true. | |
313 info = 0 | |
314 else | |
315 initv = .false. | |
316 end if | |
317 end if | |
318 c | |
319 c %---------------------------------------------% | |
320 c | Get a possibly random starting vector and | | |
321 c | force it into the range of the operator OP. | | |
322 c %---------------------------------------------% | |
323 c | |
324 10 continue | |
325 c | |
326 if (getv0) then | |
327 call dgetv0 (ido, bmat, 1, initv, n, 1, v, ldv, resid, rnorm, | |
328 & ipntr, workd, info) | |
329 c | |
330 if (ido .ne. 99) go to 9000 | |
331 c | |
332 if (rnorm .eq. zero) then | |
333 c | |
334 c %-----------------------------------------% | |
335 c | The initial vector is zero. Error exit. | | |
336 c %-----------------------------------------% | |
337 c | |
338 info = -9 | |
339 go to 1200 | |
340 end if | |
341 getv0 = .false. | |
342 ido = 0 | |
343 end if | |
344 c | |
345 c %------------------------------------------------------------% | |
346 c | Back from reverse communication: continue with update step | | |
347 c %------------------------------------------------------------% | |
348 c | |
349 if (update) go to 20 | |
350 c | |
351 c %-------------------------------------------% | |
352 c | Back from computing user specified shifts | | |
353 c %-------------------------------------------% | |
354 c | |
355 if (ushift) go to 50 | |
356 c | |
357 c %-------------------------------------% | |
358 c | Back from computing residual norm | | |
359 c | at the end of the current iteration | | |
360 c %-------------------------------------% | |
361 c | |
362 if (cnorm) go to 100 | |
363 c | |
364 c %----------------------------------------------------------% | |
365 c | Compute the first NEV steps of the Lanczos factorization | | |
366 c %----------------------------------------------------------% | |
367 c | |
368 call dsaitr (ido, bmat, n, 0, nev0, mode, resid, rnorm, v, ldv, | |
369 & h, ldh, ipntr, workd, info) | |
370 c | |
371 c %---------------------------------------------------% | |
372 c | ido .ne. 99 implies use of reverse communication | | |
373 c | to compute operations involving OP and possibly B | | |
374 c %---------------------------------------------------% | |
375 c | |
376 if (ido .ne. 99) go to 9000 | |
377 c | |
378 if (info .gt. 0) then | |
379 c | |
380 c %-----------------------------------------------------% | |
381 c | dsaitr was unable to build an Lanczos factorization | | |
382 c | of length NEV0. INFO is returned with the size of | | |
383 c | the factorization built. Exit main loop. | | |
384 c %-----------------------------------------------------% | |
385 c | |
386 np = info | |
387 mxiter = iter | |
388 info = -9999 | |
389 go to 1200 | |
390 end if | |
391 c | |
392 c %--------------------------------------------------------------% | |
393 c | | | |
394 c | M A I N LANCZOS I T E R A T I O N L O O P | | |
395 c | Each iteration implicitly restarts the Lanczos | | |
396 c | factorization in place. | | |
397 c | | | |
398 c %--------------------------------------------------------------% | |
399 c | |
400 1000 continue | |
401 c | |
402 iter = iter + 1 | |
403 c | |
404 if (msglvl .gt. 0) then | |
405 call ivout (logfil, 1, iter, ndigit, | |
406 & '_saup2: **** Start of major iteration number ****') | |
407 end if | |
408 if (msglvl .gt. 1) then | |
409 call ivout (logfil, 1, nev, ndigit, | |
410 & '_saup2: The length of the current Lanczos factorization') | |
411 call ivout (logfil, 1, np, ndigit, | |
412 & '_saup2: Extend the Lanczos factorization by') | |
413 end if | |
414 c | |
415 c %------------------------------------------------------------% | |
416 c | Compute NP additional steps of the Lanczos factorization. | | |
417 c %------------------------------------------------------------% | |
418 c | |
419 ido = 0 | |
420 20 continue | |
421 update = .true. | |
422 c | |
423 call dsaitr (ido, bmat, n, nev, np, mode, resid, rnorm, v, | |
424 & ldv, h, ldh, ipntr, workd, info) | |
425 c | |
426 c %---------------------------------------------------% | |
427 c | ido .ne. 99 implies use of reverse communication | | |
428 c | to compute operations involving OP and possibly B | | |
429 c %---------------------------------------------------% | |
430 c | |
431 if (ido .ne. 99) go to 9000 | |
432 c | |
433 if (info .gt. 0) then | |
434 c | |
435 c %-----------------------------------------------------% | |
436 c | dsaitr was unable to build an Lanczos factorization | | |
437 c | of length NEV0+NP0. INFO is returned with the size | | |
438 c | of the factorization built. Exit main loop. | | |
439 c %-----------------------------------------------------% | |
440 c | |
441 np = info | |
442 mxiter = iter | |
443 info = -9999 | |
444 go to 1200 | |
445 end if | |
446 update = .false. | |
447 c | |
448 if (msglvl .gt. 1) then | |
449 call dvout (logfil, 1, rnorm, ndigit, | |
450 & '_saup2: Current B-norm of residual for factorization') | |
451 end if | |
452 c | |
453 c %--------------------------------------------------------% | |
454 c | Compute the eigenvalues and corresponding error bounds | | |
455 c | of the current symmetric tridiagonal matrix. | | |
456 c %--------------------------------------------------------% | |
457 c | |
458 call dseigt (rnorm, kplusp, h, ldh, ritz, bounds, workl, ierr) | |
459 c | |
460 if (ierr .ne. 0) then | |
461 info = -8 | |
462 go to 1200 | |
463 end if | |
464 c | |
465 c %----------------------------------------------------% | |
466 c | Make a copy of eigenvalues and corresponding error | | |
467 c | bounds obtained from _seigt. | | |
468 c %----------------------------------------------------% | |
469 c | |
470 call dcopy(kplusp, ritz, 1, workl(kplusp+1), 1) | |
471 call dcopy(kplusp, bounds, 1, workl(2*kplusp+1), 1) | |
472 c | |
473 c %---------------------------------------------------% | |
474 c | Select the wanted Ritz values and their bounds | | |
475 c | to be used in the convergence test. | | |
476 c | The selection is based on the requested number of | | |
477 c | eigenvalues instead of the current NEV and NP to | | |
478 c | prevent possible misconvergence. | | |
479 c | * Wanted Ritz values := RITZ(NP+1:NEV+NP) | | |
480 c | * Shifts := RITZ(1:NP) := WORKL(1:NP) | | |
481 c %---------------------------------------------------% | |
482 c | |
483 nev = nev0 | |
484 np = np0 | |
485 call dsgets (ishift, which, nev, np, ritz, bounds, workl) | |
486 c | |
487 c %-------------------% | |
488 c | Convergence test. | | |
489 c %-------------------% | |
490 c | |
491 call dcopy (nev, bounds(np+1), 1, workl(np+1), 1) | |
492 call dsconv (nev, ritz(np+1), workl(np+1), tol, nconv) | |
493 c | |
494 if (msglvl .gt. 2) then | |
495 kp(1) = nev | |
496 kp(2) = np | |
497 kp(3) = nconv | |
498 call ivout (logfil, 3, kp, ndigit, | |
499 & '_saup2: NEV, NP, NCONV are') | |
500 call dvout (logfil, kplusp, ritz, ndigit, | |
501 & '_saup2: The eigenvalues of H') | |
502 call dvout (logfil, kplusp, bounds, ndigit, | |
503 & '_saup2: Ritz estimates of the current NCV Ritz values') | |
504 end if | |
505 c | |
506 c %---------------------------------------------------------% | |
507 c | Count the number of unwanted Ritz values that have zero | | |
508 c | Ritz estimates. If any Ritz estimates are equal to zero | | |
509 c | then a leading block of H of order equal to at least | | |
510 c | the number of Ritz values with zero Ritz estimates has | | |
511 c | split off. None of these Ritz values may be removed by | | |
512 c | shifting. Decrease NP the number of shifts to apply. If | | |
513 c | no shifts may be applied, then prepare to exit | | |
514 c %---------------------------------------------------------% | |
515 c | |
516 nptemp = np | |
517 do 30 j=1, nptemp | |
518 if (bounds(j) .eq. zero) then | |
519 np = np - 1 | |
520 nev = nev + 1 | |
521 end if | |
522 30 continue | |
523 c | |
524 if ( (nconv .ge. nev0) .or. | |
525 & (iter .gt. mxiter) .or. | |
526 & (np .eq. 0) ) then | |
527 c | |
528 c %------------------------------------------------% | |
529 c | Prepare to exit. Put the converged Ritz values | | |
530 c | and corresponding bounds in RITZ(1:NCONV) and | | |
531 c | BOUNDS(1:NCONV) respectively. Then sort. Be | | |
532 c | careful when NCONV > NP since we don't want to | | |
533 c | swap overlapping locations. | | |
534 c %------------------------------------------------% | |
535 c | |
536 if (which .eq. 'BE') then | |
537 c | |
538 c %-----------------------------------------------------% | |
539 c | Both ends of the spectrum are requested. | | |
540 c | Sort the eigenvalues into algebraically decreasing | | |
541 c | order first then swap low end of the spectrum next | | |
542 c | to high end in appropriate locations. | | |
543 c | NOTE: when np < floor(nev/2) be careful not to swap | | |
544 c | overlapping locations. | | |
545 c %-----------------------------------------------------% | |
546 c | |
547 wprime = 'SA' | |
548 call dsortr (wprime, .true., kplusp, ritz, bounds) | |
549 nevd2 = nev0 / 2 | |
550 nevm2 = nev0 - nevd2 | |
551 if ( nev .gt. 1 ) then | |
552 call dswap ( min(nevd2,np), ritz(nevm2+1), 1, | |
553 & ritz( max(kplusp-nevd2+1,kplusp-np+1) ), 1) | |
554 call dswap ( min(nevd2,np), bounds(nevm2+1), 1, | |
555 & bounds( max(kplusp-nevd2+1,kplusp-np+1)), 1) | |
556 end if | |
557 c | |
558 else | |
559 c | |
560 c %--------------------------------------------------% | |
561 c | LM, SM, LA, SA case. | | |
562 c | Sort the eigenvalues of H into the an order that | | |
563 c | is opposite to WHICH, and apply the resulting | | |
564 c | order to BOUNDS. The eigenvalues are sorted so | | |
565 c | that the wanted part are always within the first | | |
566 c | NEV locations. | | |
567 c %--------------------------------------------------% | |
568 c | |
569 if (which .eq. 'LM') wprime = 'SM' | |
570 if (which .eq. 'SM') wprime = 'LM' | |
571 if (which .eq. 'LA') wprime = 'SA' | |
572 if (which .eq. 'SA') wprime = 'LA' | |
573 c | |
574 call dsortr (wprime, .true., kplusp, ritz, bounds) | |
575 c | |
576 end if | |
577 c | |
578 c %--------------------------------------------------% | |
579 c | Scale the Ritz estimate of each Ritz value | | |
580 c | by 1 / max(eps23,magnitude of the Ritz value). | | |
581 c %--------------------------------------------------% | |
582 c | |
583 do 35 j = 1, nev0 | |
584 temp = max( eps23, abs(ritz(j)) ) | |
585 bounds(j) = bounds(j)/temp | |
586 35 continue | |
587 c | |
588 c %----------------------------------------------------% | |
589 c | Sort the Ritz values according to the scaled Ritz | | |
590 c | esitmates. This will push all the converged ones | | |
591 c | towards the front of ritzr, ritzi, bounds | | |
592 c | (in the case when NCONV < NEV.) | | |
593 c %----------------------------------------------------% | |
594 c | |
595 wprime = 'LA' | |
596 call dsortr(wprime, .true., nev0, bounds, ritz) | |
597 c | |
598 c %----------------------------------------------% | |
599 c | Scale the Ritz estimate back to its original | | |
600 c | value. | | |
601 c %----------------------------------------------% | |
602 c | |
603 do 40 j = 1, nev0 | |
604 temp = max( eps23, abs(ritz(j)) ) | |
605 bounds(j) = bounds(j)*temp | |
606 40 continue | |
607 c | |
608 c %--------------------------------------------------% | |
609 c | Sort the "converged" Ritz values again so that | | |
610 c | the "threshold" values and their associated Ritz | | |
611 c | estimates appear at the appropriate position in | | |
612 c | ritz and bound. | | |
613 c %--------------------------------------------------% | |
614 c | |
615 if (which .eq. 'BE') then | |
616 c | |
617 c %------------------------------------------------% | |
618 c | Sort the "converged" Ritz values in increasing | | |
619 c | order. The "threshold" values are in the | | |
620 c | middle. | | |
621 c %------------------------------------------------% | |
622 c | |
623 wprime = 'LA' | |
624 call dsortr(wprime, .true., nconv, ritz, bounds) | |
625 c | |
626 else | |
627 c | |
628 c %----------------------------------------------% | |
629 c | In LM, SM, LA, SA case, sort the "converged" | | |
630 c | Ritz values according to WHICH so that the | | |
631 c | "threshold" value appears at the front of | | |
632 c | ritz. | | |
633 c %----------------------------------------------% | |
634 | |
635 call dsortr(which, .true., nconv, ritz, bounds) | |
636 c | |
637 end if | |
638 c | |
639 c %------------------------------------------% | |
640 c | Use h( 1,1 ) as storage to communicate | | |
641 c | rnorm to _seupd if needed | | |
642 c %------------------------------------------% | |
643 c | |
644 h(1,1) = rnorm | |
645 c | |
646 if (msglvl .gt. 1) then | |
647 call dvout (logfil, kplusp, ritz, ndigit, | |
648 & '_saup2: Sorted Ritz values.') | |
649 call dvout (logfil, kplusp, bounds, ndigit, | |
650 & '_saup2: Sorted ritz estimates.') | |
651 end if | |
652 c | |
653 c %------------------------------------% | |
654 c | Max iterations have been exceeded. | | |
655 c %------------------------------------% | |
656 c | |
657 if (iter .gt. mxiter .and. nconv .lt. nev) info = 1 | |
658 c | |
659 c %---------------------% | |
660 c | No shifts to apply. | | |
661 c %---------------------% | |
662 c | |
663 if (np .eq. 0 .and. nconv .lt. nev0) info = 2 | |
664 c | |
665 np = nconv | |
666 go to 1100 | |
667 c | |
668 else if (nconv .lt. nev .and. ishift .eq. 1) then | |
669 c | |
670 c %---------------------------------------------------% | |
671 c | Do not have all the requested eigenvalues yet. | | |
672 c | To prevent possible stagnation, adjust the number | | |
673 c | of Ritz values and the shifts. | | |
674 c %---------------------------------------------------% | |
675 c | |
676 nevbef = nev | |
677 nev = nev + min (nconv, np/2) | |
678 if (nev .eq. 1 .and. kplusp .ge. 6) then | |
679 nev = kplusp / 2 | |
680 else if (nev .eq. 1 .and. kplusp .gt. 2) then | |
681 nev = 2 | |
682 end if | |
683 np = kplusp - nev | |
684 c | |
685 c %---------------------------------------% | |
686 c | If the size of NEV was just increased | | |
687 c | resort the eigenvalues. | | |
688 c %---------------------------------------% | |
689 c | |
690 if (nevbef .lt. nev) | |
691 & call dsgets (ishift, which, nev, np, ritz, bounds, | |
692 & workl) | |
693 c | |
694 end if | |
695 c | |
696 if (msglvl .gt. 0) then | |
697 call ivout (logfil, 1, nconv, ndigit, | |
698 & '_saup2: no. of "converged" Ritz values at this iter.') | |
699 if (msglvl .gt. 1) then | |
700 kp(1) = nev | |
701 kp(2) = np | |
702 call ivout (logfil, 2, kp, ndigit, | |
703 & '_saup2: NEV and NP are') | |
704 call dvout (logfil, nev, ritz(np+1), ndigit, | |
705 & '_saup2: "wanted" Ritz values.') | |
706 call dvout (logfil, nev, bounds(np+1), ndigit, | |
707 & '_saup2: Ritz estimates of the "wanted" values ') | |
708 end if | |
709 end if | |
710 | |
711 c | |
712 if (ishift .eq. 0) then | |
713 c | |
714 c %-----------------------------------------------------% | |
715 c | User specified shifts: reverse communication to | | |
716 c | compute the shifts. They are returned in the first | | |
717 c | NP locations of WORKL. | | |
718 c %-----------------------------------------------------% | |
719 c | |
720 ushift = .true. | |
721 ido = 3 | |
722 go to 9000 | |
723 end if | |
724 c | |
725 50 continue | |
726 c | |
727 c %------------------------------------% | |
728 c | Back from reverse communication; | | |
729 c | User specified shifts are returned | | |
730 c | in WORKL(1:*NP) | | |
731 c %------------------------------------% | |
732 c | |
733 ushift = .false. | |
734 c | |
735 c | |
736 c %---------------------------------------------------------% | |
737 c | Move the NP shifts to the first NP locations of RITZ to | | |
738 c | free up WORKL. This is for the non-exact shift case; | | |
739 c | in the exact shift case, dsgets already handles this. | | |
740 c %---------------------------------------------------------% | |
741 c | |
742 if (ishift .eq. 0) call dcopy (np, workl, 1, ritz, 1) | |
743 c | |
744 if (msglvl .gt. 2) then | |
745 call ivout (logfil, 1, np, ndigit, | |
746 & '_saup2: The number of shifts to apply ') | |
747 call dvout (logfil, np, workl, ndigit, | |
748 & '_saup2: shifts selected') | |
749 if (ishift .eq. 1) then | |
750 call dvout (logfil, np, bounds, ndigit, | |
751 & '_saup2: corresponding Ritz estimates') | |
752 end if | |
753 end if | |
754 c | |
755 c %---------------------------------------------------------% | |
756 c | Apply the NP0 implicit shifts by QR bulge chasing. | | |
757 c | Each shift is applied to the entire tridiagonal matrix. | | |
758 c | The first 2*N locations of WORKD are used as workspace. | | |
759 c | After dsapps is done, we have a Lanczos | | |
760 c | factorization of length NEV. | | |
761 c %---------------------------------------------------------% | |
762 c | |
763 call dsapps (n, nev, np, ritz, v, ldv, h, ldh, resid, q, ldq, | |
764 & workd) | |
765 c | |
766 c %---------------------------------------------% | |
767 c | Compute the B-norm of the updated residual. | | |
768 c | Keep B*RESID in WORKD(1:N) to be used in | | |
769 c | the first step of the next call to dsaitr. | | |
770 c %---------------------------------------------% | |
771 c | |
772 cnorm = .true. | |
773 call arscnd (t2) | |
774 if (bmat .eq. 'G') then | |
775 nbx = nbx + 1 | |
776 call dcopy (n, resid, 1, workd(n+1), 1) | |
777 ipntr(1) = n + 1 | |
778 ipntr(2) = 1 | |
779 ido = 2 | |
780 c | |
781 c %----------------------------------% | |
782 c | Exit in order to compute B*RESID | | |
783 c %----------------------------------% | |
784 c | |
785 go to 9000 | |
786 else if (bmat .eq. 'I') then | |
787 call dcopy (n, resid, 1, workd, 1) | |
788 end if | |
789 c | |
790 100 continue | |
791 c | |
792 c %----------------------------------% | |
793 c | Back from reverse communication; | | |
794 c | WORKD(1:N) := B*RESID | | |
795 c %----------------------------------% | |
796 c | |
797 if (bmat .eq. 'G') then | |
798 call arscnd (t3) | |
799 tmvbx = tmvbx + (t3 - t2) | |
800 end if | |
801 c | |
802 if (bmat .eq. 'G') then | |
803 rnorm = ddot (n, resid, 1, workd, 1) | |
804 rnorm = sqrt(abs(rnorm)) | |
805 else if (bmat .eq. 'I') then | |
806 rnorm = dnrm2(n, resid, 1) | |
807 end if | |
808 cnorm = .false. | |
809 130 continue | |
810 c | |
811 if (msglvl .gt. 2) then | |
812 call dvout (logfil, 1, rnorm, ndigit, | |
813 & '_saup2: B-norm of residual for NEV factorization') | |
814 call dvout (logfil, nev, h(1,2), ndigit, | |
815 & '_saup2: main diagonal of compressed H matrix') | |
816 call dvout (logfil, nev-1, h(2,1), ndigit, | |
817 & '_saup2: subdiagonal of compressed H matrix') | |
818 end if | |
819 c | |
820 go to 1000 | |
821 c | |
822 c %---------------------------------------------------------------% | |
823 c | | | |
824 c | E N D O F M A I N I T E R A T I O N L O O P | | |
825 c | | | |
826 c %---------------------------------------------------------------% | |
827 c | |
828 1100 continue | |
829 c | |
830 mxiter = iter | |
831 nev = nconv | |
832 c | |
833 1200 continue | |
834 ido = 99 | |
835 c | |
836 c %------------% | |
837 c | Error exit | | |
838 c %------------% | |
839 c | |
840 call arscnd (t1) | |
841 tsaup2 = t1 - t0 | |
842 c | |
843 9000 continue | |
844 return | |
845 c | |
846 c %---------------% | |
847 c | End of dsaup2 | | |
848 c %---------------% | |
849 c | |
850 end |