6
c Intermediate level interface called by znaupd .
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 )
16
c IDO, BMAT, N, WHICH, NEV, TOL, RESID: same as defined in znaupd .
17
c MODE, ISHIFT, MXITER: see the definition of IPARAM in znaupd .
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
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"
33
c Upon termination of the IRA iteration, NP contains the number
34
c of "converged" wanted Ritz values.
36
c IUPD Integer. (INPUT)
37
c IUPD .EQ. 0: use explicit restart instead implicit update.
38
c IUPD .NE. 0: use implicit update.
40
c V Complex*16 N by (NEV+NP) array. (INPUT/OUTPUT)
41
c The Arnoldi basis vectors are returned in the first NEV
44
c LDV Integer. (INPUT)
45
c Leading dimension of V exactly as declared in the calling
48
c H Complex*16 (NEV+NP) by (NEV+NP) array. (OUTPUT)
49
c H is used to store the generated upper Hessenberg matrix
51
c LDH Integer. (INPUT)
52
c Leading dimension of H exactly as declared in the calling
55
c RITZ Complex*16 array of length NEV+NP. (OUTPUT)
56
c RITZ(1:NEV) contains the computed Ritz values of OP.
58
c BOUNDS Complex*16 array of length NEV+NP. (OUTPUT)
59
c BOUNDS(1:NEV) contain the error bounds corresponding to
60
c the computed Ritz values.
62
c Q Complex*16 (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.
66
c LDQ Integer. (INPUT)
67
c Leading dimension of Q exactly as declared in the calling
70
c WORKL Complex*16 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.
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 -------------------------------------------------------------
87
c WORKD Complex*16 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 ZNAUPD .
93
c RWORK Double precision work array of length NEV+NP ( WORKSPACE)
94
c Private (replicated) array on each PE or array allocated on
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.
115
c-----------------------------------------------------------------------
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),
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.
131
c zgetv0 ARPACK initial vector generation routine.
132
c znaitr ARPACK Arnoldi factorization routine.
133
c znapps ARPACK application of implicit shifts routine.
134
c zneigh ARPACK compute Ritz values and error bounds routine.
135
c zngets ARPACK reorder Ritz values and error bounds routine.
136
c zsortc ARPACK sorting routine.
137
c ivout ARPACK utility routine that prints integers.
138
c second ARPACK utility routine for timing.
139
c zmout ARPACK utility routine that prints matrices
140
c zvout ARPACK utility routine that prints vectors.
141
c dvout ARPACK utility routine that prints vectors.
142
c dlamch LAPACK routine that determines machine constants.
143
c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
144
c zcopy Level 1 BLAS that copies one vector to another .
145
c zdotc Level 1 BLAS that computes the scalar product of two vectors.
146
c zswap Level 1 BLAS that swaps two vectors.
147
c dznrm2 Level 1 BLAS that computes the norm of a vector.
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
158
c\SCCS Information: @(#)
159
c FILE: naup2.F SID: 2.6 DATE OF SID: 06/01/00 RELEASE: 2
166
c-----------------------------------------------------------------------
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 )
173
c %----------------------------------------------------%
174
c | Include files for debugging and timing information |
175
c %----------------------------------------------------%
180
c %------------------%
181
c | Scalar Arguments |
182
c %------------------%
184
character bmat*1, which*2
185
integer ido, info, ishift, iupd, mode, ldh, ldq, ldv, mxiter,
190
c %-----------------%
191
c | Array Arguments |
192
c %-----------------%
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) )
210
parameter (one = (1.0D+0, 0.0D+0) , zero = (0.0D+0, 0.0D+0) ,
217
logical cnorm , getv0, initv , update, ushift
218
integer ierr , iter , kplusp, msglvl, nconv,
219
& nevbef, nev0 , np0 , nptemp, i ,
224
& rnorm , eps23, rtemp
227
save cnorm, getv0, initv , update, ushift,
228
& rnorm, iter , kplusp, msglvl, nconv ,
229
& nevbef, nev0 , np0 , eps23
232
c %-----------------------%
233
c | Local array arguments |
234
c %-----------------------%
238
c %----------------------%
239
c | External Subroutines |
240
c %----------------------%
242
external zcopy , zgetv0 , znaitr , zneigh , zngets , znapps ,
243
& zsortc , zswap , zmout , zvout , ivout, second
245
c %--------------------%
246
c | External functions |
247
c %--------------------%
252
& dznrm2 , dlamch , dlapy2
253
external zdotc , dznrm2 , dlamch , dlapy2
255
c %---------------------%
256
c | Intrinsic Functions |
257
c %---------------------%
259
intrinsic dimag , dble , min, max
261
c %-----------------------%
262
c | Executable Statements |
263
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 %-------------------------------------%
287
c %---------------------------------%
288
c | Get machine dependent constant. |
289
c %---------------------------------%
291
eps23 = dlamch ('Epsilon-Machine')
292
eps23 = eps23**(2.0D+0 / 3.0D+0 )
294
c %---------------------------------------%
295
c | Set flags for computing the first NEV |
296
c | steps of the Arnoldi factorization. |
297
c %---------------------------------------%
304
if (info .ne. 0) then
306
c %--------------------------------------------%
307
c | User provides the initial residual vector. |
308
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 %---------------------------------------------%
325
call zgetv0 (ido, bmat, 1, initv, n, 1, v, ldv, resid, rnorm,
326
& ipntr, workd, info)
328
if (ido .ne. 99) go to 9000
330
if (rnorm .eq. rzero) then
332
c %-----------------------------------------%
333
c | The initial vector is zero. Error exit. |
334
c %-----------------------------------------%
343
c %-----------------------------------%
344
c | Back from reverse communication : |
345
c | continue with update step |
346
c %-----------------------------------%
350
c %-------------------------------------------%
351
c | Back from computing user specified shifts |
352
c %-------------------------------------------%
356
c %-------------------------------------%
357
c | Back from computing residual norm |
358
c | at the end of the current iteration |
359
c %-------------------------------------%
363
c %----------------------------------------------------------%
364
c | Compute the first NEV steps of the Arnoldi factorization |
365
c %----------------------------------------------------------%
367
call znaitr (ido, bmat, n, 0, nev, mode, resid, rnorm, v, ldv,
368
& h, ldh, ipntr, workd, info)
370
if (ido .ne. 99) go to 9000
372
if (info .gt. 0) then
379
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. |
385
c %--------------------------------------------------------------%
391
if (msglvl .gt. 0) then
392
call ivout (logfil, 1, iter, ndigit,
393
& '_naup2: **** Start of major iteration number ****')
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 znapps . |
400
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')
411
c %-----------------------------------------------------------%
412
c | Compute NP additional steps of the Arnoldi factorization. |
413
c %-----------------------------------------------------------%
419
call znaitr (ido, bmat, n, nev, np, mode, resid, rnorm,
420
& v , ldv , h, ldh, ipntr, workd, info)
422
if (ido .ne. 99) go to 9000
424
if (info .gt. 0) then
432
if (msglvl .gt. 1) then
433
call dvout (logfil, 1, rnorm, ndigit,
434
& '_naup2: Corresponding B-norm of the residual')
437
c %--------------------------------------------------------%
438
c | Compute the eigenvalues and corresponding error bounds |
439
c | of the current upper Hessenberg matrix. |
440
c %--------------------------------------------------------%
442
call zneigh (rnorm, kplusp, h, ldh, ritz, bounds,
443
& q, ldq, workl, rwork, ierr)
445
if (ierr .ne. 0) then
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 %---------------------------------------------------%
461
c %--------------------------------------------------%
462
c | Make a copy of Ritz values and the corresponding |
463
c | Ritz estimates obtained from zneigh . |
464
c %--------------------------------------------------%
466
call zcopy (kplusp,ritz,1,workl(kplusp**2+1),1)
467
call zcopy (kplusp,bounds,1,workl(kplusp**2+kplusp+1),1)
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 %---------------------------------------------------%
477
call zngets (ishift, which, nev, np, ritz, bounds)
479
c %------------------------------------------------------------%
480
c | Convergence test: currently we use the following criteria. |
481
c | The relative accuracy of a Ritz value is considered |
484
c | error_bounds(i) .le. tol*max(eps23, magnitude_of_ritz(i)). |
486
c %------------------------------------------------------------%
491
rtemp = max( eps23, dlapy2 ( dble (ritz(np+i)),
492
& dimag (ritz(np+i)) ) )
493
if ( dlapy2 (dble (bounds(np+i)),dimag (bounds(np+i)))
494
& .le. tol*rtemp ) then
499
if (msglvl .gt. 2) then
503
call ivout (logfil, 3, kp, ndigit,
504
& '_naup2: NEV, NP, NCONV are')
505
call zvout (logfil, kplusp, ritz, ndigit,
506
& '_naup2: The eigenvalues of H')
507
call zvout (logfil, kplusp, bounds, ndigit,
508
& '_naup2: Ritz estimates of the current NCV Ritz values')
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 %---------------------------------------------------------%
523
if (bounds(j) .eq. zero) then
529
if ( (nconv .ge. nev0) .or.
530
& (iter .gt. mxiter) .or.
533
if (msglvl .gt. 4) then
534
call zvout (logfil, kplusp, workl(kplusp**2+1), ndigit,
535
& '_naup2: Eigenvalues computed by _neigh:')
536
call zvout (logfil, kplusp, workl(kplusp**2+kplusp+1),
538
& '_naup2: Ritz estimates computed by _neigh:')
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 %------------------------------------------------%
548
c %------------------------------------------%
549
c | Use h( 3,1 ) as storage to communicate |
550
c | rnorm to zneupd if needed |
551
c %------------------------------------------%
553
h(3,1) = dcmplx (rnorm,rzero)
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 %----------------------------------------------%
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'
569
call zsortc (wprime, .true., kplusp, ritz, bounds)
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 %--------------------------------------------------%
577
rtemp = max( eps23, dlapy2 ( dble (ritz(j)),
578
& dimag (ritz(j)) ) )
579
bounds(j) = bounds(j)/rtemp
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 %---------------------------------------------------%
590
call zsortc (wprime, .true., nev0, bounds, ritz)
592
c %----------------------------------------------%
593
c | Scale the Ritz estimate back to its original |
595
c %----------------------------------------------%
598
rtemp = max( eps23, dlapy2 ( dble (ritz(j)),
599
& dimag (ritz(j)) ) )
600
bounds(j) = bounds(j)*rtemp
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 %-----------------------------------------------%
609
call zsortc (which, .true., nconv, ritz, bounds)
611
if (msglvl .gt. 1) then
612
call zvout (logfil, kplusp, ritz, ndigit,
613
& '_naup2: Sorted eigenvalues')
614
call zvout (logfil, kplusp, bounds, ndigit,
615
& '_naup2: Sorted ritz estimates.')
618
c %------------------------------------%
619
c | Max iterations have been exceeded. |
620
c %------------------------------------%
622
if (iter .gt. mxiter .and. nconv .lt. nev0) info = 1
624
c %---------------------%
625
c | No shifts to apply. |
626
c %---------------------%
628
if (np .eq. 0 .and. nconv .lt. nev0) info = 2
633
else if ( (nconv .lt. nev0) .and. (ishift .eq. 1) ) then
635
c %-------------------------------------------------%
636
c | Do not have all the requested eigenvalues yet. |
637
c | To prevent possible stagnation, adjust the size |
639
c %-------------------------------------------------%
642
nev = nev + min(nconv, np/2)
643
if (nev .eq. 1 .and. kplusp .ge. 6) then
645
else if (nev .eq. 1 .and. kplusp .gt. 3) then
650
c %---------------------------------------%
651
c | If the size of NEV was just increased |
652
c | resort the eigenvalues. |
653
c %---------------------------------------%
656
& call zngets (ishift, which, nev, np, ritz, bounds)
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
666
call ivout (logfil, 2, kp, ndigit,
667
& '_naup2: NEV and NP are')
668
call zvout (logfil, nev, ritz(np+1), ndigit,
669
& '_naup2: "wanted" Ritz values ')
670
call zvout (logfil, nev, bounds(np+1), ndigit,
671
& '_naup2: Ritz estimates of the "wanted" values ')
675
if (ishift .eq. 0) then
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 %-------------------------------------------------------%
689
if ( ishift .ne. 1 ) then
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 %----------------------------------%
697
call zcopy (np, workl, 1, ritz, 1)
700
if (msglvl .gt. 2) then
701
call ivout (logfil, 1, np, ndigit,
702
& '_naup2: The number of shifts to apply ')
703
call zvout (logfil, np, ritz, ndigit,
704
& '_naup2: values of the shifts')
706
& call zvout (logfil, np, bounds, ndigit,
707
& '_naup2: Ritz estimates of the shifts')
710
c %---------------------------------------------------------%
711
c | Apply the NP implicit shifts by QR bulge chasing. |
712
c | Each shift is applied to the whole upper Hessenberg |
714
c | The first 2*N locations of WORKD are used as workspace. |
715
c %---------------------------------------------------------%
717
call znapps (n, nev, np, ritz, v, ldv,
718
& h, ldh, resid, q, ldq, workl, workd)
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 znaitr . |
724
c %---------------------------------------------%
728
if (bmat .eq. 'G') then
730
call zcopy (n, resid, 1, workd(n+1), 1)
735
c %----------------------------------%
736
c | Exit in order to compute B*RESID |
737
c %----------------------------------%
740
else if (bmat .eq. 'I') then
741
call zcopy (n, resid, 1, workd, 1)
746
c %----------------------------------%
747
c | Back from reverse communication; |
748
c | WORKD(1:N) := B*RESID |
749
c %----------------------------------%
751
if (bmat .eq. 'G') then
753
tmvbx = tmvbx + (t3 - t2)
756
if (bmat .eq. 'G') then
757
cmpnorm = zdotc (n, resid, 1, workd, 1)
758
rnorm = sqrt(dlapy2 (dble (cmpnorm),dimag (cmpnorm)))
759
else if (bmat .eq. 'I') then
760
rnorm = dznrm2 (n, resid, 1)
764
if (msglvl .gt. 2) then
765
call dvout (logfil, 1, rnorm, ndigit,
766
& '_naup2: B-norm of residual for compressed factorization')
767
call zmout (logfil, nev, nev, h, ldh, ndigit,
768
& '_naup2: Compressed upper Hessenberg matrix H')
773
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 |
777
c %---------------------------------------------------------------%