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