Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/parpack/pcnaupd.f
5209 views
1
c\BeginDoc
2
c
3
c\Name: pcnaupd
4
c
5
c Message Passing Layer: MPI
6
c
7
c\Description:
8
c Reverse communication interface for the Implicitly Restarted Arnoldi
9
c iteration. This is intended to be used to find a few eigenpairs of a
10
c complex linear operator OP with respect to a semi-inner product defined
11
c by a hermitian positive semi-definite real matrix B. B may be the identity
12
c matrix. NOTE: if both OP and B are real, then ssaupd or snaupd should
13
c be used.
14
c
15
c
16
c The computed approximate eigenvalues are called Ritz values and
17
c the corresponding approximate eigenvectors are called Ritz vectors.
18
c
19
c pcnaupd is usually called iteratively to solve one of the
20
c following problems:
21
c
22
c Mode 1: A*x = lambda*x.
23
c ===> OP = A and B = I.
24
c
25
c Mode 2: A*x = lambda*M*x, M symmetric positive definite
26
c ===> OP = inv[M]*A and B = M.
27
c ===> (If M can be factored see remark 3 below)
28
c
29
c Mode 3: A*x = lambda*M*x, M symmetric semi-definite
30
c ===> OP = inv[A - sigma*M]*M and B = M.
31
c ===> shift-and-invert mode
32
c If OP*x = amu*x, then lambda = sigma + 1/amu.
33
c
34
c
35
c NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v
36
c should be accomplished either by a direct method
37
c using a sparse matrix factorization and solving
38
c
39
c [A - sigma*M]*w = v or M*w = v,
40
c
41
c or through an iterative method for solving these
42
c systems. If an iterative method is used, the
43
c convergence test must be more stringent than
44
c the accuracy requirements for the eigenvalue
45
c approximations.
46
c
47
c\Usage:
48
c call pcnaupd
49
c ( COMM, IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
50
c IPNTR, WORKD, WORKL, LWORKL, RWORK, INFO )
51
c
52
c\Arguments
53
c COMM MPI Communicator for the processor grid. (INPUT)
54
c
55
c IDO Integer. (INPUT/OUTPUT)
56
c Reverse communication flag. IDO must be zero on the first
57
c call to pcnaupd. IDO will be set internally to
58
c indicate the type of operation to be performed. Control is
59
c then given back to the calling routine which has the
60
c responsibility to carry out the requested operation and call
61
c pcnaupd with the result. The operand is given in
62
c WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)).
63
c -------------------------------------------------------------
64
c IDO = 0: first call to the reverse communication interface
65
c IDO = -1: compute Y = OP * X where
66
c IPNTR(1) is the pointer into WORKD for X,
67
c IPNTR(2) is the pointer into WORKD for Y.
68
c This is for the initialization phase to force the
69
c starting vector into the range of OP.
70
c IDO = 1: compute Y = OP * X where
71
c IPNTR(1) is the pointer into WORKD for X,
72
c IPNTR(2) is the pointer into WORKD for Y.
73
c In mode 3, the vector B * X is already
74
c available in WORKD(ipntr(3)). It does not
75
c need to be recomputed in forming OP * X.
76
c IDO = 2: compute Y = M * X where
77
c IPNTR(1) is the pointer into WORKD for X,
78
c IPNTR(2) is the pointer into WORKD for Y.
79
c IDO = 3: compute and return the shifts in the first
80
c NP locations of WORKL.
81
c IDO = 99: done
82
c -------------------------------------------------------------
83
c After the initialization phase, when the routine is used in
84
c the "shift-and-invert" mode, the vector M * X is already
85
c available and does not need to be recomputed in forming OP*X.
86
c
87
c BMAT Character*1. (INPUT)
88
c BMAT specifies the type of the matrix B that defines the
89
c semi-inner product for the operator OP.
90
c BMAT = 'I' -> standard eigenvalue problem A*x = lambda*x
91
c BMAT = 'G' -> generalized eigenvalue problem A*x = lambda*M*x
92
c
93
c N Integer. (INPUT)
94
c Dimension of the eigenproblem.
95
c
96
c WHICH Character*2. (INPUT)
97
c 'LM' -> want the NEV eigenvalues of largest magnitude.
98
c 'SM' -> want the NEV eigenvalues of smallest magnitude.
99
c 'LR' -> want the NEV eigenvalues of largest real part.
100
c 'SR' -> want the NEV eigenvalues of smallest real part.
101
c 'LI' -> want the NEV eigenvalues of largest imaginary part.
102
c 'SI' -> want the NEV eigenvalues of smallest imaginary part.
103
c
104
c NEV Integer. (INPUT)
105
c Number of eigenvalues of OP to be computed. 0 < NEV < N-1.
106
c
107
c TOL Real scalar. (INPUT)
108
c Stopping criteria: the relative accuracy of the Ritz value
109
c is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I))
110
c where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex.
111
c DEFAULT = pslamch(comm, 'EPS') (machine precision as computed
112
c by the ScaLAPACK auxiliary subroutine pslamch).
113
c
114
c RESID Complex array of length N. (INPUT/OUTPUT)
115
c On INPUT:
116
c If INFO .EQ. 0, a random initial residual vector is used.
117
c If INFO .NE. 0, RESID contains the initial residual vector,
118
c possibly from a previous run.
119
c On OUTPUT:
120
c RESID contains the final residual vector.
121
c
122
c NCV Integer. (INPUT)
123
c Number of columns of the matrix V. NCV must satisfy the two
124
c inequalities 1 <= NCV-NEV and NCV <= N.
125
c This will indicate how many Arnoldi vectors are generated
126
c at each iteration. After the startup phase in which NEV
127
c Arnoldi vectors are generated, the algorithm generates
128
c approximately NCV-NEV Arnoldi vectors at each subsequent update
129
c iteration. Most of the cost in generating each Arnoldi vector is
130
c in the matrix-vector operation OP*x. (See remark 4 below)
131
c
132
c V Complex array N by NCV. (OUTPUT)
133
c Contains the final set of Arnoldi basis vectors.
134
c
135
c LDV Integer. (INPUT)
136
c Leading dimension of V exactly as declared in the calling program.
137
c
138
c IPARAM Integer array of length 11. (INPUT/OUTPUT)
139
c IPARAM(1) = ISHIFT: method for selecting the implicit shifts.
140
c The shifts selected at each iteration are used to filter out
141
c the components of the unwanted eigenvector.
142
c -------------------------------------------------------------
143
c ISHIFT = 0: the shifts are to be provided by the user via
144
c reverse communication. The NCV eigenvalues of
145
c the Hessenberg matrix H are returned in the part
146
c of WORKL array corresponding to RITZ.
147
c ISHIFT = 1: exact shifts with respect to the current
148
c Hessenberg matrix H. This is equivalent to
149
c restarting the iteration from the beginning
150
c after updating the starting vector with a linear
151
c combination of Ritz vectors associated with the
152
c "wanted" eigenvalues.
153
c ISHIFT = 2: other choice of internal shift to be defined.
154
c -------------------------------------------------------------
155
c
156
c IPARAM(2) = No longer referenced
157
c
158
c IPARAM(3) = MXITER
159
c On INPUT: maximum number of Arnoldi update iterations allowed.
160
c On OUTPUT: actual number of Arnoldi update iterations taken.
161
c
162
c IPARAM(4) = NB: blocksize to be used in the recurrence.
163
c The code currently works only for NB = 1.
164
c
165
c IPARAM(5) = NCONV: number of "converged" Ritz values.
166
c This represents the number of Ritz values that satisfy
167
c the convergence criterion.
168
c
169
c IPARAM(6) = IUPD
170
c No longer referenced. Implicit restarting is ALWAYS used.
171
c
172
c IPARAM(7) = MODE
173
c On INPUT determines what type of eigenproblem is being solved.
174
c Must be 1,2,3; See under \Description of pcnaupd for the
175
c four modes available.
176
c
177
c IPARAM(8) = NP
178
c When ido = 3 and the user provides shifts through reverse
179
c communication (IPARAM(1)=0), _naupd returns NP, the number
180
c of shifts the user is to provide. 0 < NP < NCV-NEV.
181
c
182
c IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,
183
c OUTPUT: NUMOP = total number of OP*x operations,
184
c NUMOPB = total number of B*x operations if BMAT='G',
185
c NUMREO = total number of steps of re-orthogonalization.
186
c
187
c IPNTR Integer array of length 14. (OUTPUT)
188
c Pointer to mark the starting locations in the WORKD and WORKL
189
c arrays for matrices/vectors used by the Arnoldi iteration.
190
c -------------------------------------------------------------
191
c IPNTR(1): pointer to the current operand vector X in WORKD.
192
c IPNTR(2): pointer to the current result vector Y in WORKD.
193
c IPNTR(3): pointer to the vector B * X in WORKD when used in
194
c the shift-and-invert mode.
195
c IPNTR(4): pointer to the next available location in WORKL
196
c that is untouched by the program.
197
c IPNTR(5): pointer to the NCV by NCV upper Hessenberg
198
c matrix H in WORKL.
199
c IPNTR(6): pointer to the ritz value array RITZ
200
c IPNTR(7): pointer to the (projected) ritz vector array Q
201
c IPNTR(8): pointer to the error BOUNDS array in WORKL.
202
c IPNTR(14): pointer to the NP shifts in WORKL. See Remark 5 below.
203
c
204
c Note: IPNTR(9:13) is only referenced by pcneupd. See Remark 2 below.
205
c
206
c IPNTR(9): pointer to the NCV RITZ values of the
207
c original system.
208
c IPNTR(10): Not Used
209
c IPNTR(11): pointer to the NCV corresponding error bounds.
210
c IPNTR(12): pointer to the NCV by NCV upper triangular
211
c Schur matrix for H.
212
c IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors
213
c of the upper Hessenberg matrix H. Only referenced by
214
c cneupd if RVEC = .TRUE. See Remark 2 below.
215
c -------------------------------------------------------------
216
c
217
c WORKD Complex work array of length 3*N. (REVERSE COMMUNICATION)
218
c Distributed array to be used in the basic Arnoldi iteration
219
c for reverse communication. The user should not use WORKD
220
c as temporary workspace during the iteration !!!!!!!!!!
221
c See Data Distribution Note below.
222
c
223
c WORKL Complex work array of length LWORKL. (OUTPUT/WORKSPACE)
224
c Private (replicated) array on each PE or array allocated on
225
c the front end. See Data Distribution Note below.
226
c
227
c LWORKL Integer. (INPUT)
228
c LWORKL must be at least 3*NCV**2 + 5*NCV.
229
c
230
c RWORK Real work array of length NCV (WORKSPACE)
231
c Private (replicated) array on each PE or array allocated on
232
c the front end.
233
c
234
c
235
c INFO Integer. (INPUT/OUTPUT)
236
c If INFO .EQ. 0, a randomly initial residual vector is used.
237
c If INFO .NE. 0, RESID contains the initial residual vector,
238
c possibly from a previous run.
239
c Error flag on output.
240
c = 0: Normal exit.
241
c = 1: Maximum number of iterations taken.
242
c All possible eigenvalues of OP has been found. IPARAM(5)
243
c returns the number of wanted converged Ritz values.
244
c = 2: No longer an informational error. Deprecated starting
245
c with release 2 of ARPACK.
246
c = 3: No shifts could be applied during a cycle of the
247
c Implicitly restarted Arnoldi iteration. One possibility
248
c is to increase the size of NCV relative to NEV.
249
c See remark 4 below.
250
c = -1: N must be positive.
251
c = -2: NEV must be positive.
252
c = -3: NCV-NEV >= 2 and less than or equal to N.
253
c = -4: The maximum number of Arnoldi update iteration
254
c must be greater than zero.
255
c = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
256
c = -6: BMAT must be one of 'I' or 'G'.
257
c = -7: Length of private work array is not sufficient.
258
c = -8: Error return from LAPACK eigenvalue calculation;
259
c = -9: Starting vector is zero.
260
c = -10: IPARAM(7) must be 1,2,3.
261
c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
262
c = -12: IPARAM(1) must be equal to 0 or 1.
263
c = -9999: Could not build an Arnoldi factorization.
264
c User input error highly likely. Please
265
c check actual array dimensions and layout.
266
c IPARAM(5) returns the size of the current Arnoldi
267
c factorization.
268
c
269
c\Remarks
270
c 1. The computed Ritz values are approximate eigenvalues of OP. The
271
c selection of WHICH should be made with this in mind when using
272
c Mode = 3. When operating in Mode = 3 setting WHICH = 'LM' will
273
c compute the NEV eigenvalues of the original problem that are
274
c closest to the shift SIGMA . After convergence, approximate eigenvalues
275
c of the original problem may be obtained with the ARPACK subroutine pcneupd.
276
c
277
c 2. If a basis for the invariant subspace corresponding to the converged Ritz
278
c values is needed, the user must call pcneupd immediately following
279
c completion of pcnaupd. This is new starting with release 2 of ARPACK.
280
c
281
c 3. If M can be factored into a Cholesky factorization M = LL`
282
c then Mode = 2 should not be selected. Instead one should use
283
c Mode = 1 with OP = inv(L)*A*inv(L`). Appropriate triangular
284
c linear systems should be solved with L and L` rather
285
c than computing inverses. After convergence, an approximate
286
c eigenvector z of the original problem is recovered by solving
287
c L`z = x where x is a Ritz vector of OP.
288
c
289
c 4. At present there is no a-priori analysis to guide the selection
290
c of NCV relative to NEV. The only formal requrement is that NCV > NEV + 1.
291
c However, it is recommended that NCV .ge. 2*NEV. If many problems of
292
c the same type are to be solved, one should experiment with increasing
293
c NCV while keeping NEV fixed for a given test problem. This will
294
c usually decrease the required number of OP*x operations but it
295
c also increases the work and storage required to maintain the orthogonal
296
c basis vectors. The optimal "cross-over" with respect to CPU time
297
c is problem dependent and must be determined empirically.
298
c See Chapter 8 of Reference 2 for further information.
299
c
300
c 5. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the
301
c NP = IPARAM(8) complex shifts in locations
302
c WORKL(IPNTR(14)), WORKL(IPNTR(14)+1), ... , WORKL(IPNTR(14)+NP).
303
c Eigenvalues of the current upper Hessenberg matrix are located in
304
c WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1). They are ordered
305
c according to the order defined by WHICH. The associated Ritz estimates
306
c are located in WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... ,
307
c WORKL(IPNTR(8)+NCV-1).
308
c
309
c-----------------------------------------------------------------------
310
c
311
c\Data Distribution Note:
312
c
313
c Fortran-D syntax:
314
c ================
315
c Complex resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
316
c decompose d1(n), d2(n,ncv)
317
c align resid(i) with d1(i)
318
c align v(i,j) with d2(i,j)
319
c align workd(i) with d1(i) range (1:n)
320
c align workd(i) with d1(i-n) range (n+1:2*n)
321
c align workd(i) with d1(i-2*n) range (2*n+1:3*n)
322
c distribute d1(block), d2(block,:)
323
c replicated workl(lworkl)
324
c
325
c Cray MPP syntax:
326
c ===============
327
c Complex resid(n), v(ldv,ncv), workd(n,3), workl(lworkl)
328
c shared resid(block), v(block,:), workd(block,:)
329
c replicated workl(lworkl)
330
c
331
c CM2/CM5 syntax:
332
c ==============
333
c
334
c-----------------------------------------------------------------------
335
c
336
c include 'ex-nonsym.doc'
337
c
338
c-----------------------------------------------------------------------
339
c
340
c\BeginLib
341
c
342
c\Local variables:
343
c xxxxxx Complex
344
c
345
c\References:
346
c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
347
c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
348
c pp 357-385.
349
c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
350
c Restarted Arnoldi Iteration", Rice University Technical Report
351
c TR95-13, Department of Computational and Applied Mathematics.
352
c 3. B.N. Parlett & Y. Saad, "_Complex_ Shift and Invert Strategies for
353
c _Real_ Matrices", Linear Algebra and its Applications, vol 88/89,
354
c pp 575-595, (1987).
355
c
356
c\Routines called:
357
c pcnaup2 Parallel ARPACK routine that implements the Implicitly Restarted
358
c Arnoldi Iteration.
359
c cstatn ARPACK routine that initializes the timing variables.
360
c pivout Parallel ARPACK utility routine that prints integers.
361
c pcvout Parallel ARPACK utility routine that prints vectors.
362
c second ARPACK utility routine for timing.
363
c pslamch ScaLAPACK routine that determines machine constants.
364
c
365
c\Author
366
c Danny Sorensen Phuong Vu
367
c Richard Lehoucq CRPC / Rice University
368
c Dept. of Computational & Houston, Texas
369
c Applied Mathematics
370
c Rice University
371
c Houston, Texas
372
c
373
c\Parallel Modifications
374
c Kristi Maschhoff
375
c
376
c\Revision history:
377
c Starting Point: Complex Serial Code FILE: naupd.F SID: 2.2
378
c
379
c\SCCS Information:
380
c FILE: naupd.F SID: 1.7 DATE OF SID: 04/10/01
381
c
382
c\Remarks
383
c
384
c\EndLib
385
c
386
c-----------------------------------------------------------------------
387
c
388
subroutine pcnaupd
389
& ( comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,
390
& iparam, ipntr, workd, workl, lworkl, rwork, info )
391
c
392
include 'mpif.h'
393
c
394
c %------------------%
395
c | MPI Variables |
396
c %------------------%
397
c
398
integer comm, myid
399
c
400
c %----------------------------------------------------%
401
c | Include files for debugging and timing information |
402
c %----------------------------------------------------%
403
c
404
include 'debug.h'
405
include 'stat.h'
406
c
407
c %------------------%
408
c | Scalar Arguments |
409
c %------------------%
410
c
411
character bmat*1, which*2
412
integer ido, info, ldv, lworkl, n, ncv, nev
413
Real
414
& tol
415
c
416
c %-----------------%
417
c | Array Arguments |
418
c %-----------------%
419
c
420
integer iparam(11), ipntr(14)
421
Complex
422
& resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
423
Real
424
& rwork(ncv)
425
c
426
c %------------%
427
c | Parameters |
428
c %------------%
429
c
430
Complex
431
& one, zero
432
parameter (one = (1.0, 0.0) , zero = (0.0, 0.0) )
433
c
434
c %---------------%
435
c | Local Scalars |
436
c %---------------%
437
c
438
integer bounds, ierr, ih, iq, ishift, iupd, iw,
439
& ldh, ldq, levec, mode, msglvl, mxiter, nb,
440
& nev0, next, np, ritz, j
441
save bounds, ih, iq, ishift, iupd, iw,
442
& ldh, ldq, levec, mode, msglvl, mxiter, nb,
443
& nev0, next, np, ritz
444
c
445
c %----------------------%
446
c | External Subroutines |
447
c %----------------------%
448
c
449
external pcnaup2, pcvout, pivout, second, cstatn
450
c
451
c %--------------------%
452
c | External Functions |
453
c %--------------------%
454
c
455
Real
456
& pslamch10
457
external pslamch10
458
c
459
c %-----------------------%
460
c | Executable Statements |
461
c %-----------------------%
462
c
463
if (ido .eq. 0) then
464
c
465
c %-------------------------------%
466
c | Initialize timing statistics |
467
c | & message level for debugging |
468
c %-------------------------------%
469
c
470
call cstatn
471
call second (t0)
472
msglvl = mcaupd
473
c
474
c %----------------%
475
c | Error checking |
476
c %----------------%
477
c
478
ierr = 0
479
ishift = iparam(1)
480
c levec = iparam(2)
481
mxiter = iparam(3)
482
c nb = iparam(4)
483
nb = 1
484
c
485
c %--------------------------------------------%
486
c | Revision 2 performs only implicit restart. |
487
c %--------------------------------------------%
488
c
489
iupd = 1
490
mode = iparam(7)
491
c
492
if (n .le. 0) then
493
ierr = -1
494
else if (nev .le. 0) then
495
ierr = -2
496
else if (ncv .le. nev) then
497
ierr = -3
498
else if (mxiter .le. 0) then
499
ierr = -4
500
else if (which .ne. 'LM' .and.
501
& which .ne. 'SM' .and.
502
& which .ne. 'LR' .and.
503
& which .ne. 'SR' .and.
504
& which .ne. 'LI' .and.
505
& which .ne. 'SI') then
506
ierr = -5
507
else if (bmat .ne. 'I' .and. bmat .ne. 'G') then
508
ierr = -6
509
else if (lworkl .lt. 3*ncv**2 + 5*ncv) then
510
ierr = -7
511
else if (mode .lt. 1 .or. mode .gt. 3) then
512
ierr = -10
513
else if (mode .eq. 1 .and. bmat .eq. 'G') then
514
ierr = -11
515
end if
516
c
517
c %------------%
518
c | Error Exit |
519
c %------------%
520
c
521
if (ierr .ne. 0) then
522
info = ierr
523
ido = 99
524
go to 9000
525
end if
526
c
527
c %------------------------%
528
c | Set default parameters |
529
c %------------------------%
530
c
531
if (nb .le. 0) nb = 1
532
if (tol .le. 0.0 ) tol = pslamch10(comm, 'EpsMach')
533
if (ishift .ne. 0 .and.
534
& ishift .ne. 1 .and.
535
& ishift .ne. 2) ishift = 1
536
c
537
c %----------------------------------------------%
538
c | NP is the number of additional steps to |
539
c | extend the length NEV Lanczos factorization. |
540
c | NEV0 is the local variable designating the |
541
c | size of the invariant subspace desired. |
542
c %----------------------------------------------%
543
c
544
np = ncv - nev
545
nev0 = nev
546
c
547
c %-----------------------------%
548
c | Zero out internal workspace |
549
c %-----------------------------%
550
c
551
do 10 j = 1, 3*ncv**2 + 5*ncv
552
workl(j) = zero
553
10 continue
554
c
555
c %-------------------------------------------------------------%
556
c | Pointer into WORKL for address of H, RITZ, BOUNDS, Q |
557
c | etc... and the remaining workspace. |
558
c | Also update pointer to be used on output. |
559
c | Memory is laid out as follows: |
560
c | workl(1:ncv*ncv) := generated Hessenberg matrix |
561
c | workl(ncv*ncv+1:ncv*ncv+ncv) := the ritz values |
562
c | workl(ncv*ncv+ncv+1:ncv*ncv+2*ncv) := error bounds |
563
c | workl(ncv*ncv+2*ncv+1:2*ncv*ncv+2*ncv) := rotation matrix Q |
564
c | workl(2*ncv*ncv+2*ncv+1:3*ncv*ncv+5*ncv) := workspace |
565
c | The final workspace is needed by subroutine pcneigh called |
566
c | by pcnaup2. Subroutine pcneigh calls LAPACK routines for |
567
c | calculating eigenvalues and the last row of the eigenvector |
568
c | matrix. |
569
c %-------------------------------------------------------------%
570
c
571
ldh = ncv
572
ldq = ncv
573
ih = 1
574
ritz = ih + ldh*ncv
575
bounds = ritz + ncv
576
iq = bounds + ncv
577
iw = iq + ldq*ncv
578
next = iw + ncv**2 + 3*ncv
579
c
580
ipntr(4) = next
581
ipntr(5) = ih
582
ipntr(6) = ritz
583
ipntr(7) = iq
584
ipntr(8) = bounds
585
ipntr(14) = iw
586
end if
587
c
588
c %-------------------------------------------------------%
589
c | Carry out the Implicitly restarted Arnoldi Iteration. |
590
c %-------------------------------------------------------%
591
c
592
call pcnaup2
593
& ( comm, ido, bmat, n, which, nev0, np, tol, resid, mode, iupd,
594
& ishift, mxiter, v, ldv, workl(ih), ldh, workl(ritz),
595
& workl(bounds), workl(iq), ldq, workl(iw),
596
& ipntr, workd, rwork, info )
597
c
598
c %--------------------------------------------------%
599
c | ido .ne. 99 implies use of reverse communication |
600
c | to compute operations involving OP. |
601
c %--------------------------------------------------%
602
c
603
if (ido .eq. 3) iparam(8) = np
604
if (ido .ne. 99) go to 9000
605
c
606
iparam(3) = mxiter
607
iparam(5) = np
608
iparam(9) = nopx
609
iparam(10) = nbx
610
iparam(11) = nrorth
611
c
612
c %------------------------------------%
613
c | Exit if there was an informational |
614
c | error within pcnaup2. |
615
c %------------------------------------%
616
c
617
if (info .lt. 0) go to 9000
618
if (info .eq. 2) info = 3
619
c
620
if (msglvl .gt. 0) then
621
call pivout (comm, logfil, 1, [mxiter], ndigit,
622
& '_naupd: Number of update iterations taken')
623
call pivout (comm, logfil, 1, [np], ndigit,
624
& '_naupd: Number of wanted "converged" Ritz values')
625
call pcvout (comm, logfil, np, workl(ritz), ndigit,
626
& '_naupd: The final Ritz values')
627
call pcvout (comm, logfil, np, workl(bounds), ndigit,
628
& '_naupd: Associated Ritz estimates')
629
end if
630
c
631
call second (t1)
632
tcaupd = t1 - t0
633
c
634
if (msglvl .gt. 0) then
635
call MPI_COMM_RANK( comm, myid, ierr )
636
if ( myid .eq. 0 ) then
637
c
638
c %--------------------------------------------------------%
639
c | Version Number & Version Date are defined in version.h |
640
c %--------------------------------------------------------%
641
c
642
write (6,1000)
643
write (6,1100) mxiter, nopx, nbx, nrorth, nitref, nrstrt,
644
& tmvopx, tmvbx, tcaupd, tcaup2, tcaitr, titref,
645
& tgetv0, tceigh, tcgets, tcapps, tcconv, trvec
646
1000 format (//,
647
& 5x, '=============================================',/
648
& 5x, '= Complex implicit Arnoldi update code =',/
649
& 5x, '= Version Number: ', ' 2.1' , 21x, ' =',/
650
& 5x, '= Version Date: ', ' 3/19/97' , 16x, ' =',/
651
& 5x, '=============================================',/
652
& 5x, '= Summary of timing statistics =',/
653
& 5x, '=============================================',//)
654
1100 format (
655
& 5x, 'Total number update iterations = ', i5,/
656
& 5x, 'Total number of OP*x operations = ', i5,/
657
& 5x, 'Total number of B*x operations = ', i5,/
658
& 5x, 'Total number of reorthogonalization steps = ', i5,/
659
& 5x, 'Total number of iterative refinement steps = ', i5,/
660
& 5x, 'Total number of restart steps = ', i5,/
661
& 5x, 'Total time in user OP*x operation = ', f12.6,/
662
& 5x, 'Total time in user B*x operation = ', f12.6,/
663
& 5x, 'Total time in Arnoldi update routine = ', f12.6,/
664
& 5x, 'Total time in p_naup2 routine = ', f12.6,/
665
& 5x, 'Total time in basic Arnoldi iteration loop = ', f12.6,/
666
& 5x, 'Total time in reorthogonalization phase = ', f12.6,/
667
& 5x, 'Total time in (re)start vector generation = ', f12.6,/
668
& 5x, 'Total time in Hessenberg eig. subproblem = ', f12.6,/
669
& 5x, 'Total time in getting the shifts = ', f12.6,/
670
& 5x, 'Total time in applying the shifts = ', f12.6,/
671
& 5x, 'Total time in convergence testing = ', f12.6,/
672
& 5x, 'Total time in computing final Ritz vectors = ', f12.6/)
673
end if
674
end if
675
c
676
9000 continue
677
c
678
return
679
c
680
c %----------------%
681
c | End of pcnaupd |
682
c %----------------%
683
c
684
end
685
686