Mercurial > octave-nkf
comparison libcruft/arpack/src/znaupd.f @ 12274:9f5d2ef078e8 release-3-4-x
import ARPACK sources to libcruft from Debian package libarpack2 2.1+parpack96.dfsg-3+b1
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 28 Jan 2011 14:04:33 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
12273:83133b5bf392 | 12274:9f5d2ef078e8 |
---|---|
1 c\BeginDoc | |
2 c | |
3 c\Name: znaupd | |
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 dsaupd or dnaupd 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 znaupd 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 hermitian 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 hermitian 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 znaupd | |
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 znaupd . 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 znaupd 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 Double precision 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 = dlamch ('EPS') (machine precision as computed | |
108 c by the LAPACK auxiliary subroutine dlamch ). | |
109 c | |
110 c RESID Complex*16 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*16 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; See under \Description of znaupd 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 zneupd . 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 zneupd if RVEC = .TRUE. See Remark 2 below. | |
211 c | |
212 c ------------------------------------------------------------- | |
213 c | |
214 c WORKD Complex*16 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*16 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 Double precision 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 incompatible. | |
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 zneupd . | |
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 zneupd immediately following | |
276 c completion of znaupd . 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 requirement 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*16 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*16 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*16 | |
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 znaup2 ARPACK routine that implements the Implicitly Restarted | |
355 c Arnoldi Iteration. | |
356 c zstatn ARPACK routine that initializes the timing variables. | |
357 c ivout ARPACK utility routine that prints integers. | |
358 c zvout ARPACK utility routine that prints vectors. | |
359 c arscnd ARPACK utility routine for timing. | |
360 c dlamch 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.8 DATE OF SID: 04/10/01 RELEASE: 2 | |
372 c | |
373 c\Remarks | |
374 c | |
375 c\EndLib | |
376 c | |
377 c----------------------------------------------------------------------- | |
378 c | |
379 subroutine znaupd | |
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 Double precision | |
397 & tol | |
398 c | |
399 c %-----------------% | |
400 c | Array Arguments | | |
401 c %-----------------% | |
402 c | |
403 integer iparam(11), ipntr(14) | |
404 Complex*16 | |
405 & resid(n), v(ldv,ncv), workd(3*n), workl(lworkl) | |
406 Double precision | |
407 & rwork(ncv) | |
408 c | |
409 c %------------% | |
410 c | Parameters | | |
411 c %------------% | |
412 c | |
413 Complex*16 | |
414 & one, zero | |
415 parameter (one = (1.0D+0, 0.0D+0) , zero = (0.0D+0, 0.0D+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 znaup2 , zvout , ivout, arscnd, zstatn | |
433 c | |
434 c %--------------------% | |
435 c | External Functions | | |
436 c %--------------------% | |
437 c | |
438 Double precision | |
439 & dlamch | |
440 external dlamch | |
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 zstatn | |
454 call arscnd (t0) | |
455 msglvl = mcaupd | |
456 c | |
457 c %----------------% | |
458 c | Error checking | | |
459 c %----------------% | |
460 c | |
461 ierr = 0 | |
462 ishift = iparam(1) | |
463 c levec = iparam(2) | |
464 mxiter = iparam(3) | |
465 c nb = iparam(4) | |
466 nb = 1 | |
467 c | |
468 c %--------------------------------------------% | |
469 c | Revision 2 performs only implicit restart. | | |
470 c %--------------------------------------------% | |
471 c | |
472 iupd = 1 | |
473 mode = iparam(7) | |
474 c | |
475 if (n .le. 0) then | |
476 ierr = -1 | |
477 else if (nev .le. 0) then | |
478 ierr = -2 | |
479 else if (ncv .le. nev .or. ncv .gt. n) then | |
480 ierr = -3 | |
481 else if (mxiter .le. 0) then | |
482 ierr = -4 | |
483 else if (which .ne. 'LM' .and. | |
484 & which .ne. 'SM' .and. | |
485 & which .ne. 'LR' .and. | |
486 & which .ne. 'SR' .and. | |
487 & which .ne. 'LI' .and. | |
488 & which .ne. 'SI') then | |
489 ierr = -5 | |
490 else if (bmat .ne. 'I' .and. bmat .ne. 'G') then | |
491 ierr = -6 | |
492 else if (lworkl .lt. 3*ncv**2 + 5*ncv) then | |
493 ierr = -7 | |
494 else if (mode .lt. 1 .or. mode .gt. 3) then | |
495 ierr = -10 | |
496 else if (mode .eq. 1 .and. bmat .eq. 'G') then | |
497 ierr = -11 | |
498 end if | |
499 c | |
500 c %------------% | |
501 c | Error Exit | | |
502 c %------------% | |
503 c | |
504 if (ierr .ne. 0) then | |
505 info = ierr | |
506 ido = 99 | |
507 go to 9000 | |
508 end if | |
509 c | |
510 c %------------------------% | |
511 c | Set default parameters | | |
512 c %------------------------% | |
513 c | |
514 if (nb .le. 0) nb = 1 | |
515 if (tol .le. 0.0D+0 ) tol = dlamch ('EpsMach') | |
516 if (ishift .ne. 0 .and. | |
517 & ishift .ne. 1 .and. | |
518 & ishift .ne. 2) ishift = 1 | |
519 c | |
520 c %----------------------------------------------% | |
521 c | NP is the number of additional steps to | | |
522 c | extend the length NEV Lanczos factorization. | | |
523 c | NEV0 is the local variable designating the | | |
524 c | size of the invariant subspace desired. | | |
525 c %----------------------------------------------% | |
526 c | |
527 np = ncv - nev | |
528 nev0 = nev | |
529 c | |
530 c %-----------------------------% | |
531 c | Zero out internal workspace | | |
532 c %-----------------------------% | |
533 c | |
534 do 10 j = 1, 3*ncv**2 + 5*ncv | |
535 workl(j) = zero | |
536 10 continue | |
537 c | |
538 c %-------------------------------------------------------------% | |
539 c | Pointer into WORKL for address of H, RITZ, BOUNDS, Q | | |
540 c | etc... and the remaining workspace. | | |
541 c | Also update pointer to be used on output. | | |
542 c | Memory is laid out as follows: | | |
543 c | workl(1:ncv*ncv) := generated Hessenberg matrix | | |
544 c | workl(ncv*ncv+1:ncv*ncv+ncv) := the ritz values | | |
545 c | workl(ncv*ncv+ncv+1:ncv*ncv+2*ncv) := error bounds | | |
546 c | workl(ncv*ncv+2*ncv+1:2*ncv*ncv+2*ncv) := rotation matrix Q | | |
547 c | workl(2*ncv*ncv+2*ncv+1:3*ncv*ncv+5*ncv) := workspace | | |
548 c | The final workspace is needed by subroutine zneigh called | | |
549 c | by znaup2 . Subroutine zneigh calls LAPACK routines for | | |
550 c | calculating eigenvalues and the last row of the eigenvector | | |
551 c | matrix. | | |
552 c %-------------------------------------------------------------% | |
553 c | |
554 ldh = ncv | |
555 ldq = ncv | |
556 ih = 1 | |
557 ritz = ih + ldh*ncv | |
558 bounds = ritz + ncv | |
559 iq = bounds + ncv | |
560 iw = iq + ldq*ncv | |
561 next = iw + ncv**2 + 3*ncv | |
562 c | |
563 ipntr(4) = next | |
564 ipntr(5) = ih | |
565 ipntr(6) = ritz | |
566 ipntr(7) = iq | |
567 ipntr(8) = bounds | |
568 ipntr(14) = iw | |
569 end if | |
570 c | |
571 c %-------------------------------------------------------% | |
572 c | Carry out the Implicitly restarted Arnoldi Iteration. | | |
573 c %-------------------------------------------------------% | |
574 c | |
575 call znaup2 | |
576 & ( ido, bmat, n, which, nev0, np, tol, resid, mode, iupd, | |
577 & ishift, mxiter, v, ldv, workl(ih), ldh, workl(ritz), | |
578 & workl(bounds), workl(iq), ldq, workl(iw), | |
579 & ipntr, workd, rwork, info ) | |
580 c | |
581 c %--------------------------------------------------% | |
582 c | ido .ne. 99 implies use of reverse communication | | |
583 c | to compute operations involving OP. | | |
584 c %--------------------------------------------------% | |
585 c | |
586 if (ido .eq. 3) iparam(8) = np | |
587 if (ido .ne. 99) go to 9000 | |
588 c | |
589 iparam(3) = mxiter | |
590 iparam(5) = np | |
591 iparam(9) = nopx | |
592 iparam(10) = nbx | |
593 iparam(11) = nrorth | |
594 c | |
595 c %------------------------------------% | |
596 c | Exit if there was an informational | | |
597 c | error within znaup2 . | | |
598 c %------------------------------------% | |
599 c | |
600 if (info .lt. 0) go to 9000 | |
601 if (info .eq. 2) info = 3 | |
602 c | |
603 if (msglvl .gt. 0) then | |
604 call ivout (logfil, 1, mxiter, ndigit, | |
605 & '_naupd: Number of update iterations taken') | |
606 call ivout (logfil, 1, np, ndigit, | |
607 & '_naupd: Number of wanted "converged" Ritz values') | |
608 call zvout (logfil, np, workl(ritz), ndigit, | |
609 & '_naupd: The final Ritz values') | |
610 call zvout (logfil, np, workl(bounds), ndigit, | |
611 & '_naupd: Associated Ritz estimates') | |
612 end if | |
613 c | |
614 call arscnd (t1) | |
615 tcaupd = t1 - t0 | |
616 c | |
617 if (msglvl .gt. 0) then | |
618 c | |
619 c %--------------------------------------------------------% | |
620 c | Version Number & Version Date are defined in version.h | | |
621 c %--------------------------------------------------------% | |
622 c | |
623 write (6,1000) | |
624 write (6,1100) mxiter, nopx, nbx, nrorth, nitref, nrstrt, | |
625 & tmvopx, tmvbx, tcaupd, tcaup2, tcaitr, titref, | |
626 & tgetv0, tceigh, tcgets, tcapps, tcconv, trvec | |
627 1000 format (//, | |
628 & 5x, '=============================================',/ | |
629 & 5x, '= Complex implicit Arnoldi update code =',/ | |
630 & 5x, '= Version Number: ', ' 2.3' , 21x, ' =',/ | |
631 & 5x, '= Version Date: ', ' 07/31/96' , 16x, ' =',/ | |
632 & 5x, '=============================================',/ | |
633 & 5x, '= Summary of timing statistics =',/ | |
634 & 5x, '=============================================',//) | |
635 1100 format ( | |
636 & 5x, 'Total number update iterations = ', i5,/ | |
637 & 5x, 'Total number of OP*x operations = ', i5,/ | |
638 & 5x, 'Total number of B*x operations = ', i5,/ | |
639 & 5x, 'Total number of reorthogonalization steps = ', i5,/ | |
640 & 5x, 'Total number of iterative refinement steps = ', i5,/ | |
641 & 5x, 'Total number of restart steps = ', i5,/ | |
642 & 5x, 'Total time in user OP*x operation = ', f12.6,/ | |
643 & 5x, 'Total time in user B*x operation = ', f12.6,/ | |
644 & 5x, 'Total time in Arnoldi update routine = ', f12.6,/ | |
645 & 5x, 'Total time in naup2 routine = ', f12.6,/ | |
646 & 5x, 'Total time in basic Arnoldi iteration loop = ', f12.6,/ | |
647 & 5x, 'Total time in reorthogonalization phase = ', f12.6,/ | |
648 & 5x, 'Total time in (re)start vector generation = ', f12.6,/ | |
649 & 5x, 'Total time in Hessenberg eig. subproblem = ', f12.6,/ | |
650 & 5x, 'Total time in getting the shifts = ', f12.6,/ | |
651 & 5x, 'Total time in applying the shifts = ', f12.6,/ | |
652 & 5x, 'Total time in convergence testing = ', f12.6,/ | |
653 & 5x, 'Total time in computing final Ritz vectors = ', f12.6/) | |
654 end if | |
655 c | |
656 9000 continue | |
657 c | |
658 return | |
659 c | |
660 c %---------------% | |
661 c | End of znaupd | | |
662 c %---------------% | |
663 c | |
664 end |