Télécharger dsaup2.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales