1c\BeginDoc 2c 3c\Name: snaup2 4c 5c\Description: 6c Intermediate level interface called by snaupd. 7c 8c\Usage: 9c call snaup2 10c ( IDO, BMAT, N, WHICH, NEV, NP, TOL, RESID, MODE, IUPD, 11c ISHIFT, MXITER, V, LDV, H, LDH, RITZR, RITZI, BOUNDS, 12c Q, LDQ, WORKL, IPNTR, WORKD, INFO ) 13c 14c\Arguments 15c 16c IDO, BMAT, N, WHICH, NEV, TOL, RESID: same as defined in snaupd. 17c MODE, ISHIFT, MXITER: see the definition of IPARAM in snaupd. 18c 19c NP Integer. (INPUT/OUTPUT) 20c Contains the number of implicit shifts to apply during 21c each Arnoldi iteration. 22c If ISHIFT=1, NP is adjusted dynamically at each iteration 23c to accelerate convergence and prevent stagnation. 24c This is also roughly equal to the number of matrix-vector 25c products (involving the operator OP) per Arnoldi iteration. 26c The logic for adjusting is contained within the current 27c subroutine. 28c If ISHIFT=0, NP is the number of shifts the user needs 29c to provide via reverse comunication. 0 < NP < NCV-NEV. 30c NP may be less than NCV-NEV for two reasons. The first, is 31c to keep complex conjugate pairs of "wanted" Ritz values 32c together. The second, is that a leading block of the current 33c upper Hessenberg matrix has split off and contains "unwanted" 34c Ritz values. 35c Upon termination of the IRA iteration, NP contains the number 36c of "converged" wanted Ritz values. 37c 38c IUPD Integer. (INPUT) 39c IUPD .EQ. 0: use explicit restart instead implicit update. 40c IUPD .NE. 0: use implicit update. 41c 42c V Real N by (NEV+NP) array. (INPUT/OUTPUT) 43c The Arnoldi basis vectors are returned in the first NEV 44c columns of V. 45c 46c LDV Integer. (INPUT) 47c Leading dimension of V exactly as declared in the calling 48c program. 49c 50c H Real (NEV+NP) by (NEV+NP) array. (OUTPUT) 51c H is used to store the generated upper Hessenberg matrix 52c 53c LDH Integer. (INPUT) 54c Leading dimension of H exactly as declared in the calling 55c program. 56c 57c RITZR, Real arrays of length NEV+NP. (OUTPUT) 58c RITZI RITZR(1:NEV) (resp. RITZI(1:NEV)) contains the real (resp. 59c imaginary) part of the computed Ritz values of OP. 60c 61c BOUNDS Real array of length NEV+NP. (OUTPUT) 62c BOUNDS(1:NEV) contain the error bounds corresponding to 63c the computed Ritz values. 64c 65c Q Real (NEV+NP) by (NEV+NP) array. (WORKSPACE) 66c Private (replicated) work array used to accumulate the 67c rotation in the shift application step. 68c 69c LDQ Integer. (INPUT) 70c Leading dimension of Q exactly as declared in the calling 71c program. 72c 73c WORKL Real work array of length at least 74c (NEV+NP)**2 + 3*(NEV+NP). (INPUT/WORKSPACE) 75c Private (replicated) array on each PE or array allocated on 76c the front end. It is used in shifts calculation, shifts 77c application and convergence checking. 78c 79c On exit, the last 3*(NEV+NP) locations of WORKL contain 80c the Ritz values (real,imaginary) and associated Ritz 81c estimates of the current Hessenberg matrix. They are 82c listed in the same order as returned from sneigh. 83c 84c If ISHIFT .EQ. O and IDO .EQ. 3, the first 2*NP locations 85c of WORKL are used in reverse communication to hold the user 86c supplied shifts. 87c 88c IPNTR Integer array of length 3. (OUTPUT) 89c Pointer to mark the starting locations in the WORKD for 90c vectors used by the Arnoldi iteration. 91c ------------------------------------------------------------- 92c IPNTR(1): pointer to the current operand vector X. 93c IPNTR(2): pointer to the current result vector Y. 94c IPNTR(3): pointer to the vector B * X when used in the 95c shift-and-invert mode. X is the current operand. 96c ------------------------------------------------------------- 97c 98c WORKD Real work array of length 3*N. (WORKSPACE) 99c Distributed array to be used in the basic Arnoldi iteration 100c for reverse communication. The user should not use WORKD 101c as temporary workspace during the iteration !!!!!!!!!! 102c See Data Distribution Note in DNAUPD. 103c 104c INFO Integer. (INPUT/OUTPUT) 105c If INFO .EQ. 0, a randomly initial residual vector is used. 106c If INFO .NE. 0, RESID contains the initial residual vector, 107c possibly from a previous run. 108c Error flag on output. 109c = 0: Normal return. 110c = 1: Maximum number of iterations taken. 111c All possible eigenvalues of OP has been found. 112c NP returns the number of converged Ritz values. 113c = 2: No shifts could be applied. 114c = -8: Error return from LAPACK eigenvalue calculation; 115c This should never happen. 116c = -9: Starting vector is zero. 117c = -9999: Could not build an Arnoldi factorization. 118c Size that was built in returned in NP. 119c 120c\EndDoc 121c 122c----------------------------------------------------------------------- 123c 124c\BeginLib 125c 126c\Local variables: 127c xxxxxx real 128c 129c\References: 130c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in 131c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), 132c pp 357-385. 133c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly 134c Restarted Arnoldi Iteration", Rice University Technical Report 135c TR95-13, Department of Computational and Applied Mathematics. 136c 137c\Routines called: 138c sgetv0 ARPACK initial vector generation routine. 139c snaitr ARPACK Arnoldi factorization routine. 140c snapps ARPACK application of implicit shifts routine. 141c snconv ARPACK convergence of Ritz values routine. 142c sneigh ARPACK compute Ritz values and error bounds routine. 143c sngets ARPACK reorder Ritz values and error bounds routine. 144c ssortc ARPACK sorting routine. 145c ivout ARPACK utility routine that prints integers. 146c second ARPACK utility routine for timing. 147c smout ARPACK utility routine that prints matrices 148c svout ARPACK utility routine that prints vectors. 149c slamch LAPACK routine that determines machine constants. 150c slapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully. 151c scopy Level 1 BLAS that copies one vector to another . 152c sdot Level 1 BLAS that computes the scalar product of two vectors. 153c snrm2 Level 1 BLAS that computes the norm of a vector. 154c sswap Level 1 BLAS that swaps two vectors. 155c 156c\Author 157c Danny Sorensen Phuong Vu 158c Richard Lehoucq CRPC / Rice University 159c Dept. of Computational & Houston, Texas 160c Applied Mathematics 161c Rice University 162c Houston, Texas 163c 164c\SCCS Information: @(#) 165c FILE: naup2.F SID: 2.4 DATE OF SID: 7/30/96 RELEASE: 2 166c 167c\Remarks 168c 1. None 169c 170c\EndLib 171c 172c----------------------------------------------------------------------- 173c 174 subroutine snaup2 175 & ( ido, bmat, n, which, nev, np, tol, resid, mode, iupd, 176 & ishift, mxiter, v, ldv, h, ldh, ritzr, ritzi, bounds, 177 & q, ldq, workl, ipntr, workd, info ) 178c 179c %----------------------------------------------------% 180c | Include files for debugging and timing information | 181c %----------------------------------------------------% 182c 183 184c 185c\SCCS Information: @(#) 186c FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 187c 188c %---------------------------------% 189c | See debug.doc for documentation | 190c %---------------------------------% 191 integer logfil, ndigit, mgetv0, 192 & msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd, 193 & mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd, 194 & mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd 195 common /debug/ 196 & logfil, ndigit, mgetv0, 197 & msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd, 198 & mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd, 199 & mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd 200 201c %--------------------------------% 202c | See stat.doc for documentation | 203c %--------------------------------% 204c 205c\SCCS Information: @(#) 206c FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 207c 208 real t0, t1, t2, t3, t4, t5 209c save t0, t1, t2, t3, t4, t5 210c 211 integer nopx, nbx, nrorth, nitref, nrstrt 212 real tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, 213 & tnaupd, tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, 214 & tcaupd, tcaup2, tcaitr, tceigh, tcgets, tcapps, tcconv, 215 & tmvopx, tmvbx, tgetv0, titref, trvec 216 common /timing/ 217 & nopx, nbx, nrorth, nitref, nrstrt, 218 & tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, 219 & tnaupd, tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, 220 & tcaupd, tcaup2, tcaitr, tceigh, tcgets, tcapps, tcconv, 221 & tmvopx, tmvbx, tgetv0, titref, trvec 222c 223c %------------------% 224c | Scalar Arguments | 225c %------------------% 226c 227 character bmat*1, which*2 228 integer ido, info, ishift, iupd, mode, ldh, ldq, ldv, mxiter, 229 & n, nev, np 230 Real 231 & tol 232c 233c %-----------------% 234c | Array Arguments | 235c %-----------------% 236c 237 integer ipntr(13) 238 Real 239 & bounds(nev+np), h(ldh,nev+np), q(ldq,nev+np), resid(n), 240 & ritzi(nev+np), ritzr(nev+np), v(ldv,nev+np), 241 & workd(3*n), workl( (nev+np)*(nev+np+3) ) 242c 243c %------------% 244c | Parameters | 245c %------------% 246c 247 Real 248 & one, zero 249 parameter (one = 1.0E+0, zero = 0.0E+0) 250c 251c %---------------% 252c | Local Scalars | 253c %---------------% 254c 255 character wprime*2 256 logical cnorm, getv0, initv, update, ushift 257 integer ierr, iter, j, kplusp, msglvl, nconv, nevbef, nev0, 258 & np0, nptemp, numcnv 259 Real 260 & rnorm, temp, eps23 261c 262c %-----------------------% 263c | Local array arguments | 264c %-----------------------% 265c 266 integer kp(4) 267 save 268c 269c %----------------------% 270c | External Subroutines | 271c %----------------------% 272c 273 external scopy, sgetv0, snaitr, snconv, sneigh, sngets, snapps, 274 & svout, ivout, second 275c 276c %--------------------% 277c | External Functions | 278c %--------------------% 279c 280 Real 281 & sdot, snrm2, slapy2, slamch 282 external sdot, snrm2, slapy2, slamch 283c 284c %---------------------% 285c | Intrinsic Functions | 286c %---------------------% 287c 288 intrinsic min, max, abs, sqrt 289c 290c %-----------------------% 291c | Executable Statements | 292c %-----------------------% 293c 294 if (ido .eq. 0) then 295c 296 call second (t0) 297c 298 msglvl = mnaup2 299c 300c %-------------------------------------% 301c | Get the machine dependent constant. | 302c %-------------------------------------% 303c 304 eps23 = slamch('Epsilon-Machine') 305 eps23 = eps23**(2.0E+0 / 3.0E+0) 306c 307 nev0 = nev 308 np0 = np 309c 310c %-------------------------------------% 311c | kplusp is the bound on the largest | 312c | Lanczos factorization built. | 313c | nconv is the current number of | 314c | "converged" eigenvlues. | 315c | iter is the counter on the current | 316c | iteration step. | 317c %-------------------------------------% 318c 319 kplusp = nev + np 320 nconv = 0 321 iter = 0 322c 323c %---------------------------------------% 324c | Set flags for computing the first NEV | 325c | steps of the Arnoldi factorization. | 326c %---------------------------------------% 327c 328 getv0 = .true. 329 update = .false. 330 ushift = .false. 331 cnorm = .false. 332c 333 if (info .ne. 0) then 334c 335c %--------------------------------------------% 336c | User provides the initial residual vector. | 337c %--------------------------------------------% 338c 339 initv = .true. 340 info = 0 341 else 342 initv = .false. 343 end if 344 end if 345c 346c %---------------------------------------------% 347c | Get a possibly random starting vector and | 348c | force it into the range of the operator OP. | 349c %---------------------------------------------% 350c 351 10 continue 352c 353 if (getv0) then 354 call sgetv0 (ido, bmat, 1, initv, n, 1, v, ldv, resid, rnorm, 355 & ipntr, workd, info) 356c 357 if (ido .ne. 99) go to 9000 358c 359 if (rnorm .eq. zero) then 360c 361c %-----------------------------------------% 362c | The initial vector is zero. Error exit. | 363c %-----------------------------------------% 364c 365 info = -9 366 go to 1100 367 end if 368 getv0 = .false. 369 ido = 0 370 end if 371c 372c %-----------------------------------% 373c | Back from reverse communication : | 374c | continue with update step | 375c %-----------------------------------% 376c 377 if (update) go to 20 378c 379c %-------------------------------------------% 380c | Back from computing user specified shifts | 381c %-------------------------------------------% 382c 383 if (ushift) go to 50 384c 385c %-------------------------------------% 386c | Back from computing residual norm | 387c | at the end of the current iteration | 388c %-------------------------------------% 389c 390 if (cnorm) go to 100 391c 392c %----------------------------------------------------------% 393c | Compute the first NEV steps of the Arnoldi factorization | 394c %----------------------------------------------------------% 395c 396 call snaitr (ido, bmat, n, 0, nev, mode, resid, rnorm, v, ldv, 397 & h, ldh, ipntr, workd, info) 398c 399c %---------------------------------------------------% 400c | ido .ne. 99 implies use of reverse communication | 401c | to compute operations involving OP and possibly B | 402c %---------------------------------------------------% 403c 404 if (ido .ne. 99) go to 9000 405c 406 if (info .gt. 0) then 407 np = info 408 mxiter = iter 409 info = -9999 410 go to 1200 411 end if 412c 413c %--------------------------------------------------------------% 414c | | 415c | M A I N ARNOLDI I T E R A T I O N L O O P | 416c | Each iteration implicitly restarts the Arnoldi | 417c | factorization in place. | 418c | | 419c %--------------------------------------------------------------% 420c 421 1000 continue 422c 423 iter = iter + 1 424c 425 if (msglvl .gt. 0) then 426 call ivout (logfil, 1, iter, ndigit, 427 & '_naup2: **** Start of major iteration number ****') 428 end if 429c 430c %-----------------------------------------------------------% 431c | Compute NP additional steps of the Arnoldi factorization. | 432c | Adjust NP since NEV might have been updated by last call | 433c | to the shift application routine snapps. | 434c %-----------------------------------------------------------% 435c 436 np = kplusp - nev 437c 438 if (msglvl .gt. 1) then 439 call ivout (logfil, 1, nev, ndigit, 440 & '_naup2: The length of the current Arnoldi factorization') 441 call ivout (logfil, 1, np, ndigit, 442 & '_naup2: Extend the Arnoldi factorization by') 443 end if 444c 445c %-----------------------------------------------------------% 446c | Compute NP additional steps of the Arnoldi factorization. | 447c %-----------------------------------------------------------% 448c 449 ido = 0 450 20 continue 451 update = .true. 452c 453 call snaitr (ido, bmat, n, nev, np, mode, resid, rnorm, v, ldv, 454 & h, ldh, ipntr, workd, info) 455c 456c %---------------------------------------------------% 457c | ido .ne. 99 implies use of reverse communication | 458c | to compute operations involving OP and possibly B | 459c %---------------------------------------------------% 460c 461 if (ido .ne. 99) go to 9000 462c 463 if (info .gt. 0) then 464 np = info 465 mxiter = iter 466 info = -9999 467 go to 1200 468 end if 469 update = .false. 470c 471 if (msglvl .gt. 1) then 472 call svout (logfil, 1, rnorm, ndigit, 473 & '_naup2: Corresponding B-norm of the residual') 474 end if 475c 476c %--------------------------------------------------------% 477c | Compute the eigenvalues and corresponding error bounds | 478c | of the current upper Hessenberg matrix. | 479c %--------------------------------------------------------% 480c 481 call sneigh (rnorm, kplusp, h, ldh, ritzr, ritzi, bounds, 482 & q, ldq, workl, ierr) 483c 484 if (ierr .ne. 0) then 485 info = -8 486 go to 1200 487 end if 488c 489c %----------------------------------------------------% 490c | Make a copy of eigenvalues and corresponding error | 491c | bounds obtained from sneigh. | 492c %----------------------------------------------------% 493c 494 call scopy(kplusp, ritzr, 1, workl(kplusp**2+1), 1) 495 call scopy(kplusp, ritzi, 1, workl(kplusp**2+kplusp+1), 1) 496 call scopy(kplusp, bounds, 1, workl(kplusp**2+2*kplusp+1), 1) 497c 498c %---------------------------------------------------% 499c | Select the wanted Ritz values and their bounds | 500c | to be used in the convergence test. | 501c | The wanted part of the spectrum and corresponding | 502c | error bounds are in the last NEV loc. of RITZR, | 503c | RITZI and BOUNDS respectively. The variables NEV | 504c | and NP may be updated if the NEV-th wanted Ritz | 505c | value has a non zero imaginary part. In this case | 506c | NEV is increased by one and NP decreased by one. | 507c | NOTE: The last two arguments of sngets are no | 508c | longer used as of version 2.1. | 509c %---------------------------------------------------% 510c 511 nev = nev0 512 np = np0 513 numcnv = nev 514 call sngets (ishift, which, nev, np, ritzr, ritzi, 515 & bounds, workl, workl(np+1)) 516 if (nev .eq. nev0+1) numcnv = nev0+1 517c 518c %-------------------% 519c | Convergence test. | 520c %-------------------% 521c 522 call scopy (nev, bounds(np+1), 1, workl(2*np+1), 1) 523 call snconv (nev, ritzr(np+1), ritzi(np+1), workl(2*np+1), 524 & tol, nconv) 525c 526 if (msglvl .gt. 2) then 527 kp(1) = nev 528 kp(2) = np 529 kp(3) = numcnv 530 kp(4) = nconv 531 call ivout (logfil, 4, kp, ndigit, 532 & '_naup2: NEV, NP, NUMCNV, NCONV are') 533 call svout (logfil, kplusp, ritzr, ndigit, 534 & '_naup2: Real part of the eigenvalues of H') 535 call svout (logfil, kplusp, ritzi, ndigit, 536 & '_naup2: Imaginary part of the eigenvalues of H') 537 call svout (logfil, kplusp, bounds, ndigit, 538 & '_naup2: Ritz estimates of the current NCV Ritz values') 539 end if 540c 541c %---------------------------------------------------------% 542c | Count the number of unwanted Ritz values that have zero | 543c | Ritz estimates. If any Ritz estimates are equal to zero | 544c | then a leading block of H of order equal to at least | 545c | the number of Ritz values with zero Ritz estimates has | 546c | split off. None of these Ritz values may be removed by | 547c | shifting. Decrease NP the number of shifts to apply. If | 548c | no shifts may be applied, then prepare to exit | 549c %---------------------------------------------------------% 550c 551 nptemp = np 552 do 30 j=1, nptemp 553 if (bounds(j) .eq. zero) then 554 np = np - 1 555 nev = nev + 1 556 end if 557 30 continue 558c 559 if ( (nconv .ge. numcnv) .or. 560 & (iter .gt. mxiter) .or. 561 & (np .eq. 0) ) then 562c 563 if (msglvl .gt. 4) then 564 call svout(logfil, kplusp, workl(kplusp**2+1), ndigit, 565 & '_naup2: Real part of the eig computed by _neigh:') 566 call svout(logfil, kplusp, workl(kplusp**2+kplusp+1), 567 & ndigit, 568 & '_naup2: Imag part of the eig computed by _neigh:') 569 call svout(logfil, kplusp, workl(kplusp**2+kplusp*2+1), 570 & ndigit, 571 & '_naup2: Ritz eistmates computed by _neigh:') 572 end if 573c 574c %------------------------------------------------% 575c | Prepare to exit. Put the converged Ritz values | 576c | and corresponding bounds in RITZ(1:NCONV) and | 577c | BOUNDS(1:NCONV) respectively. Then sort. Be | 578c | careful when NCONV > NP | 579c %------------------------------------------------% 580c 581c %------------------------------------------% 582c | Use h( 3,1 ) as storage to communicate | 583c | rnorm to _neupd if needed | 584c %------------------------------------------% 585 586 h(3,1) = rnorm 587c 588c %----------------------------------------------% 589c | To be consistent with sngets, we first do a | 590c | pre-processing sort in order to keep complex | 591c | conjugate pairs together. This is similar | 592c | to the pre-processing sort used in sngets | 593c | except that the sort is done in the opposite | 594c | order. | 595c %----------------------------------------------% 596c 597 if (which .eq. 'LM') wprime = 'SR' 598 if (which .eq. 'SM') wprime = 'LR' 599 if (which .eq. 'LR') wprime = 'SM' 600 if (which .eq. 'SR') wprime = 'LM' 601 if (which .eq. 'LI') wprime = 'SM' 602 if (which .eq. 'SI') wprime = 'LM' 603c 604 call ssortc (wprime, .true., kplusp, ritzr, ritzi, bounds) 605c 606c %----------------------------------------------% 607c | Now sort Ritz values so that converged Ritz | 608c | values appear within the first NEV locations | 609c | of ritzr, ritzi and bounds, and the most | 610c | desired one appears at the front. | 611c %----------------------------------------------% 612c 613 if (which .eq. 'LM') wprime = 'SM' 614 if (which .eq. 'SM') wprime = 'LM' 615 if (which .eq. 'LR') wprime = 'SR' 616 if (which .eq. 'SR') wprime = 'LR' 617 if (which .eq. 'LI') wprime = 'SI' 618 if (which .eq. 'SI') wprime = 'LI' 619c 620 call ssortc(wprime, .true., kplusp, ritzr, ritzi, bounds) 621c 622c %--------------------------------------------------% 623c | Scale the Ritz estimate of each Ritz value | 624c | by 1 / max(eps23,magnitude of the Ritz value). | 625c %--------------------------------------------------% 626c 627 do 35 j = 1, nev0 628 temp = max(eps23,slapy2(ritzr(j), 629 & ritzi(j))) 630 bounds(j) = bounds(j)/temp 631 35 continue 632c 633c %----------------------------------------------------% 634c | Sort the Ritz values according to the scaled Ritz | 635c | esitmates. This will push all the converged ones | 636c | towards the front of ritzr, ritzi, bounds | 637c | (in the case when NCONV < NEV.) | 638c %----------------------------------------------------% 639c 640 wprime = 'LR' 641 call ssortc(wprime, .true., nev0, bounds, ritzr, ritzi) 642c 643c %----------------------------------------------% 644c | Scale the Ritz estimate back to its original | 645c | value. | 646c %----------------------------------------------% 647c 648 do 40 j = 1, nev0 649 temp = max(eps23, slapy2(ritzr(j), 650 & ritzi(j))) 651 bounds(j) = bounds(j)*temp 652 40 continue 653c 654c %------------------------------------------------% 655c | Sort the converged Ritz values again so that | 656c | the "threshold" value appears at the front of | 657c | ritzr, ritzi and bound. | 658c %------------------------------------------------% 659c 660 call ssortc(which, .true., nconv, ritzr, ritzi, bounds) 661c 662 if (msglvl .gt. 1) then 663 call svout (logfil, kplusp, ritzr, ndigit, 664 & '_naup2: Sorted real part of the eigenvalues') 665 call svout (logfil, kplusp, ritzi, ndigit, 666 & '_naup2: Sorted imaginary part of the eigenvalues') 667 call svout (logfil, kplusp, bounds, ndigit, 668 & '_naup2: Sorted ritz estimates.') 669 end if 670c 671c %------------------------------------% 672c | Max iterations have been exceeded. | 673c %------------------------------------% 674c 675 if (iter .gt. mxiter .and. nconv .lt. numcnv) info = 1 676c 677c %---------------------% 678c | No shifts to apply. | 679c %---------------------% 680c 681 if (np .eq. 0 .and. nconv .lt. numcnv) info = 2 682c 683 np = nconv 684 go to 1100 685c 686 else if ( (nconv .lt. numcnv) .and. (ishift .eq. 1) ) then 687c 688c %-------------------------------------------------% 689c | Do not have all the requested eigenvalues yet. | 690c | To prevent possible stagnation, adjust the size | 691c | of NEV. | 692c %-------------------------------------------------% 693c 694 nevbef = nev 695 nev = nev + min(nconv, np/2) 696 if (nev .eq. 1 .and. kplusp .ge. 6) then 697 nev = kplusp / 2 698 else if (nev .eq. 1 .and. kplusp .gt. 3) then 699 nev = 2 700 end if 701 np = kplusp - nev 702c 703c %---------------------------------------% 704c | If the size of NEV was just increased | 705c | resort the eigenvalues. | 706c %---------------------------------------% 707c 708 if (nevbef .lt. nev) 709 & call sngets (ishift, which, nev, np, ritzr, ritzi, 710 & bounds, workl, workl(np+1)) 711c 712 end if 713c 714 if (msglvl .gt. 0) then 715 call ivout (logfil, 1, nconv, ndigit, 716 & '_naup2: no. of "converged" Ritz values at this iter.') 717 if (msglvl .gt. 1) then 718 kp(1) = nev 719 kp(2) = np 720 call ivout (logfil, 2, kp, ndigit, 721 & '_naup2: NEV and NP are') 722 call svout (logfil, nev, ritzr(np+1), ndigit, 723 & '_naup2: "wanted" Ritz values -- real part') 724 call svout (logfil, nev, ritzi(np+1), ndigit, 725 & '_naup2: "wanted" Ritz values -- imag part') 726 call svout (logfil, nev, bounds(np+1), ndigit, 727 & '_naup2: Ritz estimates of the "wanted" values ') 728 end if 729 end if 730c 731 if (ishift .eq. 0) then 732c 733c %-------------------------------------------------------% 734c | User specified shifts: reverse comminucation to | 735c | compute the shifts. They are returned in the first | 736c | 2*NP locations of WORKL. | 737c %-------------------------------------------------------% 738c 739 ushift = .true. 740 ido = 3 741 go to 9000 742 end if 743c 744 50 continue 745c 746c %------------------------------------% 747c | Back from reverse communication; | 748c | User specified shifts are returned | 749c | in WORKL(1:2*NP) | 750c %------------------------------------% 751c 752 ushift = .false. 753c 754 if ( ishift .eq. 0 ) then 755c 756c %----------------------------------% 757c | Move the NP shifts from WORKL to | 758c | RITZR, RITZI to free up WORKL | 759c | for non-exact shift case. | 760c %----------------------------------% 761c 762 call scopy (np, workl, 1, ritzr, 1) 763 call scopy (np, workl(np+1), 1, ritzi, 1) 764 end if 765c 766 if (msglvl .gt. 2) then 767 call ivout (logfil, 1, np, ndigit, 768 & '_naup2: The number of shifts to apply ') 769 call svout (logfil, np, ritzr, ndigit, 770 & '_naup2: Real part of the shifts') 771 call svout (logfil, np, ritzi, ndigit, 772 & '_naup2: Imaginary part of the shifts') 773 if ( ishift .eq. 1 ) 774 & call svout (logfil, np, bounds, ndigit, 775 & '_naup2: Ritz estimates of the shifts') 776 end if 777c 778c %---------------------------------------------------------% 779c | Apply the NP implicit shifts by QR bulge chasing. | 780c | Each shift is applied to the whole upper Hessenberg | 781c | matrix H. | 782c | The first 2*N locations of WORKD are used as workspace. | 783c %---------------------------------------------------------% 784c 785 call snapps (n, nev, np, ritzr, ritzi, v, ldv, 786 & h, ldh, resid, q, ldq, workl, workd) 787c 788c %---------------------------------------------% 789c | Compute the B-norm of the updated residual. | 790c | Keep B*RESID in WORKD(1:N) to be used in | 791c | the first step of the next call to snaitr. | 792c %---------------------------------------------% 793c 794 cnorm = .true. 795 call second (t2) 796 if (bmat .eq. 'G') then 797 nbx = nbx + 1 798 call scopy (n, resid, 1, workd(n+1), 1) 799 ipntr(1) = n + 1 800 ipntr(2) = 1 801 ido = 2 802c 803c %----------------------------------% 804c | Exit in order to compute B*RESID | 805c %----------------------------------% 806c 807 go to 9000 808 else if (bmat .eq. 'I') then 809 call scopy (n, resid, 1, workd, 1) 810 end if 811c 812 100 continue 813c 814c %----------------------------------% 815c | Back from reverse communication; | 816c | WORKD(1:N) := B*RESID | 817c %----------------------------------% 818c 819 if (bmat .eq. 'G') then 820 call second (t3) 821 tmvbx = tmvbx + (t3 - t2) 822 end if 823c 824 if (bmat .eq. 'G') then 825 rnorm = sdot (n, resid, 1, workd, 1) 826 rnorm = sqrt(abs(rnorm)) 827 else if (bmat .eq. 'I') then 828 rnorm = snrm2(n, resid, 1) 829 end if 830 cnorm = .false. 831c 832 if (msglvl .gt. 2) then 833 call svout (logfil, 1, rnorm, ndigit, 834 & '_naup2: B-norm of residual for compressed factorization') 835 call smout (logfil, nev, nev, h, ldh, ndigit, 836 & '_naup2: Compressed upper Hessenberg matrix H') 837 end if 838c 839 go to 1000 840c 841c %---------------------------------------------------------------% 842c | | 843c | E N D O F M A I N I T E R A T I O N L O O P | 844c | | 845c %---------------------------------------------------------------% 846c 847 1100 continue 848c 849 mxiter = iter 850 nev = numcnv 851c 852 1200 continue 853 ido = 99 854c 855c %------------% 856c | Error Exit | 857c %------------% 858c 859 call second (t1) 860 tnaup2 = t1 - t0 861c 862 9000 continue 863c 864c %---------------% 865c | End of snaup2 | 866c %---------------% 867c 868 return 869 end 870