1c\BeginDoc 2c 3c\Name: dnaup2 4c 5c\Description: 6c Intermediate level interface called by dnaupd. 7c 8c\Usage: 9c call dnaup2 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 dnaupd. 17c MODE, ISHIFT, MXITER: see the definition of IPARAM in dnaupd. 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 Double precision 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 Double precision (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, Double precision 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 Double precision array of length NEV+NP. (OUTPUT) 62c BOUNDS(1:NEV) contain the error bounds corresponding to 63c the computed Ritz values. 64c 65c Q Double precision (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 Double precision 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 dneigh. 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 Double precision 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 dgetv0 ARPACK initial vector generation routine. 139c dnaitr ARPACK Arnoldi factorization routine. 140c dnapps ARPACK application of implicit shifts routine. 141c dnconv ARPACK convergence of Ritz values routine. 142c dneigh ARPACK compute Ritz values and error bounds routine. 143c dngets ARPACK reorder Ritz values and error bounds routine. 144c dsortc ARPACK sorting routine. 145c ivout ARPACK utility routine that prints integers. 146c second ARPACK utility routine for timing. 147c dmout ARPACK utility routine that prints matrices 148c dvout ARPACK utility routine that prints vectors. 149c dlamch LAPACK routine that determines machine constants. 150c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully. 151c dcopy Level 1 BLAS that copies one vector to another . 152c ddot Level 1 BLAS that computes the scalar product of two vectors. 153c dnrm2 Level 1 BLAS that computes the norm of a vector. 154c dswap 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 dnaup2 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 183c 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 200c %--------------------------------% 201c | See stat.doc for documentation | 202c %--------------------------------% 203c 204c\SCCS Information: @(#) 205c FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 206c 207 real t0, t1, t2, t3, t4, t5 208c save t0, t1, t2, t3, t4, t5 209c 210 integer nopx, nbx, nrorth, nitref, nrstrt 211 real tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, 212 & tnaupd, tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, 213 & tcaupd, tcaup2, tcaitr, tceigh, tcgets, tcapps, tcconv, 214 & tmvopx, tmvbx, tgetv0, titref, trvec 215 common /timing/ 216 & nopx, nbx, nrorth, nitref, nrstrt, 217 & tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, 218 & tnaupd, tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, 219 & tcaupd, tcaup2, tcaitr, tceigh, tcgets, tcapps, tcconv, 220 & tmvopx, tmvbx, tgetv0, titref, trvec 221c %------------------% 222c | Scalar Arguments | 223c %------------------% 224c 225 character bmat*1, which*2 226 integer ido, info, ishift, iupd, mode, ldh, ldq, ldv, mxiter, 227 & n, nev, np 228 Double precision 229 & tol 230c 231c %-----------------% 232c | Array Arguments | 233c %-----------------% 234c 235 integer ipntr(13) 236 Double precision 237 & bounds(nev+np), h(ldh,nev+np), q(ldq,nev+np), resid(n), 238 & ritzi(nev+np), ritzr(nev+np), v(ldv,nev+np), 239 & workd(3*n), workl( (nev+np)*(nev+np+3) ) 240c 241c %------------% 242c | Parameters | 243c %------------% 244c 245 Double precision 246 & one, zero 247 parameter (one = 1.0D+0, zero = 0.0D+0) 248c 249c %---------------% 250c | Local Scalars | 251c %---------------% 252c 253 character wprime*2 254 logical cnorm, getv0, initv, update, ushift 255 integer ierr, iter, j, kplusp, msglvl, nconv, nevbef, nev0, 256 & np0, nptemp, numcnv 257 Double precision 258 & rnorm, temp, eps23 259c 260c %-----------------------% 261c | Local array arguments | 262c %-----------------------% 263c 264 integer kp(4) 265 save 266 267c 268c %----------------------% 269c | External Subroutines | 270c %----------------------% 271c 272 external dcopy, dgetv0, dnaitr, dnconv, dneigh, dngets, dnapps, 273 & dvout, ivout, second 274c 275c %--------------------% 276c | External Functions | 277c %--------------------% 278c 279 Double precision 280 & ddot, dnrm2, dlapy2, dlamch 281 external ddot, dnrm2, dlapy2, dlamch 282c 283c %---------------------% 284c | Intrinsic Functions | 285c %---------------------% 286c 287 intrinsic min, max, abs, sqrt 288c 289c %-----------------------% 290c | Executable Statements | 291c %-----------------------% 292c 293 if (ido .eq. 0) then 294c 295 call second (t0) 296c 297 msglvl = mnaup2 298c 299c %-------------------------------------% 300c | Get the machine dependent constant. | 301c %-------------------------------------% 302c 303 eps23 = dlamch('Epsilon-Machine') 304 eps23 = eps23**(2.0D+0 / 3.0D+0) 305c 306 nev0 = nev 307 np0 = np 308c 309c %-------------------------------------% 310c | kplusp is the bound on the largest | 311c | Lanczos factorization built. | 312c | nconv is the current number of | 313c | "converged" eigenvlues. | 314c | iter is the counter on the current | 315c | iteration step. | 316c %-------------------------------------% 317c 318 kplusp = nev + np 319 nconv = 0 320 iter = 0 321c 322c %---------------------------------------% 323c | Set flags for computing the first NEV | 324c | steps of the Arnoldi factorization. | 325c %---------------------------------------% 326c 327 getv0 = .true. 328 update = .false. 329 ushift = .false. 330 cnorm = .false. 331c 332 if (info .ne. 0) then 333c 334c %--------------------------------------------% 335c | User provides the initial residual vector. | 336c %--------------------------------------------% 337c 338 initv = .true. 339 info = 0 340 else 341 initv = .false. 342 end if 343 end if 344c 345c %---------------------------------------------% 346c | Get a possibly random starting vector and | 347c | force it into the range of the operator OP. | 348c %---------------------------------------------% 349c 350 10 continue 351c 352 if (getv0) then 353 call dgetv0 (ido, bmat, 1, initv, n, 1, v, ldv, resid, rnorm, 354 & ipntr, workd, info) 355c 356 if (ido .ne. 99) go to 9000 357c 358 if (rnorm .eq. zero) then 359c 360c %-----------------------------------------% 361c | The initial vector is zero. Error exit. | 362c %-----------------------------------------% 363c 364 info = -9 365 go to 1100 366 end if 367 getv0 = .false. 368 ido = 0 369 end if 370c 371c %-----------------------------------% 372c | Back from reverse communication : | 373c | continue with update step | 374c %-----------------------------------% 375c 376 if (update) go to 20 377c 378c %-------------------------------------------% 379c | Back from computing user specified shifts | 380c %-------------------------------------------% 381c 382 if (ushift) go to 50 383c 384c %-------------------------------------% 385c | Back from computing residual norm | 386c | at the end of the current iteration | 387c %-------------------------------------% 388c 389 if (cnorm) go to 100 390c 391c %----------------------------------------------------------% 392c | Compute the first NEV steps of the Arnoldi factorization | 393c %----------------------------------------------------------% 394c 395 call dnaitr (ido, bmat, n, 0, nev, mode, resid, rnorm, v, ldv, 396 & h, ldh, ipntr, workd, info) 397c 398c %---------------------------------------------------% 399c | ido .ne. 99 implies use of reverse communication | 400c | to compute operations involving OP and possibly B | 401c %---------------------------------------------------% 402c 403 if (ido .ne. 99) go to 9000 404c 405 if (info .gt. 0) then 406 np = info 407 mxiter = iter 408 info = -9999 409 go to 1200 410 end if 411c 412c %--------------------------------------------------------------% 413c | | 414c | M A I N ARNOLDI I T E R A T I O N L O O P | 415c | Each iteration implicitly restarts the Arnoldi | 416c | factorization in place. | 417c | | 418c %--------------------------------------------------------------% 419c 420 1000 continue 421c 422 iter = iter + 1 423c 424 if (msglvl .gt. 0) then 425 call ivout (logfil, 1, iter, ndigit, 426 & '_naup2: **** Start of major iteration number ****') 427 end if 428c 429c %-----------------------------------------------------------% 430c | Compute NP additional steps of the Arnoldi factorization. | 431c | Adjust NP since NEV might have been updated by last call | 432c | to the shift application routine dnapps. | 433c %-----------------------------------------------------------% 434c 435 np = kplusp - nev 436c 437 if (msglvl .gt. 1) then 438 call ivout (logfil, 1, nev, ndigit, 439 & '_naup2: The length of the current Arnoldi factorization') 440 call ivout (logfil, 1, np, ndigit, 441 & '_naup2: Extend the Arnoldi factorization by') 442 end if 443c 444c %-----------------------------------------------------------% 445c | Compute NP additional steps of the Arnoldi factorization. | 446c %-----------------------------------------------------------% 447c 448 ido = 0 449 20 continue 450 update = .true. 451c 452 call dnaitr (ido, bmat, n, nev, np, mode, resid, rnorm, v, ldv, 453 & h, ldh, ipntr, workd, info) 454c 455c %---------------------------------------------------% 456c | ido .ne. 99 implies use of reverse communication | 457c | to compute operations involving OP and possibly B | 458c %---------------------------------------------------% 459c 460 if (ido .ne. 99) go to 9000 461c 462 if (info .gt. 0) then 463 np = info 464 mxiter = iter 465 info = -9999 466 go to 1200 467 end if 468 update = .false. 469c 470 if (msglvl .gt. 1) then 471 call dvout (logfil, 1, rnorm, ndigit, 472 & '_naup2: Corresponding B-norm of the residual') 473 end if 474c 475c %--------------------------------------------------------% 476c | Compute the eigenvalues and corresponding error bounds | 477c | of the current upper Hessenberg matrix. | 478c %--------------------------------------------------------% 479c 480 call dneigh (rnorm, kplusp, h, ldh, ritzr, ritzi, bounds, 481 & q, ldq, workl, ierr) 482c 483 if (ierr .ne. 0) then 484 info = -8 485 go to 1200 486 end if 487c 488c %----------------------------------------------------% 489c | Make a copy of eigenvalues and corresponding error | 490c | bounds obtained from dneigh. | 491c %----------------------------------------------------% 492c 493 call dcopy(kplusp, ritzr, 1, workl(kplusp**2+1), 1) 494 call dcopy(kplusp, ritzi, 1, workl(kplusp**2+kplusp+1), 1) 495 call dcopy(kplusp, bounds, 1, workl(kplusp**2+2*kplusp+1), 1) 496c 497c %---------------------------------------------------% 498c | Select the wanted Ritz values and their bounds | 499c | to be used in the convergence test. | 500c | The wanted part of the spectrum and corresponding | 501c | error bounds are in the last NEV loc. of RITZR, | 502c | RITZI and BOUNDS respectively. The variables NEV | 503c | and NP may be updated if the NEV-th wanted Ritz | 504c | value has a non zero imaginary part. In this case | 505c | NEV is increased by one and NP decreased by one. | 506c | NOTE: The last two arguments of dngets are no | 507c | longer used as of version 2.1. | 508c %---------------------------------------------------% 509c 510 nev = nev0 511 np = np0 512 numcnv = nev 513 call dngets (ishift, which, nev, np, ritzr, ritzi, 514 & bounds, workl, workl(np+1)) 515 if (nev .eq. nev0+1) numcnv = nev0+1 516c 517c %-------------------% 518c | Convergence test. | 519c %-------------------% 520c 521 call dcopy (nev, bounds(np+1), 1, workl(2*np+1), 1) 522 call dnconv (nev, ritzr(np+1), ritzi(np+1), workl(2*np+1), 523 & tol, nconv) 524c 525 if (msglvl .gt. 2) then 526 kp(1) = nev 527 kp(2) = np 528 kp(3) = numcnv 529 kp(4) = nconv 530 call ivout (logfil, 4, kp, ndigit, 531 & '_naup2: NEV, NP, NUMCNV, NCONV are') 532 call dvout (logfil, kplusp, ritzr, ndigit, 533 & '_naup2: Real part of the eigenvalues of H') 534 call dvout (logfil, kplusp, ritzi, ndigit, 535 & '_naup2: Imaginary part of the eigenvalues of H') 536 call dvout (logfil, kplusp, bounds, ndigit, 537 & '_naup2: Ritz estimates of the current NCV Ritz values') 538 end if 539c 540c %---------------------------------------------------------% 541c | Count the number of unwanted Ritz values that have zero | 542c | Ritz estimates. If any Ritz estimates are equal to zero | 543c | then a leading block of H of order equal to at least | 544c | the number of Ritz values with zero Ritz estimates has | 545c | split off. None of these Ritz values may be removed by | 546c | shifting. Decrease NP the number of shifts to apply. If | 547c | no shifts may be applied, then prepare to exit | 548c %---------------------------------------------------------% 549c 550 nptemp = np 551 do 30 j=1, nptemp 552 if (bounds(j) .eq. zero) then 553 np = np - 1 554 nev = nev + 1 555 end if 556 30 continue 557c 558 if ( (nconv .ge. numcnv) .or. 559 & (iter .gt. mxiter) .or. 560 & (np .eq. 0) ) then 561c 562 if (msglvl .gt. 4) then 563 call dvout(logfil, kplusp, workl(kplusp**2+1), ndigit, 564 & '_naup2: Real part of the eig computed by _neigh:') 565 call dvout(logfil, kplusp, workl(kplusp**2+kplusp+1), 566 & ndigit, 567 & '_naup2: Imag part of the eig computed by _neigh:') 568 call dvout(logfil, kplusp, workl(kplusp**2+kplusp*2+1), 569 & ndigit, 570 & '_naup2: Ritz eistmates computed by _neigh:') 571 end if 572c 573c %------------------------------------------------% 574c | Prepare to exit. Put the converged Ritz values | 575c | and corresponding bounds in RITZ(1:NCONV) and | 576c | BOUNDS(1:NCONV) respectively. Then sort. Be | 577c | careful when NCONV > NP | 578c %------------------------------------------------% 579c 580c %------------------------------------------% 581c | Use h( 3,1 ) as storage to communicate | 582c | rnorm to _neupd if needed | 583c %------------------------------------------% 584 585 h(3,1) = rnorm 586c 587c %----------------------------------------------% 588c | To be consistent with dngets, we first do a | 589c | pre-processing sort in order to keep complex | 590c | conjugate pairs together. This is similar | 591c | to the pre-processing sort used in dngets | 592c | except that the sort is done in the opposite | 593c | order. | 594c %----------------------------------------------% 595c 596 if (which .eq. 'LM') wprime = 'SR' 597 if (which .eq. 'SM') wprime = 'LR' 598 if (which .eq. 'LR') wprime = 'SM' 599 if (which .eq. 'SR') wprime = 'LM' 600 if (which .eq. 'LI') wprime = 'SM' 601 if (which .eq. 'SI') wprime = 'LM' 602c 603 call dsortc (wprime, .true., kplusp, ritzr, ritzi, bounds) 604c 605c %----------------------------------------------% 606c | Now sort Ritz values so that converged Ritz | 607c | values appear within the first NEV locations | 608c | of ritzr, ritzi and bounds, and the most | 609c | desired one appears at the front. | 610c %----------------------------------------------% 611c 612 if (which .eq. 'LM') wprime = 'SM' 613 if (which .eq. 'SM') wprime = 'LM' 614 if (which .eq. 'LR') wprime = 'SR' 615 if (which .eq. 'SR') wprime = 'LR' 616 if (which .eq. 'LI') wprime = 'SI' 617 if (which .eq. 'SI') wprime = 'LI' 618c 619 call dsortc(wprime, .true., kplusp, ritzr, ritzi, bounds) 620c 621c %--------------------------------------------------% 622c | Scale the Ritz estimate of each Ritz value | 623c | by 1 / max(eps23,magnitude of the Ritz value). | 624c %--------------------------------------------------% 625c 626 do 35 j = 1, nev0 627 temp = max(eps23,dlapy2(ritzr(j), 628 & ritzi(j))) 629 bounds(j) = bounds(j)/temp 630 35 continue 631c 632c %----------------------------------------------------% 633c | Sort the Ritz values according to the scaled Ritz | 634c | esitmates. This will push all the converged ones | 635c | towards the front of ritzr, ritzi, bounds | 636c | (in the case when NCONV < NEV.) | 637c %----------------------------------------------------% 638c 639 wprime = 'LR' 640 call dsortc(wprime, .true., nev0, bounds, ritzr, ritzi) 641c 642c %----------------------------------------------% 643c | Scale the Ritz estimate back to its original | 644c | value. | 645c %----------------------------------------------% 646c 647 do 40 j = 1, nev0 648 temp = max(eps23, dlapy2(ritzr(j), 649 & ritzi(j))) 650 bounds(j) = bounds(j)*temp 651 40 continue 652c 653c %------------------------------------------------% 654c | Sort the converged Ritz values again so that | 655c | the "threshold" value appears at the front of | 656c | ritzr, ritzi and bound. | 657c %------------------------------------------------% 658c 659 call dsortc(which, .true., nconv, ritzr, ritzi, bounds) 660c 661 if (msglvl .gt. 1) then 662 call dvout (logfil, kplusp, ritzr, ndigit, 663 & '_naup2: Sorted real part of the eigenvalues') 664 call dvout (logfil, kplusp, ritzi, ndigit, 665 & '_naup2: Sorted imaginary part of the eigenvalues') 666 call dvout (logfil, kplusp, bounds, ndigit, 667 & '_naup2: Sorted ritz estimates.') 668 end if 669c 670c %------------------------------------% 671c | Max iterations have been exceeded. | 672c %------------------------------------% 673c 674 if (iter .gt. mxiter .and. nconv .lt. numcnv) info = 1 675c 676c %---------------------% 677c | No shifts to apply. | 678c %---------------------% 679c 680 if (np .eq. 0 .and. nconv .lt. numcnv) info = 2 681c 682 np = nconv 683 go to 1100 684c 685 else if ( (nconv .lt. numcnv) .and. (ishift .eq. 1) ) then 686c 687c %-------------------------------------------------% 688c | Do not have all the requested eigenvalues yet. | 689c | To prevent possible stagnation, adjust the size | 690c | of NEV. | 691c %-------------------------------------------------% 692c 693 nevbef = nev 694 nev = nev + min(nconv, np/2) 695 if (nev .eq. 1 .and. kplusp .ge. 6) then 696 nev = kplusp / 2 697 else if (nev .eq. 1 .and. kplusp .gt. 3) then 698 nev = 2 699 end if 700 np = kplusp - nev 701c 702c %---------------------------------------% 703c | If the size of NEV was just increased | 704c | resort the eigenvalues. | 705c %---------------------------------------% 706c 707 if (nevbef .lt. nev) 708 & call dngets (ishift, which, nev, np, ritzr, ritzi, 709 & bounds, workl, workl(np+1)) 710c 711 end if 712c 713 if (msglvl .gt. 0) then 714 call ivout (logfil, 1, nconv, ndigit, 715 & '_naup2: no. of "converged" Ritz values at this iter.') 716 if (msglvl .gt. 1) then 717 kp(1) = nev 718 kp(2) = np 719 call ivout (logfil, 2, kp, ndigit, 720 & '_naup2: NEV and NP are') 721 call dvout (logfil, nev, ritzr(np+1), ndigit, 722 & '_naup2: "wanted" Ritz values -- real part') 723 call dvout (logfil, nev, ritzi(np+1), ndigit, 724 & '_naup2: "wanted" Ritz values -- imag part') 725 call dvout (logfil, nev, bounds(np+1), ndigit, 726 & '_naup2: Ritz estimates of the "wanted" values ') 727 end if 728 end if 729c 730 if (ishift .eq. 0) then 731c 732c %-------------------------------------------------------% 733c | User specified shifts: reverse comminucation to | 734c | compute the shifts. They are returned in the first | 735c | 2*NP locations of WORKL. | 736c %-------------------------------------------------------% 737c 738 ushift = .true. 739 ido = 3 740 go to 9000 741 end if 742c 743 50 continue 744c 745c %------------------------------------% 746c | Back from reverse communication; | 747c | User specified shifts are returned | 748c | in WORKL(1:2*NP) | 749c %------------------------------------% 750c 751 ushift = .false. 752c 753 if ( ishift .eq. 0 ) then 754c 755c %----------------------------------% 756c | Move the NP shifts from WORKL to | 757c | RITZR, RITZI to free up WORKL | 758c | for non-exact shift case. | 759c %----------------------------------% 760c 761 call dcopy (np, workl, 1, ritzr, 1) 762 call dcopy (np, workl(np+1), 1, ritzi, 1) 763 end if 764c 765 if (msglvl .gt. 2) then 766 call ivout (logfil, 1, np, ndigit, 767 & '_naup2: The number of shifts to apply ') 768 call dvout (logfil, np, ritzr, ndigit, 769 & '_naup2: Real part of the shifts') 770 call dvout (logfil, np, ritzi, ndigit, 771 & '_naup2: Imaginary part of the shifts') 772 if ( ishift .eq. 1 ) 773 & call dvout (logfil, np, bounds, ndigit, 774 & '_naup2: Ritz estimates of the shifts') 775 end if 776c 777c %---------------------------------------------------------% 778c | Apply the NP implicit shifts by QR bulge chasing. | 779c | Each shift is applied to the whole upper Hessenberg | 780c | matrix H. | 781c | The first 2*N locations of WORKD are used as workspace. | 782c %---------------------------------------------------------% 783c 784 call dnapps (n, nev, np, ritzr, ritzi, v, ldv, 785 & h, ldh, resid, q, ldq, workl, workd) 786c 787c %---------------------------------------------% 788c | Compute the B-norm of the updated residual. | 789c | Keep B*RESID in WORKD(1:N) to be used in | 790c | the first step of the next call to dnaitr. | 791c %---------------------------------------------% 792c 793 cnorm = .true. 794 call second (t2) 795 if (bmat .eq. 'G') then 796 nbx = nbx + 1 797 call dcopy (n, resid, 1, workd(n+1), 1) 798 ipntr(1) = n + 1 799 ipntr(2) = 1 800 ido = 2 801c 802c %----------------------------------% 803c | Exit in order to compute B*RESID | 804c %----------------------------------% 805c 806 go to 9000 807 else if (bmat .eq. 'I') then 808 call dcopy (n, resid, 1, workd, 1) 809 end if 810c 811 100 continue 812c 813c %----------------------------------% 814c | Back from reverse communication; | 815c | WORKD(1:N) := B*RESID | 816c %----------------------------------% 817c 818 if (bmat .eq. 'G') then 819 call second (t3) 820 tmvbx = tmvbx + (t3 - t2) 821 end if 822c 823 if (bmat .eq. 'G') then 824 rnorm = ddot (n, resid, 1, workd, 1) 825 rnorm = sqrt(abs(rnorm)) 826 else if (bmat .eq. 'I') then 827 rnorm = dnrm2(n, resid, 1) 828 end if 829 cnorm = .false. 830c 831 if (msglvl .gt. 2) then 832 call dvout (logfil, 1, rnorm, ndigit, 833 & '_naup2: B-norm of residual for compressed factorization') 834 call dmout (logfil, nev, nev, h, ldh, ndigit, 835 & '_naup2: Compressed upper Hessenberg matrix H') 836 end if 837c 838 go to 1000 839c 840c %---------------------------------------------------------------% 841c | | 842c | E N D O F M A I N I T E R A T I O N L O O P | 843c | | 844c %---------------------------------------------------------------% 845c 846 1100 continue 847c 848 mxiter = iter 849 nev = numcnv 850c 851 1200 continue 852 ido = 99 853c 854c %------------% 855c | Error Exit | 856c %------------% 857c 858 call second (t1) 859 tnaup2 = t1 - t0 860c 861 9000 continue 862c 863c %---------------% 864c | End of dnaup2 | 865c %---------------% 866c 867 return 868 end 869