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