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