Télécharger dsaup2.eso

Retour à la liste

Numérotation des lignes :

  1. C DSAUP2 SOURCE BP208322 20/02/06 21:15:28 10512
  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. -> deleted by BP in 2020
  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, ITRAK, info )
  184. c
  185. c %----------------------------------------------------%
  186. c | Include files for debugging and timing information |
  187. c -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 real*8 T0,T1,T2,T3
  201. c
  202. c %-----------------%
  203. c | Array Arguments |
  204. c %-----------------%
  205. c
  206. integer ipntr(3)
  207. INTEGER ITRAK(5)
  208. REAL*8
  209. & bounds(nev+np), h(ldh,2), q(ldq,nev+np), resid(n),
  210. & ritz(nev+np), v(ldv,nev+np), workd(3*n),
  211. & workl(3*(nev+np))
  212. c
  213. c %------------%
  214. c | Parameters |
  215. c %------------%
  216. c
  217. REAL*8
  218. & one, zero
  219. parameter (one = 1.0D+0, zero = 0.0D+0)
  220. c
  221. c %---------------%
  222. c | Local Scalars |
  223. c %---------------%
  224. c
  225. character wprime*2
  226. logical cnorm, getv0, initv, update, ushift
  227. integer ierr, iter, j, kplusp, msglvl, nconv, nevbef, nev0,
  228. c & np0, nptemp, nevd2, nevm2, kp(3)
  229. & np0, nptemp, nevd2, nevm2
  230. REAL*8
  231. & rnorm, temp, eps23
  232. save cnorm, getv0, initv, update, ushift,
  233. c & iter, kplusp, msglvl, nconv, nev0, np0,
  234. & iter, kplusp, nconv, nev0, np0,
  235. & rnorm, eps23
  236. parameter (msglvl=0)
  237. c
  238. c %----------------------%
  239. c | External Subroutines |
  240. c %----------------------%
  241. c
  242. c
  243. c %--------------------%
  244. c | External Functions |
  245. c %--------------------%
  246. c
  247. REAL*8
  248. external ddot, dnrm2, dlamch
  249. c
  250. c %---------------------%
  251. **c | Intrinsic Functions |
  252. **c %---------------------%
  253. **c
  254. ** intrinsic min
  255. **c
  256. **c %-----------------------%
  257. **c | Executable Statements |
  258. c %-----------------------%
  259. c T0=0.D0
  260. c T1=0.D0
  261. c T2=0.D0
  262. c T3=0.D0
  263. NOPX =ITRAK(1)
  264. NBX =ITRAK(2)
  265. NRORTH=ITRAK(3)
  266. NITREF=ITRAK(4)
  267. NRSTRT=ITRAK(5)
  268. c
  269. if (ido .eq. 0) then
  270. c
  271. c %-------------------------------%
  272. c | Initialize timing statistics |
  273. c | & message level for debugging |
  274. c %-------------------------------%
  275. c
  276. * call arscnd (t0)
  277. c msglvl = msaup2
  278. c
  279. c %---------------------------------%
  280. c | Set machine dependent constant. |
  281. c %---------------------------------%
  282. c
  283. eps23 = dlamch('Epsilon-Machine')
  284. eps23 = eps23**(2.0D+0/3.0D+0)
  285. c
  286. c %-------------------------------------%
  287. c | nev0 and np0 are integer variables |
  288. c | hold the initial values of NEV & NP |
  289. c %-------------------------------------%
  290. c
  291. nev0 = nev
  292. np0 = np
  293. c
  294. c %-------------------------------------%
  295. c | kplusp is the bound on the largest |
  296. c | Lanczos factorization built. |
  297. c | nconv is the current number of |
  298. c | "converged" eigenvlues. |
  299. c | iter is the counter on the current |
  300. c | iteration step. |
  301. c %-------------------------------------%
  302. c
  303. kplusp = nev0 + np0
  304. nconv = 0
  305. iter = 0
  306. c
  307. c %--------------------------------------------%
  308. c | Set flags for computing the first NEV steps |
  309. c | of the Lanczos factorization. |
  310. c %--------------------------------------------%
  311. c
  312. getv0 = .true.
  313. update = .false.
  314. ushift = .false.
  315. cnorm = .false.
  316. c
  317. if (info .ne. 0) then
  318. c
  319. c %--------------------------------------------%
  320. c | User provides the initial residual vector. |
  321. c %--------------------------------------------%
  322. c
  323. initv = .true.
  324. info = 0
  325. else
  326. initv = .false.
  327. end if
  328. end if
  329. c
  330. c %---------------------------------------------%
  331. c | Get a possibly random starting vector and |
  332. c | force it into the range of the operator OP. |
  333. c %---------------------------------------------%
  334. c
  335. 10 continue
  336. c
  337. if (getv0) then
  338. call dgetv0 (ido, bmat, 1, initv, n, 1, v, ldv, resid, rnorm,
  339. & ipntr, workd, nbx,nopx, info)
  340. c
  341. if (ido .ne. 99) go to 9000
  342. c
  343. if (rnorm .eq. zero) then
  344. c
  345. c %-----------------------------------------%
  346. c | The initial vector is zero. Error exit. |
  347. c %-----------------------------------------%
  348. c
  349. info = -9
  350. go to 1200
  351. end if
  352. getv0 = .false.
  353. ido = 0
  354. end if
  355. c
  356. c %------------------------------------------------------------%
  357. c | Back from reverse communication: continue with update step |
  358. c %------------------------------------------------------------%
  359. c
  360. if (update) go to 20
  361. c
  362. c %-------------------------------------------%
  363. c | Back from computing user specified shifts |
  364. c %-------------------------------------------%
  365. c
  366. if (ushift) go to 50
  367. c
  368. c %-------------------------------------%
  369. c | Back from computing residual norm |
  370. c | at the end of the current iteration |
  371. c %-------------------------------------%
  372. c
  373. if (cnorm) go to 100
  374. c
  375. c %----------------------------------------------------------%
  376. c | Compute the first NEV steps of the Lanczos factorization |
  377. c %----------------------------------------------------------%
  378. c
  379. call dsaitr (ido, bmat, n, 0, nev0, mode, resid, rnorm, v, ldv,
  380. & h, ldh, ipntr, workd, ITRAK, info)
  381. c
  382. c %---------------------------------------------------%
  383. c | ido .ne. 99 implies use of reverse communication |
  384. c | to compute operations involving OP and possibly B |
  385. c %---------------------------------------------------%
  386. c
  387. if (ido .ne. 99) go to 9000
  388. c
  389. if (info .gt. 0) then
  390. c
  391. c %-----------------------------------------------------%
  392. c | dsaitr was unable to build an Lanczos factorization |
  393. c | of length NEV0. INFO is returned with the size of |
  394. c | the factorization built. Exit main loop. |
  395. c %-----------------------------------------------------%
  396. c
  397. np = info
  398. mxiter = iter
  399. info = -9999
  400. go to 1200
  401. end if
  402. c
  403. c %--------------------------------------------------------------%
  404. c | |
  405. c | M A I N LANCZOS I T E R A T I O N L O O P |
  406. c | Each iteration implicitly restarts the Lanczos |
  407. c | factorization in place. |
  408. c | |
  409. c %--------------------------------------------------------------%
  410. c
  411. 1000 continue
  412. c
  413. iter = iter + 1
  414. c
  415. if (msglvl .gt. 0) then
  416. call ivout ( 1, iter, ndigit,
  417. & '_saup2: **** Start of major iteration number ****')
  418. end if
  419. c if (msglvl .gt. 1) then
  420. c call ivout ( 1, nev, ndigit,
  421. c & '_saup2: The length of the current Lanczos factorization')
  422. c call ivout ( 1, np, ndigit,
  423. c & '_saup2: Extend the Lanczos factorization by')
  424. c end if
  425. c c
  426. c %------------------------------------------------------------%
  427. c | Compute NP additional steps of the Lanczos factorization. |
  428. c %------------------------------------------------------------%
  429. c
  430. ido = 0
  431. 20 continue
  432. update = .true.
  433. c
  434. call dsaitr (ido, bmat, n, nev, np, mode, resid, rnorm, v,
  435. & ldv, h, ldh, ipntr, workd, ITRAK, info)
  436. c
  437. c %---------------------------------------------------%
  438. c | ido .ne. 99 implies use of reverse communication |
  439. c | to compute operations involving OP and possibly B |
  440. c %---------------------------------------------------%
  441. c
  442. if (ido .ne. 99) go to 9000
  443. c
  444. if (info .gt. 0) then
  445. c
  446. c %-----------------------------------------------------%
  447. c | dsaitr was unable to build an Lanczos factorization |
  448. c | of length NEV0+NP0. INFO is returned with the size |
  449. c | of the factorization built. Exit main loop. |
  450. c %-----------------------------------------------------%
  451. c
  452. np = info
  453. mxiter = iter
  454. info = -9999
  455. go to 1200
  456. end if
  457. update = .false.
  458. c
  459. c if (msglvl .gt. 1) then
  460. c call dvout ( 1, rnorm, ndigit,
  461. c & '_saup2: Current B-norm of residual for factorization')
  462. c end if
  463. c
  464. c %--------------------------------------------------------%
  465. c | Compute the eigenvalues and corresponding error bounds |
  466. c | of the current symmetric tridiagonal matrix. |
  467. c %--------------------------------------------------------%
  468. c
  469. call dseigt (rnorm, kplusp, h, ldh, ritz, bounds, workl, ierr)
  470. c
  471. if (ierr .ne. 0) then
  472. info = -8
  473. go to 1200
  474. end if
  475. c
  476. c %----------------------------------------------------%
  477. c | Make a copy of eigenvalues and corresponding error |
  478. c | bounds obtained from _seigt. |
  479. c %----------------------------------------------------%
  480. c
  481. call dcopy(kplusp, ritz, 1, workl(kplusp+1), 1)
  482. call dcopy(kplusp, bounds, 1, workl(2*kplusp+1), 1)
  483. c
  484. c %---------------------------------------------------%
  485. c | Select the wanted Ritz values and their bounds |
  486. c | to be used in the convergence test. |
  487. c | The selection is based on the requested number of |
  488. c | eigenvalues instead of the current NEV and NP to |
  489. c | prevent possible misconvergence. |
  490. c | * Wanted Ritz values := RITZ(NP+1:NEV+NP) |
  491. c | * Shifts := RITZ(1:NP) := WORKL(1:NP) |
  492. c %---------------------------------------------------%
  493. c
  494. nev = nev0
  495. np = np0
  496. call dsgets (ishift, which, nev, np, ritz, bounds, workl)
  497. c
  498. c %-------------------%
  499. c | Convergence test. |
  500. c %-------------------%
  501. c
  502. call dcopy (nev, bounds(np+1), 1, workl(np+1), 1)
  503. call dsconv (nev, ritz(np+1), workl(np+1), tol, nconv)
  504. c
  505. c if (msglvl .gt. 2) then
  506. c kp(1) = nev
  507. c kp(2) = np
  508. c kp(3) = nconv
  509. c call ivout ( 3, kp, ndigit,
  510. c & '_saup2: NEV, NP, NCONV are')
  511. c call dvout ( kplusp, ritz, ndigit,
  512. c & '_saup2: The eigenvalues of H')
  513. c call dvout ( kplusp, bounds, ndigit,
  514. c & '_saup2: Ritz estimates of the current NCV Ritz values')
  515. c end if
  516. c
  517. c %---------------------------------------------------------%
  518. c | Count the number of unwanted Ritz values that have zero |
  519. c | Ritz estimates. If any Ritz estimates are equal to zero |
  520. c | then a leading block of H of order equal to at least |
  521. c | the number of Ritz values with zero Ritz estimates has |
  522. c | split off. None of these Ritz values may be removed by |
  523. c | shifting. Decrease NP the number of shifts to apply. If |
  524. c | no shifts may be applied, then prepare to exit |
  525. c %---------------------------------------------------------%
  526. c
  527. nptemp = np
  528. do 30 j=1, nptemp
  529. if (bounds(j) .eq. zero) then
  530. np = np - 1
  531. nev = nev + 1
  532. end if
  533. 30 continue
  534. c
  535. if ( (nconv .ge. nev0) .or.
  536. & (iter .gt. mxiter) .or.
  537. & (np .eq. 0) ) then
  538. c
  539. c %------------------------------------------------%
  540. c | Prepare to exit. Put the converged Ritz values |
  541. c | and corresponding bounds in RITZ(1:NCONV) and |
  542. c | BOUNDS(1:NCONV) respectively. Then sort. Be |
  543. c | careful when NCONV > NP since we don't want to |
  544. c | swap overlapping locations. |
  545. c %------------------------------------------------%
  546. c
  547. if (which .eq. 'BE') then
  548. c
  549. c %-----------------------------------------------------%
  550. c | Both ends of the spectrum are requested. |
  551. c | Sort the eigenvalues into algebraically decreasing |
  552. c | order first then swap low end of the spectrum next |
  553. c | to high end in appropriate locations. |
  554. c | NOTE: when np < floor(nev/2) be careful not to swap |
  555. c | overlapping locations. |
  556. c %-----------------------------------------------------%
  557. c
  558. wprime = 'SA'
  559. call dsortr (wprime, .true., kplusp, ritz, bounds)
  560. nevd2 = nev0 / 2
  561. nevm2 = nev0 - nevd2
  562. if ( nev .gt. 1 ) then
  563. np = kplusp - nev0
  564. call dswap ( min(nevd2,np), ritz(nevm2+1), 1,
  565. & ritz( max(kplusp-nevd2+1,kplusp-np+1) ), 1)
  566. call dswap ( min(nevd2,np), bounds(nevm2+1), 1,
  567. & bounds( max(kplusp-nevd2+1,kplusp-np+1)), 1)
  568. end if
  569. c
  570. else
  571. c
  572. c %--------------------------------------------------%
  573. c | LM, SM, LA, SA case. |
  574. c | Sort the eigenvalues of H into the an order that |
  575. c | is opposite to WHICH, and apply the resulting |
  576. c | order to BOUNDS. The eigenvalues are sorted so |
  577. c | that the wanted part are always within the first |
  578. c | NEV locations. |
  579. c %--------------------------------------------------%
  580. c
  581. if (which .eq. 'LM') wprime = 'SM'
  582. if (which .eq. 'SM') wprime = 'LM'
  583. if (which .eq. 'LA') wprime = 'SA'
  584. if (which .eq. 'SA') wprime = 'LA'
  585. c
  586. call dsortr (wprime, .true., kplusp, ritz, bounds)
  587. c
  588. end if
  589. c
  590. c %--------------------------------------------------%
  591. c | Scale the Ritz estimate of each Ritz value |
  592. c | by 1 / max(eps23,magnitude of the Ritz value). |
  593. c %--------------------------------------------------%
  594. c
  595. do 35 j = 1, nev0
  596. temp = max( eps23, abs(ritz(j)) )
  597. bounds(j) = bounds(j)/temp
  598. 35 continue
  599. c
  600. c %----------------------------------------------------%
  601. c | Sort the Ritz values according to the scaled Ritz |
  602. c | esitmates. This will push all the converged ones |
  603. c | towards the front of ritzr, ritzi, bounds |
  604. c | (in the case when NCONV < NEV.) |
  605. c %----------------------------------------------------%
  606. c
  607. wprime = 'LA'
  608. call dsortr(wprime, .true., nev0, bounds, ritz)
  609. c
  610. c %----------------------------------------------%
  611. c | Scale the Ritz estimate back to its original |
  612. c | value. |
  613. c %----------------------------------------------%
  614. c
  615. do 40 j = 1, nev0
  616. temp = max( eps23, abs(ritz(j)) )
  617. bounds(j) = bounds(j)*temp
  618. 40 continue
  619. c
  620. c %--------------------------------------------------%
  621. c | Sort the "converged" Ritz values again so that |
  622. c | the "threshold" values and their associated Ritz |
  623. c | estimates appear at the appropriate position in |
  624. c | ritz and bound. |
  625. c %--------------------------------------------------%
  626. c
  627. if (which .eq. 'BE') then
  628. c
  629. c %------------------------------------------------%
  630. c | Sort the "converged" Ritz values in increasing |
  631. c | order. The "threshold" values are in the |
  632. c | middle. |
  633. c %------------------------------------------------%
  634. c
  635. wprime = 'LA'
  636. call dsortr(wprime, .true., nconv, ritz, bounds)
  637. c
  638. else
  639. c
  640. c %----------------------------------------------%
  641. c | In LM, SM, LA, SA case, sort the "converged" |
  642. c | Ritz values according to WHICH so that the |
  643. c | "threshold" value appears at the front of |
  644. c | ritz. |
  645. c %----------------------------------------------%
  646.  
  647. call dsortr(which, .true., nconv, ritz, bounds)
  648. c
  649. end if
  650. c
  651. c %------------------------------------------%
  652. c | Use h( 1,1 ) as storage to communicate |
  653. c | rnorm to _seupd if needed |
  654. c %------------------------------------------%
  655. c
  656. h(1,1) = rnorm
  657. c
  658. c if (msglvl .gt. 1) then
  659. c call dvout ( kplusp, ritz, ndigit,
  660. c & '_saup2: Sorted Ritz values.')
  661. c call dvout ( kplusp, bounds, ndigit,
  662. c & '_saup2: Sorted ritz estimates.')
  663. c end if
  664. c c
  665. c %------------------------------------%
  666. c | Max iterations have been exceeded. |
  667. c %------------------------------------%
  668. c
  669. if (iter .gt. mxiter .and. nconv .lt. nev) info = 1
  670. c
  671. c %---------------------%
  672. c | No shifts to apply. |
  673. c %---------------------%
  674. c
  675. if (np .eq. 0 .and. nconv .lt. nev0) info = 2
  676. c
  677. np = nconv
  678. go to 1100
  679. c
  680. else if (nconv .lt. nev .and. ishift .eq. 1) then
  681. c
  682. c %---------------------------------------------------%
  683. c | Do not have all the requested eigenvalues yet. |
  684. c | To prevent possible stagnation, adjust the number |
  685. c | of Ritz values and the shifts. |
  686. c %---------------------------------------------------%
  687. c
  688. nevbef = nev
  689. nev = nev + min (nconv, np/2)
  690. if (nev .eq. 1 .and. kplusp .ge. 6) then
  691. nev = kplusp / 2
  692. else if (nev .eq. 1 .and. kplusp .gt. 2) then
  693. nev = 2
  694. end if
  695. np = kplusp - nev
  696. c
  697. c %---------------------------------------%
  698. c | If the size of NEV was just increased |
  699. c | resort the eigenvalues. |
  700. c %---------------------------------------%
  701. c
  702. if (nevbef .lt. nev)
  703. & call dsgets (ishift, which, nev, np, ritz, bounds,
  704. & workl)
  705. c
  706. end if
  707. c
  708. if (msglvl .gt. 0) then
  709. call ivout ( 1, nconv, ndigit,
  710. & '_saup2: no. of "converged" Ritz values at this iter.')
  711. c if (msglvl .gt. 1) then
  712. c kp(1) = nev
  713. c kp(2) = np
  714. c call ivout ( 2, kp, ndigit,
  715. c & '_saup2: NEV and NP are')
  716. c call dvout ( nev, ritz(np+1), ndigit,
  717. c & '_saup2: "wanted" Ritz values.')
  718. c call dvout ( nev, bounds(np+1), ndigit,
  719. c & '_saup2: Ritz estimates of the "wanted" values ')
  720. c end if
  721. end if
  722.  
  723. c
  724. if (ishift .eq. 0) then
  725. c
  726. c %-----------------------------------------------------%
  727. c | User specified shifts: reverse communication to |
  728. c | compute the shifts. They are returned in the first |
  729. c | NP locations of WORKL. |
  730. c %-----------------------------------------------------%
  731. c
  732. ushift = .true.
  733. ido = 3
  734. go to 9000
  735. end if
  736. c
  737. 50 continue
  738. c
  739. c %------------------------------------%
  740. c | Back from reverse communication; |
  741. c | User specified shifts are returned |
  742. c | in WORKL(1:*NP) |
  743. c %------------------------------------%
  744. c
  745. ushift = .false.
  746. c
  747. c
  748. c %---------------------------------------------------------%
  749. c | Move the NP shifts to the first NP locations of RITZ to |
  750. c | free up WORKL. This is for the non-exact shift case; |
  751. c | in the exact shift case, dsgets already handles this. |
  752. c %---------------------------------------------------------%
  753. c
  754. if (ishift .eq. 0) call dcopy (np, workl, 1, ritz, 1)
  755. c
  756. c if (msglvl .gt. 2) then
  757. c call ivout ( 1, np, ndigit,
  758. c & '_saup2: The number of shifts to apply ')
  759. c call dvout ( np, workl, ndigit,
  760. c & '_saup2: shifts selected')
  761. c if (ishift .eq. 1) then
  762. c call dvout ( np, bounds, ndigit,
  763. c & '_saup2: corresponding Ritz estimates')
  764. c end if
  765. c end if
  766. c
  767. c %---------------------------------------------------------%
  768. c | Apply the NP0 implicit shifts by QR bulge chasing. |
  769. c | Each shift is applied to the entire tridiagonal matrix. |
  770. c | The first 2*N locations of WORKD are used as workspace. |
  771. c | After dsapps is done, we have a Lanczos |
  772. c | factorization of length NEV. |
  773. c %---------------------------------------------------------%
  774. c
  775. call dsapps (n, nev, np, ritz, v, ldv, h, ldh, resid, q, ldq,
  776. & workd)
  777. c
  778. c %---------------------------------------------%
  779. c | Compute the B-norm of the updated residual. |
  780. c | Keep B*RESID in WORKD(1:N) to be used in |
  781. c | the first step of the next call to dsaitr. |
  782. c %---------------------------------------------%
  783. c
  784. cnorm = .true.
  785. * call arscnd (t2)
  786. if (bmat .eq. 'G') then
  787. nbx = nbx + 1
  788. call dcopy (n, resid, 1, workd(n+1), 1)
  789. ipntr(1) = n + 1
  790. ipntr(2) = 1
  791. ido = 2
  792. c
  793. c %----------------------------------%
  794. c | Exit in order to compute B*RESID |
  795. c %----------------------------------%
  796. c
  797. go to 9000
  798. else if (bmat .eq. 'I') then
  799. call dcopy (n, resid, 1, workd, 1)
  800. end if
  801. c
  802. 100 continue
  803. c
  804. c %----------------------------------%
  805. c | Back from reverse communication; |
  806. c | WORKD(1:N) := B*RESID |
  807. c %----------------------------------%
  808. c
  809. if (bmat .eq. 'G') then
  810. * call arscnd (t3)
  811. c tmvbx = tmvbx + (t3 - t2)
  812. end if
  813. c
  814. if (bmat .eq. 'G') then
  815. rnorm = ddot (n, resid, 1, workd, 1)
  816. rnorm = sqrt(abs(rnorm))
  817. else if (bmat .eq. 'I') then
  818. rnorm = dnrm2(n, resid, 1)
  819. end if
  820. cnorm = .false.
  821. 130 continue
  822. c
  823. c if (msglvl .gt. 2) then
  824. c call dvout ( 1, rnorm, ndigit,
  825. c & '_saup2: B-norm of residual for NEV factorization')
  826. c call dvout ( nev, h(1,2), ndigit,
  827. c & '_saup2: main diagonal of compressed H matrix')
  828. c call dvout ( nev-1, h(2,1), ndigit,
  829. c & '_saup2: subdiagonal of compressed H matrix')
  830. c end if
  831. c
  832. go to 1000
  833. c
  834. c %---------------------------------------------------------------%
  835. c | |
  836. c | E N D O F M A I N I T E R A T I O N L O O P |
  837. c | |
  838. c %---------------------------------------------------------------%
  839. c
  840. 1100 continue
  841. c
  842. mxiter = iter
  843. nev = nconv
  844. c
  845. 1200 continue
  846. ido = 99
  847. c
  848. c %------------%
  849. c | Error exit |
  850. c %------------%
  851. c
  852. * call arscnd (t1)
  853. c tsaup2 = t1 - t0
  854. c
  855. 9000 continue
  856. ITRAK(1)=NOPX
  857. ITRAK(2)=NBX
  858. ITRAK(3)=NRORTH
  859. ITRAK(4)=NITREF
  860. ITRAK(5)=NRSTRT
  861. c
  862. c %---------------%
  863. c | End of dsaup2 |
  864. c %---------------%
  865. c
  866. end
  867.  
  868.  
  869.  
  870.  

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