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