1
c-----------------------------------------------------------------------
8
c Reverse communication interface for the Implicitly Restarted Arnoldi
9
c Iteration. For symmetric problems this reduces to a variant of the Lanczos
10
c method. This method has been designed to compute approximations to a
11
c few eigenpairs of a linear operator OP that is real and symmetric
12
c with respect to a real positive semi-definite symmetric matrix B,
17
c Another way to express this condition is
19
c < x,OPy > = < OPx,y > where < z,w > = z`Bw .
21
c In the standard eigenproblem B is the identity matrix.
22
c ( A` denotes transpose of A)
24
c The computed approximate eigenvalues are called Ritz values and
25
c the corresponding approximate eigenvectors are called Ritz vectors.
27
c dsaupd is usually called iteratively to solve one of the
30
c Mode 1: A*x = lambda*x, A symmetric
31
c ===> OP = A and B = I.
33
c Mode 2: A*x = lambda*M*x, A symmetric, M symmetric positive definite
34
c ===> OP = inv[M]*A and B = M.
35
c ===> (If M can be factored see remark 3 below)
37
c Mode 3: K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
38
c ===> OP = (inv[K - sigma*M])*M and B = M.
39
c ===> Shift-and-Invert mode
41
c Mode 4: K*x = lambda*KG*x, K symmetric positive semi-definite,
42
c KG symmetric indefinite
43
c ===> OP = (inv[K - sigma*KG])*K and B = K.
46
c Mode 5: A*x = lambda*M*x, A symmetric, M symmetric positive semi-definite
47
c ===> OP = inv[A - sigma*M]*[A + sigma*M] and B = M.
48
c ===> Cayley transformed mode
50
c NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v
51
c should be accomplished either by a direct method
52
c using a sparse matrix factorization and solving
54
c [A - sigma*M]*w = v or M*w = v,
56
c or through an iterative method for solving these
57
c systems. If an iterative method is used, the
58
c convergence test must be more stringent than
59
c the accuracy requirements for the eigenvalue
64
c ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
65
c IPNTR, WORKD, WORKL, LWORKL, INFO )
68
c IDO Integer. (INPUT/OUTPUT)
69
c Reverse communication flag. IDO must be zero on the first
70
c call to dsaupd . IDO will be set internally to
71
c indicate the type of operation to be performed. Control is
72
c then given back to the calling routine which has the
73
c responsibility to carry out the requested operation and call
74
c dsaupd with the result. The operand is given in
75
c WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)).
76
c (If Mode = 2 see remark 5 below)
77
c -------------------------------------------------------------
78
c IDO = 0: first call to the reverse communication interface
79
c IDO = -1: compute Y = OP * X where
80
c IPNTR(1) is the pointer into WORKD for X,
81
c IPNTR(2) is the pointer into WORKD for Y.
82
c This is for the initialization phase to force the
83
c starting vector into the range of OP.
84
c IDO = 1: compute Y = OP * X where
85
c IPNTR(1) is the pointer into WORKD for X,
86
c IPNTR(2) is the pointer into WORKD for Y.
87
c In mode 3,4 and 5, the vector B * X is already
88
c available in WORKD(ipntr(3)). It does not
89
c need to be recomputed in forming OP * X.
90
c IDO = 2: compute Y = B * X where
91
c IPNTR(1) is the pointer into WORKD for X,
92
c IPNTR(2) is the pointer into WORKD for Y.
93
c IDO = 3: compute the IPARAM(8) shifts where
94
c IPNTR(11) is the pointer into WORKL for
95
c placing the shifts. See remark 6 below.
97
c -------------------------------------------------------------
99
c BMAT Character*1. (INPUT)
100
c BMAT specifies the type of the matrix B that defines the
101
c semi-inner product for the operator OP.
102
c B = 'I' -> standard eigenvalue problem A*x = lambda*x
103
c B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
106
c Dimension of the eigenproblem.
108
c WHICH Character*2. (INPUT)
109
c Specify which of the Ritz values of OP to compute.
111
c 'LA' - compute the NEV largest (algebraic) eigenvalues.
112
c 'SA' - compute the NEV smallest (algebraic) eigenvalues.
113
c 'LM' - compute the NEV largest (in magnitude) eigenvalues.
114
c 'SM' - compute the NEV smallest (in magnitude) eigenvalues.
115
c 'BE' - compute NEV eigenvalues, half from each end of the
116
c spectrum. When NEV is odd, compute one more from the
117
c high end than from the low end.
118
c (see remark 1 below)
120
c NEV Integer. (INPUT)
121
c Number of eigenvalues of OP to be computed. 0 < NEV < N.
123
c TOL Double precision scalar. (INPUT)
124
c Stopping criterion: the relative accuracy of the Ritz value
125
c is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I)).
126
c If TOL .LE. 0. is passed a default is set:
127
c DEFAULT = DLAMCH ('EPS') (machine precision as computed
128
c by the LAPACK auxiliary subroutine DLAMCH ).
130
c RESID Double precision array of length N. (INPUT/OUTPUT)
132
c If INFO .EQ. 0, a random initial residual vector is used.
133
c If INFO .NE. 0, RESID contains the initial residual vector,
134
c possibly from a previous run.
136
c RESID contains the final residual vector.
138
c NCV Integer. (INPUT)
139
c Number of columns of the matrix V (less than or equal to N).
140
c This will indicate how many Lanczos vectors are generated
141
c at each iteration. After the startup phase in which NEV
142
c Lanczos vectors are generated, the algorithm generates
143
c NCV-NEV Lanczos vectors at each subsequent update iteration.
144
c Most of the cost in generating each Lanczos vector is in the
145
c matrix-vector product OP*x. (See remark 4 below).
147
c V Double precision N by NCV array. (OUTPUT)
148
c The NCV columns of V contain the Lanczos basis vectors.
150
c LDV Integer. (INPUT)
151
c Leading dimension of V exactly as declared in the calling
154
c IPARAM Integer array of length 11. (INPUT/OUTPUT)
155
c IPARAM(1) = ISHIFT: method for selecting the implicit shifts.
156
c The shifts selected at each iteration are used to restart
157
c the Arnoldi iteration in an implicit fashion.
158
c -------------------------------------------------------------
159
c ISHIFT = 0: the shifts are provided by the user via
160
c reverse communication. The NCV eigenvalues of
161
c the current tridiagonal matrix T are returned in
162
c the part of WORKL array corresponding to RITZ.
163
c See remark 6 below.
164
c ISHIFT = 1: exact shifts with respect to the reduced
165
c tridiagonal matrix T. This is equivalent to
166
c restarting the iteration with a starting vector
167
c that is a linear combination of Ritz vectors
168
c associated with the "wanted" Ritz values.
169
c -------------------------------------------------------------
172
c No longer referenced. See remark 2 below.
175
c On INPUT: maximum number of Arnoldi update iterations allowed.
176
c On OUTPUT: actual number of Arnoldi update iterations taken.
178
c IPARAM(4) = NB: blocksize to be used in the recurrence.
179
c The code currently works only for NB = 1.
181
c IPARAM(5) = NCONV: number of "converged" Ritz values.
182
c This represents the number of Ritz values that satisfy
183
c the convergence criterion.
186
c No longer referenced. Implicit restarting is ALWAYS used.
189
c On INPUT determines what type of eigenproblem is being solved.
190
c Must be 1,2,3,4,5; See under \Description of dsaupd for the
191
c five modes available.
194
c When ido = 3 and the user provides shifts through reverse
195
c communication (IPARAM(1)=0), dsaupd returns NP, the number
196
c of shifts the user is to provide. 0 < NP <=NCV-NEV. See Remark
199
c IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,
200
c OUTPUT: NUMOP = total number of OP*x operations,
201
c NUMOPB = total number of B*x operations if BMAT='G',
202
c NUMREO = total number of steps of re-orthogonalization.
204
c IPNTR Integer array of length 11. (OUTPUT)
205
c Pointer to mark the starting locations in the WORKD and WORKL
206
c arrays for matrices/vectors used by the Lanczos iteration.
207
c -------------------------------------------------------------
208
c IPNTR(1): pointer to the current operand vector X in WORKD.
209
c IPNTR(2): pointer to the current result vector Y in WORKD.
210
c IPNTR(3): pointer to the vector B * X in WORKD when used in
211
c the shift-and-invert mode.
212
c IPNTR(4): pointer to the next available location in WORKL
213
c that is untouched by the program.
214
c IPNTR(5): pointer to the NCV by 2 tridiagonal matrix T in WORKL.
215
c IPNTR(6): pointer to the NCV RITZ values array in WORKL.
216
c IPNTR(7): pointer to the Ritz estimates in array WORKL associated
217
c with the Ritz values located in RITZ in WORKL.
218
c IPNTR(11): pointer to the NP shifts in WORKL. See Remark 6 below.
220
c Note: IPNTR(8:10) is only referenced by dseupd . See Remark 2.
221
c IPNTR(8): pointer to the NCV RITZ values of the original system.
222
c IPNTR(9): pointer to the NCV corresponding error bounds.
223
c IPNTR(10): pointer to the NCV by NCV matrix of eigenvectors
224
c of the tridiagonal matrix T. Only referenced by
225
c dseupd if RVEC = .TRUE. See Remarks.
226
c -------------------------------------------------------------
228
c WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION)
229
c Distributed array to be used in the basic Arnoldi iteration
230
c for reverse communication. The user should not use WORKD
231
c as temporary workspace during the iteration. Upon termination
232
c WORKD(1:N) contains B*RESID(1:N). If the Ritz vectors are desired
233
c subroutine dseupd uses this output.
234
c See Data Distribution Note below.
236
c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)
237
c Private (replicated) array on each PE or array allocated on
238
c the front end. See Data Distribution Note below.
240
c LWORKL Integer. (INPUT)
241
c LWORKL must be at least NCV**2 + 8*NCV .
243
c INFO Integer. (INPUT/OUTPUT)
244
c If INFO .EQ. 0, a randomly initial residual vector is used.
245
c If INFO .NE. 0, RESID contains the initial residual vector,
246
c possibly from a previous run.
247
c Error flag on output.
249
c = 1: Maximum number of iterations taken.
250
c All possible eigenvalues of OP has been found. IPARAM(5)
251
c returns the number of wanted converged Ritz values.
252
c = 2: No longer an informational error. Deprecated starting
253
c with release 2 of ARPACK.
254
c = 3: No shifts could be applied during a cycle of the
255
c Implicitly restarted Arnoldi iteration. One possibility
256
c is to increase the size of NCV relative to NEV.
257
c See remark 4 below.
258
c = -1: N must be positive.
259
c = -2: NEV must be positive.
260
c = -3: NCV must be greater than NEV and less than or equal to N.
261
c = -4: The maximum number of Arnoldi update iterations allowed
262
c must be greater than zero.
263
c = -5: WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.
264
c = -6: BMAT must be one of 'I' or 'G'.
265
c = -7: Length of private work array WORKL is not sufficient.
266
c = -8: Error return from trid. eigenvalue calculation;
267
c Informatinal error from LAPACK routine dsteqr .
268
c = -9: Starting vector is zero.
269
c = -10: IPARAM(7) must be 1,2,3,4,5.
270
c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
271
c = -12: IPARAM(1) must be equal to 0 or 1.
272
c = -13: NEV and WHICH = 'BE' are incompatable.
273
c = -9999: Could not build an Arnoldi factorization.
274
c IPARAM(5) returns the size of the current Arnoldi
275
c factorization. The user is advised to check that
276
c enough workspace and array storage has been allocated.
280
c 1. The converged Ritz values are always returned in ascending
281
c algebraic order. The computed Ritz values are approximate
282
c eigenvalues of OP. The selection of WHICH should be made
283
c with this in mind when Mode = 3,4,5. After convergence,
284
c approximate eigenvalues of the original problem may be obtained
285
c with the ARPACK subroutine dseupd .
287
c 2. If the Ritz vectors corresponding to the converged Ritz values
288
c are needed, the user must call dseupd immediately following completion
289
c of dsaupd . This is new starting with version 2.1 of ARPACK.
291
c 3. If M can be factored into a Cholesky factorization M = LL`
292
c then Mode = 2 should not be selected. Instead one should use
293
c Mode = 1 with OP = inv(L)*A*inv(L`). Appropriate triangular
294
c linear systems should be solved with L and L` rather
295
c than computing inverses. After convergence, an approximate
296
c eigenvector z of the original problem is recovered by solving
297
c L`z = x where x is a Ritz vector of OP.
299
c 4. At present there is no a-priori analysis to guide the selection
300
c of NCV relative to NEV. The only formal requrement is that NCV > NEV.
301
c However, it is recommended that NCV .ge. 2*NEV. If many problems of
302
c the same type are to be solved, one should experiment with increasing
303
c NCV while keeping NEV fixed for a given test problem. This will
304
c usually decrease the required number of OP*x operations but it
305
c also increases the work and storage required to maintain the orthogonal
306
c basis vectors. The optimal "cross-over" with respect to CPU time
307
c is problem dependent and must be determined empirically.
309
c 5. If IPARAM(7) = 2 then in the Reverse commuication interface the user
310
c must do the following. When IDO = 1, Y = OP * X is to be computed.
311
c When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user
312
c must overwrite X with A*X. Y is then the solution to the linear set
313
c of equations B*Y = A*X.
315
c 6. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the
316
c NP = IPARAM(8) shifts in locations:
318
c 2 WORKL(IPNTR(11)+1)
322
c NP WORKL(IPNTR(11)+NP-1).
324
c The eigenvalues of the current tridiagonal matrix are located in
325
c WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1). They are in the
326
c order defined by WHICH. The associated Ritz estimates are located in
327
c WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1).
329
c-----------------------------------------------------------------------
331
c\Data Distribution Note:
335
c REAL RESID(N), V(LDV,NCV), WORKD(3*N), WORKL(LWORKL)
336
c DECOMPOSE D1(N), D2(N,NCV)
337
c ALIGN RESID(I) with D1(I)
338
c ALIGN V(I,J) with D2(I,J)
339
c ALIGN WORKD(I) with D1(I) range (1:N)
340
c ALIGN WORKD(I) with D1(I-N) range (N+1:2*N)
341
c ALIGN WORKD(I) with D1(I-2*N) range (2*N+1:3*N)
342
c DISTRIBUTE D1(BLOCK), D2(BLOCK,:)
343
c REPLICATED WORKL(LWORKL)
347
c REAL RESID(N), V(LDV,NCV), WORKD(N,3), WORKL(LWORKL)
348
c SHARED RESID(BLOCK), V(BLOCK,:), WORKD(BLOCK,:)
349
c REPLICATED WORKL(LWORKL)
355
c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
356
c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
358
c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
359
c Restarted Arnoldi Iteration", Rice University Technical Report
360
c TR95-13, Department of Computational and Applied Mathematics.
361
c 3. B.N. Parlett, "The Symmetric Eigenvalue Problem". Prentice-Hall,
363
c 4. B.N. Parlett, B. Nour-Omid, "Towards a Black Box Lanczos Program",
364
c Computer Physics Communications, 53 (1989), pp 169-179.
365
c 5. B. Nour-Omid, B.N. Parlett, T. Ericson, P.S. Jensen, "How to
366
c Implement the Spectral Transformation", Math. Comp., 48 (1987),
368
c 6. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos
369
c Algorithm for Solving Sparse Symmetric Generalized Eigenproblems",
370
c SIAM J. Matr. Anal. Apps., January (1993).
371
c 7. L. Reichel, W.B. Gragg, "Algorithm 686: FORTRAN Subroutines
372
c for Updating the QR decomposition", ACM TOMS, December 1990,
373
c Volume 16 Number 4, pp 369-377.
374
c 8. R.B. Lehoucq, D.C. Sorensen, "Implementation of Some Spectral
375
c Transformations in a k-Step Arnoldi Method". In Preparation.
378
c dsaup2 ARPACK routine that implements the Implicitly Restarted
380
c dstats ARPACK routine that initialize timing and other statistics
382
c ivout ARPACK utility routine that prints integers.
383
c second ARPACK utility routine for timing.
384
c dvout ARPACK utility routine that prints vectors.
385
c dlamch LAPACK routine that determines machine constants.
388
c Danny Sorensen Phuong Vu
389
c Richard Lehoucq CRPC / Rice University
390
c Dept. of Computational & Houston, Texas
391
c Applied Mathematics
396
c 12/15/93: Version ' 2.4'
398
c\SCCS Information: @(#)
399
c FILE: saupd.F SID: 2.8 DATE OF SID: 04/10/01 RELEASE: 2
406
c-----------------------------------------------------------------------
409
& ( ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam,
410
& ipntr, workd, workl, lworkl, info )
412
c %----------------------------------------------------%
413
c | Include files for debugging and timing information |
414
c %----------------------------------------------------%
419
c %------------------%
420
c | Scalar Arguments |
421
c %------------------%
423
character bmat*1, which*2
424
integer ido, info, ldv, lworkl, n, ncv, nev
428
c %-----------------%
429
c | Array Arguments |
430
c %-----------------%
432
integer iparam(11), ipntr(11)
434
& resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
442
parameter (one = 1.0D+0 , zero = 0.0D+0 )
448
integer bounds, ierr, ih, iq, ishift, iupd, iw,
449
& ldh, ldq, msglvl, mxiter, mode, nb,
450
& nev0, next, np, ritz, j
451
save bounds, ierr, ih, iq, ishift, iupd, iw,
452
& ldh, ldq, msglvl, mxiter, mode, nb,
453
& nev0, next, np, ritz
455
c %----------------------%
456
c | External Subroutines |
457
c %----------------------%
459
external dsaup2 , dvout , ivout, second, dstats
461
c %--------------------%
462
c | External Functions |
463
c %--------------------%
469
c %-----------------------%
470
c | Executable Statements |
471
c %-----------------------%
475
c %-------------------------------%
476
c | Initialize timing statistics |
477
c | & message level for debugging |
478
c %-------------------------------%
490
c %--------------------------------------------%
491
c | Revision 2 performs only implicit restart. |
492
c %--------------------------------------------%
503
else if (nev .le. 0) then
505
else if (ncv .le. nev .or. ncv .gt. n) then
509
c %----------------------------------------------%
510
c | NP is the number of additional steps to |
511
c | extend the length NEV Lanczos factorization. |
512
c %----------------------------------------------%
516
if (mxiter .le. 0) ierr = -4
517
if (which .ne. 'LM' .and.
518
& which .ne. 'SM' .and.
519
& which .ne. 'LA' .and.
520
& which .ne. 'SA' .and.
521
& which .ne. 'BE') ierr = -5
522
if (bmat .ne. 'I' .and. bmat .ne. 'G') ierr = -6
524
if (lworkl .lt. ncv**2 + 8*ncv) ierr = -7
525
if (mode .lt. 1 .or. mode .gt. 5) then
527
else if (mode .eq. 1 .and. bmat .eq. 'G') then
529
else if (ishift .lt. 0 .or. ishift .gt. 1) then
531
else if (nev .eq. 1 .and. which .eq. 'BE') then
539
if (ierr .ne. 0) then
545
c %------------------------%
546
c | Set default parameters |
547
c %------------------------%
549
if (nb .le. 0) nb = 1
550
if (tol .le. zero) tol = dlamch ('EpsMach')
552
c %----------------------------------------------%
553
c | NP is the number of additional steps to |
554
c | extend the length NEV Lanczos factorization. |
555
c | NEV0 is the local variable designating the |
556
c | size of the invariant subspace desired. |
557
c %----------------------------------------------%
562
c %-----------------------------%
563
c | Zero out internal workspace |
564
c %-----------------------------%
566
do 10 j = 1, ncv**2 + 8*ncv
570
c %-------------------------------------------------------%
571
c | Pointer into WORKL for address of H, RITZ, BOUNDS, Q |
572
c | etc... and the remaining workspace. |
573
c | Also update pointer to be used on output. |
574
c | Memory is laid out as follows: |
575
c | workl(1:2*ncv) := generated tridiagonal matrix |
576
c | workl(2*ncv+1:2*ncv+ncv) := ritz values |
577
c | workl(3*ncv+1:3*ncv+ncv) := computed error bounds |
578
c | workl(4*ncv+1:4*ncv+ncv*ncv) := rotation matrix Q |
579
c | workl(4*ncv+ncv*ncv+1:7*ncv+ncv*ncv) := workspace |
580
c %-------------------------------------------------------%
598
c %-------------------------------------------------------%
599
c | Carry out the Implicitly restarted Lanczos Iteration. |
600
c %-------------------------------------------------------%
603
& ( ido, bmat, n, which, nev0, np, tol, resid, mode, iupd,
604
& ishift, mxiter, v, ldv, workl(ih), ldh, workl(ritz),
605
& workl(bounds), workl(iq), ldq, workl(iw), ipntr, workd,
608
c %--------------------------------------------------%
609
c | ido .ne. 99 implies use of reverse communication |
610
c | to compute operations involving OP or shifts. |
611
c %--------------------------------------------------%
613
if (ido .eq. 3) iparam(8) = np
614
if (ido .ne. 99) go to 9000
622
c %------------------------------------%
623
c | Exit if there was an informational |
624
c | error within dsaup2 . |
625
c %------------------------------------%
627
if (info .lt. 0) go to 9000
628
if (info .eq. 2) info = 3
630
if (msglvl .gt. 0) then
631
call ivout (logfil, 1, mxiter, ndigit,
632
& '_saupd: number of update iterations taken')
633
call ivout (logfil, 1, np, ndigit,
634
& '_saupd: number of "converged" Ritz values')
635
call dvout (logfil, np, workl(Ritz), ndigit,
636
& '_saupd: final Ritz values')
637
call dvout (logfil, np, workl(Bounds), ndigit,
638
& '_saupd: corresponding error bounds')
644
if (msglvl .gt. 0) then
646
c %--------------------------------------------------------%
647
c | Version Number & Version Date are defined in version.h |
648
c %--------------------------------------------------------%
651
write (6,1100) mxiter, nopx, nbx, nrorth, nitref, nrstrt,
652
& tmvopx, tmvbx, tsaupd, tsaup2, tsaitr, titref,
653
& tgetv0, tseigt, tsgets, tsapps, tsconv
655
& 5x, '==========================================',/
656
& 5x, '= Symmetric implicit Arnoldi update code =',/
657
& 5x, '= Version Number:', ' 2.4' , 19x, ' =',/
658
& 5x, '= Version Date: ', ' 07/31/96' , 14x, ' =',/
659
& 5x, '==========================================',/
660
& 5x, '= Summary of timing statistics =',/
661
& 5x, '==========================================',//)
663
& 5x, 'Total number update iterations = ', i5,/
664
& 5x, 'Total number of OP*x operations = ', i5,/
665
& 5x, 'Total number of B*x operations = ', i5,/
666
& 5x, 'Total number of reorthogonalization steps = ', i5,/
667
& 5x, 'Total number of iterative refinement steps = ', i5,/
668
& 5x, 'Total number of restart steps = ', i5,/
669
& 5x, 'Total time in user OP*x operation = ', f12.6,/
670
& 5x, 'Total time in user B*x operation = ', f12.6,/
671
& 5x, 'Total time in Arnoldi update routine = ', f12.6,/
672
& 5x, 'Total time in saup2 routine = ', f12.6,/
673
& 5x, 'Total time in basic Arnoldi iteration loop = ', f12.6,/
674
& 5x, 'Total time in reorthogonalization phase = ', f12.6,/
675
& 5x, 'Total time in (re)start vector generation = ', f12.6,/
676
& 5x, 'Total time in trid eigenvalue subproblem = ', f12.6,/
677
& 5x, 'Total time in getting the shifts = ', f12.6,/
678
& 5x, 'Total time in applying the shifts = ', f12.6,/
679
& 5x, 'Total time in convergence testing = ', f12.6)