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