Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
mydsaupd.f
1 c @HEADER
2 c *****************************************************************************
3 c Anasazi: Block Eigensolvers Package
4 c
5 c Copyright 2004 NTESS and the Anasazi contributors.
6 c SPDX-License-Identifier: BSD-3-Clause
7 c *****************************************************************************
8 c @HEADER
9 
10 c This software is a result of the research described in the report
11 c
12 c "A comparison of algorithms for modal analysis in the absence
13 c of a sparse direct method", P. Arbenz, R. Lehoucq, and U. Hetmaniuk,
14 c Sandia National Laboratories, Technical report SAND2003-1028J.
15 c
16 c It is based on the Epetra, AztecOO, and ML packages defined in the Trilinos
17 c framework ( http://trilinos.org/ ).
18 c-----------------------------------------------------------------------
19 c\BeginDoc
20 c
21 c\Name: mydsaupd
22 c
23 c\Description:
24 c
25 c Reverse communication interface for the Implicitly Restarted Arnoldi
26 c Iteration. For symmetric problems this reduces to a variant of the Lanczos
27 c method. This method has been designed to compute approximations to a
28 c few eigenpairs of a linear operator OP that is real and symmetric
29 c with respect to a real positive semi-definite symmetric matrix B,
30 c i.e.
31 c
32 c B*OP = (OP`)*B.
33 c
34 c Another way to express this condition is
35 c
36 c < x,OPy > = < OPx,y > where < z,w > = z`Bw .
37 c
38 c In the standard eigenproblem B is the identity matrix.
39 c ( A` denotes transpose of A)
40 c
41 c The computed approximate eigenvalues are called Ritz values and
42 c the corresponding approximate eigenvectors are called Ritz vectors.
43 c
44 c mydsaupd is usually called iteratively to solve one of the
45 c following problems:
46 c
47 c Mode 1: A*x = lambda*x, A symmetric
48 c ===> OP = A and B = I.
49 c
50 c Mode 2: A*x = lambda*M*x, A symmetric, M symmetric positive definite
51 c ===> OP = inv[M]*A and B = M.
52 c ===> (If M can be factored see remark 3 below)
53 c
54 c Mode 3: K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
55 c ===> OP = (inv[K - sigma*M])*M and B = M.
56 c ===> Shift-and-Invert mode
57 c
58 c Mode 4: K*x = lambda*KG*x, K symmetric positive semi-definite,
59 c KG symmetric indefinite
60 c ===> OP = (inv[K - sigma*KG])*K and B = K.
61 c ===> Buckling mode
62 c
63 c Mode 5: A*x = lambda*M*x, A symmetric, M symmetric positive semi-definite
64 c ===> OP = inv[A - sigma*M]*[A + sigma*M] and B = M.
65 c ===> Cayley transformed mode
66 c
67 c NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v
68 c should be accomplished either by a direct method
69 c using a sparse matrix factorization and solving
70 c
71 c [A - sigma*M]*w = v or M*w = v,
72 c
73 c or through an iterative method for solving these
74 c systems. If an iterative method is used, the
75 c convergence test must be more stringent than
76 c the accuracy requirements for the eigenvalue
77 c approximations.
78 c
79 c\Usage:
80 c call mydsaupd
81 c ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
82 c IPNTR, WORKD, WORKL, LWORKL, INFO, VERBOSE )
83 c
84 c\Arguments
85 c IDO Integer. (INPUT/OUTPUT)
86 c Reverse communication flag. IDO must be zero on the first
87 c call to dsaupd . IDO will be set internally to
88 c indicate the type of operation to be performed. Control is
89 c then given back to the calling routine which has the
90 c responsibility to carry out the requested operation and call
91 c dsaupd with the result. The operand is given in
92 c WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)).
93 c (If Mode = 2 see remark 5 below)
94 c -------------------------------------------------------------
95 c IDO = 0: first call to the reverse communication interface
96 c IDO = -1: compute Y = OP * X where
97 c IPNTR(1) is the pointer into WORKD for X,
98 c IPNTR(2) is the pointer into WORKD for Y.
99 c This is for the initialization phase to force the
100 c starting vector into the range of OP.
101 c IDO = 1: compute Y = OP * X where
102 c IPNTR(1) is the pointer into WORKD for X,
103 c IPNTR(2) is the pointer into WORKD for Y.
104 c In mode 3,4 and 5, the vector B * X is already
105 c available in WORKD(ipntr(3)). It does not
106 c need to be recomputed in forming OP * X.
107 c IDO = 2: compute Y = B * X where
108 c IPNTR(1) is the pointer into WORKD for X,
109 c IPNTR(2) is the pointer into WORKD for Y.
110 c IDO = 3: compute the IPARAM(8) shifts where
111 c IPNTR(11) is the pointer into WORKL for
112 c placing the shifts. See remark 6 below.
113 c IDO = 4: user checks convergence. See remark 7 below.
114 c IDO = 99: done
115 c -------------------------------------------------------------
116 c
117 c BMAT Character*1. (INPUT)
118 c BMAT specifies the type of the matrix B that defines the
119 c semi-inner product for the operator OP.
120 c B = 'I' -> standard eigenvalue problem A*x = lambda*x
121 c B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
122 c
123 c N Integer. (INPUT)
124 c Dimension of the eigenproblem.
125 c
126 c WHICH Character*2. (INPUT)
127 c Specify which of the Ritz values of OP to compute.
128 c
129 c 'LA' - compute the NEV largest (algebraic) eigenvalues.
130 c 'SA' - compute the NEV smallest (algebraic) eigenvalues.
131 c 'LM' - compute the NEV largest (in magnitude) eigenvalues.
132 c 'SM' - compute the NEV smallest (in magnitude) eigenvalues.
133 c 'BE' - compute NEV eigenvalues, half from each end of the
134 c spectrum. When NEV is odd, compute one more from the
135 c high end than from the low end.
136 c (see remark 1 below)
137 c
138 c NEV Integer. (INPUT)
139 c Number of eigenvalues of OP to be computed. 0 < NEV < N.
140 c
141 c TOL Double precision scalar. (INPUT)
142 c Stopping criterion: the relative accuracy of the Ritz value
143 c is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I)).
144 c If TOL .LE. 0. is passed a default is set:
145 c DEFAULT = DLAMCH ('EPS') (machine precision as computed
146 c by the LAPACK auxiliary subroutine DLAMCH ).
147 c Not referenced if IPARAM(4) is not equal to zero.
148 c See remark 7 below.
149 c
150 c RESID Double precision array of length N. (INPUT/OUTPUT)
151 c On INPUT:
152 c If INFO .EQ. 0, a random initial residual vector is used.
153 c If INFO .NE. 0, RESID contains the initial residual vector,
154 c possibly from a previous run.
155 c On OUTPUT:
156 c RESID contains the final residual vector.
157 c
158 c NCV Integer. (INPUT)
159 c Number of columns of the matrix V (less than or equal to N).
160 c This will indicate how many Lanczos vectors are generated
161 c at each iteration. After the startup phase in which NEV
162 c Lanczos vectors are generated, the algorithm generates
163 c NCV-NEV Lanczos vectors at each subsequent update iteration.
164 c Most of the cost in generating each Lanczos vector is in the
165 c matrix-vector product OP*x. (See remark 4 below).
166 c
167 c V Double precision N by NCV array. (OUTPUT)
168 c The NCV columns of V contain the Lanczos basis vectors.
169 c
170 c LDV Integer. (INPUT)
171 c Leading dimension of V exactly as declared in the calling
172 c program.
173 c
174 c IPARAM Integer array of length 11. (INPUT/OUTPUT)
175 c IPARAM(1) = ISHIFT: method for selecting the implicit shifts.
176 c The shifts selected at each iteration are used to restart
177 c the Arnoldi iteration in an implicit fashion.
178 c -------------------------------------------------------------
179 c ISHIFT = 0: the shifts are provided by the user via
180 c reverse communication. The NCV eigenvalues of
181 c the current tridiagonal matrix T are returned in
182 c the part of WORKL array corresponding to RITZ.
183 c See remark 6 below.
184 c ISHIFT = 1: exact shifts with respect to the reduced
185 c tridiagonal matrix T. This is equivalent to
186 c restarting the iteration with a starting vector
187 c that is a linear combination of Ritz vectors
188 c associated with the "wanted" Ritz values.
189 c -------------------------------------------------------------
190 c
191 c IPARAM(2) = LEVEC
192 c No longer referenced. See remark 2 below.
193 c
194 c IPARAM(3) = MXITER
195 c On INPUT: maximum number of Arnoldi update iterations allowed.
196 c On OUTPUT: actual number of Arnoldi update iterations taken.
197 c
198 c IPARAM(4) = USRCHK: a non-zero value denotes that the user will
199 c check convergence. See remark 7 below.
200 c
201 c IPARAM(5) = NCONV: number of "converged" Ritz values.
202 c If USRCHK is equal to zero, then, upon exit (IDO=99)
203 c NCONV is the number of Ritz values that satisfy
204 c the convergence criterion as determined by dsaupd.f
205 c If USRCHK is not equal to zero, then the user must set
206 c NCONV to be the number of Ritz values that satisfy the
207 c convergence criterion.
208 c
209 c IPARAM(6) = IUPD
210 c No longer referenced. Implicit restarting is ALWAYS used.
211 c
212 c IPARAM(7) = MODE
213 c On INPUT determines what type of eigenproblem is being solved.
214 c Must be 1,2,3,4,5; See under \Description of dsaupd for the
215 c five modes available.
216 c
217 c IPARAM(8) = NP
218 c When ido = 3 and the user provides shifts through reverse
219 c communication (IPARAM(1)=0), dsaupd returns NP, the number
220 c of shifts the user is to provide. 0 < NP <=NCV-NEV. See Remark
221 c 6 below.
222 c
223 c IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,
224 c OUTPUT: NUMOP = total number of OP*x operations,
225 c NUMOPB = total number of B*x operations if BMAT='G',
226 c NUMREO = total number of steps of re-orthogonalization.
227 c
228 c IPNTR Integer array of length 11. (OUTPUT)
229 c Pointer to mark the starting locations in the WORKD and WORKL
230 c arrays for matrices/vectors used by the Lanczos iteration.
231 c -------------------------------------------------------------
232 c IPNTR(1): pointer to the current operand vector X in WORKD.
233 c IPNTR(2): pointer to the current result vector Y in WORKD.
234 c IPNTR(3): pointer to the vector B * X in WORKD when used in
235 c the shift-and-invert mode.
236 c IPNTR(4): pointer to the next available location in WORKL
237 c that is untouched by the program.
238 c IPNTR(5): pointer to the NCV by 2 tridiagonal matrix T in WORKL.
239 c IPNTR(6): pointer to the NCV RITZ values array in WORKL.
240 c IPNTR(7): pointer to the Ritz estimates in array WORKL associated
241 c with the Ritz values located in RITZ in WORKL.
242 c IPNTR(11): pointer to the NP shifts in WORKL. See Remark 6 below.
243 c
244 c Note: IPNTR(8:10) is only referenced by dseupd . See Remark 2.
245 c IPNTR(8): pointer to the NCV RITZ values of the original system.
246 c IPNTR(9): pointer to the NCV corresponding error bounds.
247 c IPNTR(10): pointer to the NCV by NCV matrix of eigenvectors
248 c of the tridiagonal matrix T. Only referenced by
249 c dseupd if RVEC = .TRUE. See Remarks.
250 c -------------------------------------------------------------
251 c
252 c WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION)
253 c Distributed array to be used in the basic Arnoldi iteration
254 c for reverse communication. The user should not use WORKD
255 c as temporary workspace during the iteration. Upon termination
256 c WORKD(1:N) contains B*RESID(1:N). If the Ritz vectors are desired
257 c subroutine dseupd uses this output.
258 c See Data Distribution Note below.
259 c
260 c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)
261 c Private (replicated) array on each PE or array allocated on
262 c the front end. See Data Distribution Note below.
263 c
264 c LWORKL Integer. (INPUT)
265 c LWORKL must be at least NCV**2 + 8*NCV .
266 c
267 c INFO Integer. (INPUT/OUTPUT)
268 c If INFO .EQ. 0, a randomly initial residual vector is used.
269 c If INFO .NE. 0, RESID contains the initial residual vector,
270 c possibly from a previous run.
271 c Error flag on output.
272 c = 0: Normal exit.
273 c = 1: Maximum number of iterations taken.
274 c All possible eigenvalues of OP has been found. IPARAM(5)
275 c returns the number of wanted converged Ritz values.
276 c = 2: No longer an informational error. Deprecated starting
277 c with release 2 of ARPACK.
278 c = 3: No shifts could be applied during a cycle of the
279 c Implicitly restarted Arnoldi iteration. One possibility
280 c is to increase the size of NCV relative to NEV.
281 c See remark 4 below.
282 c = -1: N must be positive.
283 c = -2: NEV must be positive.
284 c = -3: NCV must be greater than NEV and less than or equal to N.
285 c = -4: The maximum number of Arnoldi update iterations allowed
286 c must be greater than zero.
287 c = -5: WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.
288 c = -6: BMAT must be one of 'I' or 'G'.
289 c = -7: Length of private work array WORKL is not sufficient.
290 c = -8: Error return from trid. eigenvalue calculation;
291 c Informatinal error from LAPACK routine dsteqr .
292 c = -9: Starting vector is zero.
293 c = -10: IPARAM(7) must be 1,2,3,4,5.
294 c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
295 c = -12: IPARAM(1) must be equal to 0 or 1.
296 c = -13: NEV and WHICH = 'BE' are incompatable.
297 c = -14: IPARAM(4) is nonzero and the user returned the number of
298 c converged Ritz values less than zero or greater and ncv.
299 c = -9999: Could not build an Arnoldi factorization.
300 c IPARAM(5) returns the size of the current Arnoldi
301 c factorization. The user is advised to check that
302 c enough workspace and array storage has been allocated.
303 c
304 c VERBOSE : Flag to print out the status of the computation per iteration
305 c
306 c\Remarks
307 c 1. The converged Ritz values are always returned in ascending
308 c algebraic order. The computed Ritz values are approximate
309 c eigenvalues of OP. The selection of WHICH should be made
310 c with this in mind when Mode = 3,4,5. After convergence,
311 c approximate eigenvalues of the original problem may be obtained
312 c with the ARPACK subroutine dseupd .
313 c
314 c 2. If the Ritz vectors corresponding to the converged Ritz values
315 c are needed, the user must call dseupd immediately following completion
316 c of dsaupd . This is new starting with version 2.1 of ARPACK.
317 c
318 c 3. If M can be factored into a Cholesky factorization M = LL`
319 c then Mode = 2 should not be selected. Instead one should use
320 c Mode = 1 with OP = inv(L)*A*inv(L`). Appropriate triangular
321 c linear systems should be solved with L and L` rather
322 c than computing inverses. After convergence, an approximate
323 c eigenvector z of the original problem is recovered by solving
324 c L`z = x where x is a Ritz vector of OP.
325 c
326 c 4. At present there is no a-priori analysis to guide the selection
327 c of NCV relative to NEV. The only formal requrement is that NCV > NEV.
328 c However, it is recommended that NCV .ge. 2*NEV. If many problems of
329 c the same type are to be solved, one should experiment with increasing
330 c NCV while keeping NEV fixed for a given test problem. This will
331 c usually decrease the required number of OP*x operations but it
332 c also increases the work and storage required to maintain the orthogonal
333 c basis vectors. The optimal "cross-over" with respect to CPU time
334 c is problem dependent and must be determined empirically.
335 c
336 c 5. If IPARAM(7) = 2 then in the Reverse commuication interface the user
337 c must do the following. When IDO = 1, Y = OP * X is to be computed.
338 c When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user
339 c must overwrite X with A*X. Y is then the solution to the linear set
340 c of equations B*Y = A*X.
341 c
342 c 6. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the
343 c NP = IPARAM(8) shifts in locations:
344 c 1 WORKL(IPNTR(11))
345 c 2 WORKL(IPNTR(11)+1)
346 c .
347 c .
348 c .
349 c NP WORKL(IPNTR(11)+NP-1).
350 c
351 c The eigenvalues of the current tridiagonal matrix are located in
352 c WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1). They are in the
353 c order defined by WHICH. The associated Ritz estimates are located in
354 c WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1).
355 c
356 c 7. When IPARAM(4)=USRCHK is not equal to zero and IDO=4, then the
357 c user will set IPARAM(5)=NCONV with the number of Ritz values that
358 c satisfy the users' convergence criterion before dsaupd.f is
359 c re-entered. Note that IPNTR(5) contains the pointer into WORKL that
360 c contains the NCV by 2 tridiagonal matrix T.
361 c
362 c-----------------------------------------------------------------------
363 c
364 c\Data Distribution Note:
365 c
366 c Fortran-D syntax:
367 c ================
368 c REAL RESID(N), V(LDV,NCV), WORKD(3*N), WORKL(LWORKL)
369 c DECOMPOSE D1(N), D2(N,NCV)
370 c ALIGN RESID(I) with D1(I)
371 c ALIGN V(I,J) with D2(I,J)
372 c ALIGN WORKD(I) with D1(I) range (1:N)
373 c ALIGN WORKD(I) with D1(I-N) range (N+1:2*N)
374 c ALIGN WORKD(I) with D1(I-2*N) range (2*N+1:3*N)
375 c DISTRIBUTE D1(BLOCK), D2(BLOCK,:)
376 c REPLICATED WORKL(LWORKL)
377 c
378 c Cray MPP syntax:
379 c ===============
380 c REAL RESID(N), V(LDV,NCV), WORKD(N,3), WORKL(LWORKL)
381 c SHARED RESID(BLOCK), V(BLOCK,:), WORKD(BLOCK,:)
382 c REPLICATED WORKL(LWORKL)
383 c
384 c
385 c\BeginLib
386 c
387 c\References:
388 c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
389 c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
390 c pp 357-385.
391 c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
392 c Restarted Arnoldi Iteration", Rice University Technical Report
393 c TR95-13, Department of Computational and Applied Mathematics.
394 c 3. B.N. Parlett, "The Symmetric Eigenvalue Problem". Prentice-Hall,
395 c 1980.
396 c 4. B.N. Parlett, B. Nour-Omid, "Towards a Black Box Lanczos Program",
397 c Computer Physics Communications, 53 (1989), pp 169-179.
398 c 5. B. Nour-Omid, B.N. Parlett, T. Ericson, P.S. Jensen, "How to
399 c Implement the Spectral Transformation", Math. Comp., 48 (1987),
400 c pp 663-673.
401 c 6. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos
402 c Algorithm for Solving Sparse Symmetric Generalized Eigenproblems",
403 c SIAM J. Matr. Anal. Apps., January (1993).
404 c 7. L. Reichel, W.B. Gragg, "Algorithm 686: FORTRAN Subroutines
405 c for Updating the QR decomposition", ACM TOMS, December 1990,
406 c Volume 16 Number 4, pp 369-377.
407 c 8. R.B. Lehoucq, D.C. Sorensen, "Implementation of Some Spectral
408 c Transformations in a k-Step Arnoldi Method". In Preparation.
409 c
410 c\Routines called:
411 c mydsaup2 ARPACK routine that implements the Implicitly Restarted
412 c Arnoldi Iteration.
413 c dstats ARPACK routine that initialize timing and other statistics
414 c variables.
415 c ivout ARPACK utility routine that prints integers.
416 c second ARPACK utility routine for timing.
417 c dvout ARPACK utility routine that prints vectors.
418 c dlamch LAPACK routine that determines machine constants.
419 c
420 c\Authors
421 c Danny Sorensen Phuong Vu
422 c Richard Lehoucq CRPC / Rice University
423 c Dept. of Computational & Houston, Texas
424 c Applied Mathematics
425 c Rice University
426 c Houston, Texas
427 c
428 c\Revision history:
429 c 12/15/93: Version ' 2.4'
430 c
431 c\SCCS Information: @(#)
432 c FILE: saupd.F SID: 2.8 DATE OF SID: 04/10/01 RELEASE: 2
433 c
434 c\Remarks
435 c 1. None
436 c
437 c\EndLib
438 c
439 c-----------------------------------------------------------------------
440 c
441  subroutine mydsaupd
442  & ( ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam,
443  & ipntr, workd, workl, lworkl, info, verbose )
444 c
445 c %----------------------------------------------------%
446 c | Include files for debugging and timing information |
447 c %----------------------------------------------------%
448 c
449  include 'debug.h'
450  include 'stat.h'
451 c
452 c %------------------%
453 c | Scalar Arguments |
454 c %------------------%
455 c
456  character bmat*1, which*2
457  integer ido, info, ldv, lworkl, n, ncv, nev
458  integer verbose
459  Double precision
460  & tol
461 c
462 c %-----------------%
463 c | Array Arguments |
464 c %-----------------%
465 c
466  integer iparam(11), ipntr(11)
467  Double precision
468  & resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
469 c
470 c %------------%
471 c | Parameters |
472 c %------------%
473 c
474  Double precision
475  & one, zero
476  parameter(one = 1.0d+0 , zero = 0.0d+0 )
477 c
478 c %---------------%
479 c | Local Scalars |
480 c %---------------%
481 c
482  integer bounds, ierr, ih, iq, ishift, iupd, iw,
483  & ldh, ldq, msglvl, mxiter, mode, nb,
484  & nev0, next, np, ritz, j, usrchk, unconv
485  save bounds, ierr, ih, iq, ishift, iupd, iw,
486  & ldh, ldq, msglvl, mxiter, mode, nb,
487  & nev0, next, np, ritz, usrchk, unconv
488 c
489 c %----------------------%
490 c | External Subroutines |
491 c %----------------------%
492 c
493  external mydsaup2 , dvout , ivout, second, dstats
494 c
495 c %--------------------%
496 c | External Functions |
497 c %--------------------%
498 c
499  Double precision
500  & dlamch
501  external dlamch
502 c
503 c %-----------------------%
504 c | Executable Statements |
505 c %-----------------------%
506 c
507  if (ido .eq. 0) then
508 c
509 c %-------------------------------%
510 c | Initialize timing statistics |
511 c | & message level for debugging |
512 c %-------------------------------%
513 c
514  call dstats
515  call second(t0)
516  msglvl = msaupd
517 c
518  ierr = 0
519  ishift = iparam(1)
520  mxiter = iparam(3)
521 c nb = iparam(4)
522  nb = 1
523  usrchk = iparam(4)
524  unconv = 0
525 c
526 c %--------------------------------------------%
527 c | Revision 2 performs only implicit restart. |
528 c %--------------------------------------------%
529 c
530  iupd = 1
531  mode = iparam(7)
532 c
533 c %----------------%
534 c | Error checking |
535 c %----------------%
536 c
537  if (n .le. 0) then
538  ierr = -1
539  else if (nev .le. 0) then
540  ierr = -2
541  else if (ncv .le. nev .or. ncv .gt. n) then
542  ierr = -3
543  end if
544 c
545 c %----------------------------------------------%
546 c | NP is the number of additional steps to |
547 c | extend the length NEV Lanczos factorization. |
548 c %----------------------------------------------%
549 c
550  np = ncv - nev
551 c
552  if (mxiter .le. 0) ierr = -4
553  if (which .ne. 'LM' .and.
554  & which .ne. 'SM' .and.
555  & which .ne. 'LA' .and.
556  & which .ne. 'SA' .and.
557  & which .ne. 'BE') ierr = -5
558  if (bmat .ne. 'I' .and. bmat .ne. 'G') ierr = -6
559 c
560  if (lworkl .lt. ncv**2 + 8*ncv) ierr = -7
561  if (mode .lt. 1 .or. mode .gt. 5) then
562  ierr = -10
563  else if (mode .eq. 1 .and. bmat .eq. 'G') then
564  ierr = -11
565  else if (ishift .lt. 0 .or. ishift .gt. 1) then
566  ierr = -12
567  else if (nev .eq. 1 .and. which .eq. 'BE') then
568  ierr = -13
569  end if
570 c
571 c %------------%
572 c | Error Exit |
573 c %------------%
574 c
575  if (ierr .ne. 0) then
576  info = ierr
577  ido = 99
578  go to 9000
579  end if
580 c
581 c %------------------------%
582 c | Set default parameters |
583 c %------------------------%
584 c
585  if (nb .le. 0) nb = 1
586  if (tol .le. zero .and. usrchk .eq. 0)
587  & tol = dlamch('EpsMach')
588 c
589 c %----------------------------------------------%
590 c | NP is the number of additional steps to |
591 c | extend the length NEV Lanczos factorization. |
592 c | NEV0 is the local variable designating the |
593 c | size of the invariant subspace desired. |
594 c %----------------------------------------------%
595 c
596  np = ncv - nev
597  nev0 = nev
598 c
599 c %-----------------------------%
600 c | Zero out internal workspace |
601 c %-----------------------------%
602 c
603  do 10 j = 1, ncv**2 + 8*ncv
604  workl(j) = zero
605  10 continue
606 c
607 c %-------------------------------------------------------%
608 c | Pointer into WORKL for address of H, RITZ, BOUNDS, Q |
609 c | etc... and the remaining workspace. |
610 c | Also update pointer to be used on output. |
611 c | Memory is laid out as follows: |
612 c | workl(1:2*ncv) := generated tridiagonal matrix |
613 c | workl(2*ncv+1:2*ncv+ncv) := ritz values |
614 c | workl(3*ncv+1:3*ncv+ncv) := computed error bounds |
615 c | workl(4*ncv+1:4*ncv+ncv*ncv) := rotation matrix Q |
616 c | workl(4*ncv+ncv*ncv+1:7*ncv+ncv*ncv) := workspace |
617 c %-------------------------------------------------------%
618 c
619  ldh = ncv
620  ldq = ncv
621  ih = 1
622  ritz = ih + 2*ldh
623  bounds = ritz + ncv
624  iq = bounds + ncv
625  iw = iq + ncv**2
626  next = iw + 3*ncv
627 c
628  ipntr(4) = next
629  ipntr(5) = ih
630  ipntr(6) = ritz
631  ipntr(7) = bounds
632  ipntr(11) = iw
633  end if
634 c
635  if (ido .eq. 4) then
636 c
637 c %------------------------------------------------------------%
638 c | Returning from reverse communication; iparam(5) contains |
639 c | the number of Ritz values that satisfy the users criterion |
640 c %------------------------------------------------------------%
641 c
642  unconv = iparam(5)
643  if (unconv .lt. 0 .or. unconv .gt. ncv ) then
644  ierr = -14
645  info = ierr
646  ido = 99
647  go to 9000
648  end if
649  end if
650 c
651 c %-------------------------------------------------------%
652 c | Carry out the Implicitly restarted Lanczos Iteration. |
653 c %-------------------------------------------------------%
654 c
655  call mydsaup2
656  & ( ido, bmat, n, which, nev0, np, tol, resid, mode, iupd,
657  & ishift, mxiter, v, ldv, workl(ih), ldh, workl(ritz),
658  & workl(bounds), workl(iq), ldq, workl(iw), ipntr, workd,
659  & info, verbose, usrchk, unconv )
660 c
661 c %--------------------------------------------------%
662 c | ido .ne. 99 implies use of reverse communication |
663 c | to compute operations involving OP or shifts. |
664 c %--------------------------------------------------%
665 c
666  if (ido .eq. 3) iparam(8) = np
667  if (ido .ne. 99) go to 9000
668 c
669  iparam(3) = mxiter
670  iparam(5) = np
671  iparam(9) = nopx
672  iparam(10) = nbx
673  iparam(11) = nrorth
674 c
675 c %------------------------------------%
676 c | Exit if there was an informational |
677 c | error within mydsaup2 . |
678 c %------------------------------------%
679 c
680  if (info .lt. 0) go to 9000
681  if (info .eq. 2) info = 3
682 c
683  if (msglvl .gt. 0) then
684  call ivout(logfil, 1, mxiter, ndigit,
685  & '_saupd: number of update iterations taken')
686  call ivout(logfil, 1, np, ndigit,
687  & '_saupd: number of "converged" Ritz values')
688  call dvout(logfil, np, workl(ritz), ndigit,
689  & '_saupd: final Ritz values')
690  call dvout(logfil, np, workl(bounds), ndigit,
691  & '_saupd: corresponding error bounds')
692  end if
693 c
694  call second(t1)
695  tsaupd = t1 - t0
696 c
697  if (msglvl .gt. 0) then
698 c
699 c %--------------------------------------------------------%
700 c | Version Number & Version Date are defined in version.h |
701 c %--------------------------------------------------------%
702 c
703  write (6,1000)
704  write (6,1100) mxiter, nopx, nbx, nrorth, nitref, nrstrt,
705  & tmvopx, tmvbx, tsaupd, tsaup2, tsaitr, titref,
706  & tgetv0, tseigt, tsgets, tsapps, tsconv
707  1000 format (//,
708  & 5x, '==========================================',/
709  & 5x, '= Symmetric implicit Arnoldi update code =',/
710  & 5x, '= Version Number:', ' 2.4' , 19x, ' =',/
711  & 5x, '= Version Date: ', ' 07/31/96' , 14x, ' =',/
712  & 5x, '==========================================',/
713  & 5x, '= Summary of timing statistics =',/
714  & 5x, '==========================================',//)
715  1100 format (
716  & 5x, 'Total number update iterations = ', i5,/
717  & 5x, 'Total number of OP*x operations = ', i5,/
718  & 5x, 'Total number of B*x operations = ', i5,/
719  & 5x, 'Total number of reorthogonalization steps = ', i5,/
720  & 5x, 'Total number of iterative refinement steps = ', i5,/
721  & 5x, 'Total number of restart steps = ', i5,/
722  & 5x, 'Total time in user OP*x operation = ', f12.6,/
723  & 5x, 'Total time in user B*x operation = ', f12.6,/
724  & 5x, 'Total time in Arnoldi update routine = ', f12.6,/
725  & 5x, 'Total time in saup2 routine = ', f12.6,/
726  & 5x, 'Total time in basic Arnoldi iteration loop = ', f12.6,/
727  & 5x, 'Total time in reorthogonalization phase = ', f12.6,/
728  & 5x, 'Total time in (re)start vector generation = ', f12.6,/
729  & 5x, 'Total time in trid eigenvalue subproblem = ', f12.6,/
730  & 5x, 'Total time in getting the shifts = ', f12.6,/
731  & 5x, 'Total time in applying the shifts = ', f12.6,/
732  & 5x, 'Total time in convergence testing = ', f12.6)
733  end if
734 c
735  9000 continue
736 c
737  return
738 c
739 c %---------------%
740 c | End of dsaupd |
741 c %---------------%
742 c
743  end
Generated on Fri Jul 25 2025 09:18:08 for Anasazi by   doxygen 1.8.5