Télécharger dseupd.eso

Retour à la liste

Numérotation des lignes :

  1. C DSEUPD SOURCE BP208322 20/02/06 21:15:31 10512
  2. c\BeginDoc
  3. c
  4. c\Name: dseupd
  5. c
  6. c\Description:
  7. c
  8. c This subroutine returns the converged approximations to eigenvalues
  9. c of A*z = lambda*B*z and (optionally):
  10. c
  11. c (1) the corresponding approximate eigenvectors,
  12. c
  13. c (2) an orthonormal (Lanczos) basis for the associated approximate
  14. c invariant subspace,
  15. c
  16. c (3) Both.
  17. c
  18. c There is negligible additional cost to obtain eigenvectors. An orthonormal
  19. c (Lanczos) basis is always computed. There is an additional storage cost
  20. c of n*nev if both are requested (in this case a separate array Z must be
  21. c supplied).
  22. c
  23. c These quantities are obtained from the Lanczos factorization computed
  24. c by DSAUPD for the linear operator OP prescribed by the MODE selection
  25. c (see IPARAM(7) in DSAUPD documentation.) DSAUPD must be called before
  26. c this routine is called. These approximate eigenvalues and vectors are
  27. c commonly called Ritz values and Ritz vectors respectively. They are
  28. c referred to as such in the comments that follow. The computed orthonormal
  29. c basis for the invariant subspace corresponding to these Ritz values is
  30. c referred to as a Lanczos basis.
  31. c
  32. c See documentation in the header of the subroutine DSAUPD for a definition
  33. c of OP as well as other terms and the relation of computed Ritz values
  34. c and vectors of OP with respect to the given problem A*z = lambda*B*z.
  35. c
  36. c The approximate eigenvalues of the original problem are returned in
  37. c ascending algebraic order. The user may elect to call this routine
  38. c once for each desired Ritz vector and store it peripherally if desired.
  39. c There is also the option of computing a selected set of these vectors
  40. c with a single call.
  41. c
  42. c\Usage:
  43. c call dseupd
  44. c ( RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, BMAT, N, WHICH, NEV, TOL,
  45. c RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD, WORKL, LWORKL, INFO )
  46. c
  47. c RVEC LOGICAL (INPUT)
  48. c Specifies whether Ritz vectors corresponding to the Ritz value
  49. c approximations to the eigenproblem A*z = lambda*B*z are computed.
  50. c
  51. c RVEC = .FALSE. Compute Ritz values only.
  52. c
  53. c RVEC = .TRUE. Compute Ritz vectors.
  54. c
  55. c HOWMNY Character*1 (INPUT)
  56. c Specifies how many Ritz vectors are wanted and the form of Z
  57. c the matrix of Ritz vectors. See remark 1 below.
  58. c = 'A': compute NEV Ritz vectors;
  59. c = 'S': compute some of the Ritz vectors, specified
  60. c by the logical array SELECT.
  61. c
  62. c SELECT Logical array of dimension NCV. (INPUT/WORKSPACE)
  63. c If HOWMNY = 'S', SELECT specifies the Ritz vectors to be
  64. c computed. To select the Ritz vector corresponding to a
  65. c Ritz value D(j), SELECT(j) must be set to .TRUE..
  66. c If HOWMNY = 'A' , SELECT is used as a workspace for
  67. c reordering the Ritz values.
  68. c
  69. c D REAL*8 array of dimension NEV. (OUTPUT)
  70. c On exit, D contains the Ritz value approximations to the
  71. c eigenvalues of A*z = lambda*B*z. The values are returned
  72. c in ascending order. If IPARAM(7) = 3,4,5 then D represents
  73. c the Ritz values of OP computed by dsaupd transformed to
  74. c those of the original eigensystem A*z = lambda*B*z. If
  75. c IPARAM(7) = 1,2 then the Ritz values of OP are the same
  76. c as the those of A*z = lambda*B*z.
  77. c
  78. c Z REAL*8 N by NEV array if HOWMNY = 'A'. (OUTPUT)
  79. c On exit, Z contains the B-orthonormal Ritz vectors of the
  80. c eigensystem A*z = lambda*B*z corresponding to the Ritz
  81. c value approximations.
  82. c If RVEC = .FALSE. then Z is not referenced.
  83. c NOTE: The array Z may be set equal to first NEV columns of the
  84. c Arnoldi/Lanczos basis array V computed by DSAUPD .
  85. c
  86. c LDZ Integer. (INPUT)
  87. c The leading dimension of the array Z. If Ritz vectors are
  88. c desired, then LDZ .ge. max( 1, N ). In any case, LDZ .ge. 1.
  89. c
  90. c SIGMA Double precision (INPUT)
  91. c If IPARAM(7) = 3,4,5 represents the shift. Not referenced if
  92. c IPARAM(7) = 1 or 2.
  93. c
  94. c
  95. c **** The remaining arguments MUST be the same as for the ****
  96. c **** call to DSAUPD that was just completed. ****
  97. c
  98. c NOTE: The remaining arguments
  99. c
  100. c BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR,
  101. c WORKD, WORKL, LWORKL, INFO
  102. c
  103. c must be passed directly to DSEUPD following the last call
  104. c to DSAUPD . These arguments MUST NOT BE MODIFIED between
  105. c the the last call to DSAUPD and the call to DSEUPD .
  106. c
  107. c Two of these parameters (WORKL, INFO) are also output parameters:
  108. c
  109. c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)
  110. c WORKL(1:4*ncv) contains information obtained in
  111. c dsaupd . They are not changed by dseupd .
  112. c WORKL(4*ncv+1:ncv*ncv+8*ncv) holds the
  113. c untransformed Ritz values, the computed error estimates,
  114. c and the associated eigenvector matrix of H.
  115. c
  116. c Note: IPNTR(8:10) contains the pointer into WORKL for addresses
  117. c of the above information computed by dseupd .
  118. c -------------------------------------------------------------
  119. c IPNTR(8): pointer to the NCV RITZ values of the original system.
  120. c IPNTR(9): pointer to the NCV corresponding error bounds.
  121. c IPNTR(10): pointer to the NCV by NCV matrix of eigenvectors
  122. c of the tridiagonal matrix T. Only referenced by
  123. c dseupd if RVEC = .TRUE. See Remarks.
  124. c -------------------------------------------------------------
  125. c
  126. c INFO Integer. (OUTPUT)
  127. c Error flag on output.
  128. c = 0: Normal exit.
  129. c = -1: N must be positive.
  130. c = -2: NEV must be positive.
  131. c = -3: NCV must be greater than NEV and less than or equal to N.
  132. c = -5: WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.
  133. c = -6: BMAT must be one of 'I' or 'G'.
  134. c = -7: Length of private work WORKL array is not sufficient.
  135. c = -8: Error return from trid. eigenvalue calculation;
  136. c Information error from LAPACK routine dsteqr .
  137. c = -9: Starting vector is zero.
  138. c = -10: IPARAM(7) must be 1,2,3,4,5.
  139. c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
  140. c = -12: NEV and WHICH = 'BE' are incompatible.
  141. c = -14: DSAUPD did not find any eigenvalues to sufficient
  142. c accuracy.
  143. c = -15: HOWMNY must be one of 'A' or 'S' if RVEC = .true.
  144. c = -16: HOWMNY = 'S' not yet implemented
  145. c = -17: DSEUPD got a different count of the number of converged
  146. c Ritz values than DSAUPD got. This indicates the user
  147. c probably made an error in passing data from DSAUPD to
  148. c DSEUPD or that the data was modified before entering
  149. c DSEUPD .
  150. c
  151. c\BeginLib
  152. c
  153. c\References:
  154. c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
  155. c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
  156. c pp 357-385.
  157. c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
  158. c Restarted Arnoldi Iteration", Rice University Technical Report
  159. c TR95-13, Department of Computational and Applied Mathematics.
  160. c 3. B.N. Parlett, "The Symmetric Eigenvalue Problem". Prentice-Hall,
  161. c 1980.
  162. c 4. B.N. Parlett, B. Nour-Omid, "Towards a Black Box Lanczos Program",
  163. c Computer Physics Communications, 53 (1989), pp 169-179.
  164. c 5. B. Nour-Omid, B.N. Parlett, T. Ericson, P.S. Jensen, "How to
  165. c Implement the Spectral Transformation", Math. Comp., 48 (1987),
  166. c pp 663-673.
  167. c 6. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos
  168. c Algorithm for Solving Sparse Symmetric Generalized Eigenproblems",
  169. c SIAM J. Matr. Anal. Apps., January (1993).
  170. c 7. L. Reichel, W.B. Gragg, "Algorithm 686: FORTRAN Subroutines
  171. c for Updating the QR decomposition", ACM TOMS, December 1990,
  172. c Volume 16 Number 4, pp 369-377.
  173. c
  174. c\Remarks
  175. c 1. The converged Ritz values are always returned in increasing
  176. c (algebraic) order.
  177. c
  178. c 2. Currently only HOWMNY = 'A' is implemented. It is included at this
  179. c stage for the user who wants to incorporate it.
  180. c
  181. c\Routines called:
  182. c dsesrt ARPACK routine that sorts an array X, and applies the
  183. c corresponding permutation to a matrix A.
  184. c dsortr dsortr ARPACK sorting routine.
  185. c ivout ARPACK utility routine that prints integers.
  186. c dvout ARPACK utility routine that prints vectors.
  187. c dgeqr2 LAPACK routine that computes the QR factorization of
  188. c a matrix.
  189. c dlacpy LAPACK matrix copy routine.
  190. c dlamch LAPACK routine that determines machine constants.
  191. c dorm2r LAPACK routine that applies an orthogonal matrix in
  192. c factored form.
  193. c dsteqr LAPACK routine that computes eigenvalues and eigenvectors
  194. c of a tridiagonal matrix.
  195. c dger Level 2 BLAS rank one update to a matrix.
  196. c dcopy Level 1 BLAS that copies one vector to another .
  197. c dnrm2 Level 1 BLAS that computes the norm of a vector.
  198. c dscal Level 1 BLAS that scales a vector.
  199. c dswap Level 1 BLAS that swaps the contents of two vectors.
  200.  
  201. c\Authors
  202. c Danny Sorensen Phuong Vu
  203. c Richard Lehoucq CRPC / Rice University
  204. c Chao Yang Houston, Texas
  205. c Dept. of Computational &
  206. c Applied Mathematics
  207. c Rice University
  208. c Houston, Texas
  209. c
  210. c\Revision history:
  211. c 12/15/93: Version ' 2.1'
  212. c
  213. c\SCCS Information: @(#)
  214. c FILE: seupd.F SID: 2.11 DATE OF SID: 04/10/01 RELEASE: 2
  215. c
  216. c\EndLib
  217. c
  218. c-----------------------------------------------------------------------
  219. subroutine dseupd (rvec , howmny, select, d ,
  220. & z , ldz , sigma , bmat ,
  221. & n , which , nev , tol ,
  222. & resid , ncv , v , ldv ,
  223. & iparam, ipntr , workd , workl,
  224. & lworkl, info )
  225. c
  226. c %----------------------------------------------------%
  227. c | Include files for debugging and timing information |
  228. c -INC TARTRAK
  229. c %----------------------------------------------------%
  230. c
  231. c
  232. c %------------------%
  233. c | Scalar Arguments |
  234. c %------------------%
  235. c
  236. character bmat, howmny, which*2
  237. logical rvec
  238. integer info, ldz, ldv, lworkl, n, ncv, nev
  239. REAL*8
  240. & sigma, tol
  241. c
  242. c %-----------------%
  243. c | Array Arguments |
  244. c %-----------------%
  245. c
  246. cbp integer iparam(7), ipntr(11)
  247. integer iparam(11), ipntr(11)
  248. logical select(ncv)
  249. REAL*8
  250. & d(nev) , resid(n) , v(ldv,ncv),
  251. & z(ldz, nev), workd(2*n), workl(lworkl)
  252. c
  253. c %------------%
  254. c | Parameters |
  255. c %------------%
  256. c
  257. REAL*8
  258. & one, zero
  259. parameter (one = 1.0D+0 , zero = 0.0D+0 )
  260. c
  261. c %---------------%
  262. c | Local Scalars |
  263. c %---------------%
  264. c
  265. character type*6
  266. integer bounds , ierr , ih , ihb , ihd ,
  267. & iq , iw , j , k , ldh ,
  268. & ldq , mode , msglvl, nconv , next ,
  269. & ritz , irz , ibd , np , ishift,
  270. & leftptr, rghtptr, numcnv, jj
  271. parameter (msglvl=0)
  272. REAL*8
  273. & bnorm2 , rnorm, temp, temp1, eps23
  274. logical reord
  275. c
  276. c %----------------------%
  277. c | External Subroutines |
  278. c %----------------------%
  279. c
  280. external dcopy , dger , dgeqr2 , dlacpy , dorm2r , dscal ,
  281. c
  282. c %--------------------%
  283. c | External Functions |
  284. c %--------------------%
  285. c
  286. REAL*8
  287. external dnrm2 , dlamch
  288. c
  289. c %---------------------%
  290. **c | Intrinsic Functions |
  291. **c %---------------------%
  292. **c
  293. ** intrinsic min
  294. **c
  295. **c %-----------------------%
  296. **c | Executable Statements |
  297. c %-----------------------%
  298. c
  299. c %------------------------%
  300. c | Set default parameters |
  301. c %------------------------%
  302. c
  303. c msglvl = mseupd
  304. mode = iparam(7)
  305. nconv = iparam(5)
  306. info = 0
  307. c
  308. c %--------------%
  309. c | Quick return |
  310. c %--------------%
  311. c
  312. if (nconv .eq. 0) go to 9000
  313. ierr = 0
  314. c
  315. if (nconv .le. 0) ierr = -14
  316. if (n .le. 0) ierr = -1
  317. if (nev .le. 0) ierr = -2
  318. if (ncv .le. nev .or. ncv .gt. n) ierr = -3
  319. if (which .ne. 'LM' .and.
  320. & which .ne. 'SM' .and.
  321. & which .ne. 'LA' .and.
  322. & which .ne. 'SA' .and.
  323. & which .ne. 'BE') ierr = -5
  324. if (bmat .ne. 'I' .and. bmat .ne. 'G') ierr = -6
  325. if ( (howmny .ne. 'A' .and.
  326. & howmny .ne. 'P' .and.
  327. & howmny .ne. 'S') .and. rvec )
  328. & ierr = -15
  329. if (rvec .and. howmny .eq. 'S') ierr = -16
  330. c
  331. if (rvec .and. lworkl .lt. ncv**2+8*ncv) ierr = -7
  332. c
  333. if (mode .eq. 1 .or. mode .eq. 2) then
  334. type = 'REGULR'
  335. else if (mode .eq. 3 ) then
  336. type = 'SHIFTI'
  337. else if (mode .eq. 4 ) then
  338. type = 'BUCKLE'
  339. else if (mode .eq. 5 ) then
  340. type = 'CAYLEY'
  341. else
  342. ierr = -10
  343. end if
  344. if (mode .eq. 1 .and. bmat .eq. 'G') ierr = -11
  345. if (nev .eq. 1 .and. which .eq. 'BE') ierr = -12
  346. c
  347. c %------------%
  348. c | Error Exit |
  349. c %------------%
  350. c
  351. if (ierr .ne. 0) then
  352. info = ierr
  353. go to 9000
  354. end if
  355. c
  356. c %-------------------------------------------------------%
  357. c | Pointer into WORKL for address of H, RITZ, BOUNDS, Q |
  358. c | etc... and the remaining workspace. |
  359. c | Also update pointer to be used on output. |
  360. c | Memory is laid out as follows: |
  361. c | workl(1:2*ncv) := generated tridiagonal matrix H |
  362. c | The subdiagonal is stored in workl(2:ncv). |
  363. c | The dead spot is workl(1) but upon exiting |
  364. c | dsaupd stores the B-norm of the last residual |
  365. c | vector in workl(1). We use this !!! |
  366. c | workl(2*ncv+1:2*ncv+ncv) := ritz values |
  367. c | The wanted values are in the first NCONV spots. |
  368. c | workl(3*ncv+1:3*ncv+ncv) := computed Ritz estimates |
  369. c | The wanted values are in the first NCONV spots. |
  370. c | NOTE: workl(1:4*ncv) is set by dsaupd and is not |
  371. c | modified by dseupd . |
  372. c %-------------------------------------------------------%
  373. c
  374. c %-------------------------------------------------------%
  375. c | The following is used and set by dseupd . |
  376. c | workl(4*ncv+1:4*ncv+ncv) := used as workspace during |
  377. c | computation of the eigenvectors of H. Stores |
  378. c | the diagonal of H. Upon EXIT contains the NCV |
  379. c | Ritz values of the original system. The first |
  380. c | NCONV spots have the wanted values. If MODE = |
  381. c | 1 or 2 then will equal workl(2*ncv+1:3*ncv). |
  382. c | workl(5*ncv+1:5*ncv+ncv) := used as workspace during |
  383. c | computation of the eigenvectors of H. Stores |
  384. c | the subdiagonal of H. Upon EXIT contains the |
  385. c | NCV corresponding Ritz estimates of the |
  386. c | original system. The first NCONV spots have the |
  387. c | wanted values. If MODE = 1,2 then will equal |
  388. c | workl(3*ncv+1:4*ncv). |
  389. c | workl(6*ncv+1:6*ncv+ncv*ncv) := orthogonal Q that is |
  390. c | the eigenvector matrix for H as returned by |
  391. c | dsteqr . Not referenced if RVEC = .False. |
  392. c | Ordering follows that of workl(4*ncv+1:5*ncv) |
  393. c | workl(6*ncv+ncv*ncv+1:6*ncv+ncv*ncv+2*ncv) := |
  394. c | Workspace. Needed by dsteqr and by dseupd . |
  395. c | GRAND total of NCV*(NCV+8) locations. |
  396. c %-------------------------------------------------------%
  397. c
  398. c
  399. ih = ipntr(5)
  400. ritz = ipntr(6)
  401. bounds = ipntr(7)
  402. ldh = ncv
  403. ldq = ncv
  404. ihd = bounds + ldh
  405. ihb = ihd + ldh
  406. iq = ihb + ldh
  407. iw = iq + ldh*ncv
  408. next = iw + 2*ncv
  409. ipntr(4) = next
  410. ipntr(8) = ihd
  411. ipntr(9) = ihb
  412. ipntr(10) = iq
  413. c
  414. c %----------------------------------------%
  415. c | irz points to the Ritz values computed |
  416. c | by _seigt before exiting _saup2. |
  417. c | ibd points to the Ritz estimates |
  418. c | computed by _seigt before exiting |
  419. c | _saup2. |
  420. c %----------------------------------------%
  421. c
  422. irz = ipntr(11)+ncv
  423. ibd = irz+ncv
  424. c
  425. c
  426. c %---------------------------------%
  427. c | Set machine dependent constant. |
  428. c %---------------------------------%
  429. c
  430. eps23 = dlamch ('Epsilon-Machine')
  431. eps23 = eps23**(2.0D+0 / 3.0D+0 )
  432. c
  433. c %---------------------------------------%
  434. c | RNORM is B-norm of the RESID(1:N). |
  435. c | BNORM2 is the 2 norm of B*RESID(1:N). |
  436. c | Upon exit of dsaupd WORKD(1:N) has |
  437. c | B*RESID(1:N). |
  438. c %---------------------------------------%
  439. c
  440. rnorm = workl(ih)
  441. if (bmat .eq. 'I') then
  442. bnorm2 = rnorm
  443. else if (bmat .eq. 'G') then
  444. bnorm2 = dnrm2 (n, workd, 1)
  445. end if
  446. c
  447. c if (msglvl .gt. 2) then
  448. c call dvout ( ncv, workl(irz), ndigit,
  449. c & '_seupd: Ritz values passed in from _SAUPD.')
  450. c call dvout ( ncv, workl(ibd), ndigit,
  451. c & '_seupd: Ritz estimates passed in from _SAUPD.')
  452. c end if
  453. c
  454. if (rvec) then
  455. c
  456. reord = .false.
  457. c
  458. c %---------------------------------------------------%
  459. c | Use the temporary bounds array to store indices |
  460. c | These will be used to mark the select array later |
  461. c %---------------------------------------------------%
  462. c
  463. do 10 j = 1,ncv
  464. workl(bounds+j-1) = j
  465. select(j) = .false.
  466. 10 continue
  467. c
  468. c %-------------------------------------%
  469. c | Select the wanted Ritz values. |
  470. c | Sort the Ritz values so that the |
  471. c | wanted ones appear at the tailing |
  472. c | NEV positions of workl(irr) and |
  473. c | workl(iri). Move the corresponding |
  474. c | error estimates in workl(bound) |
  475. c | accordingly. |
  476. c %-------------------------------------%
  477. c
  478. np = ncv - nev
  479. ishift = 0
  480. call dsgets (ishift, which , nev ,
  481. & np , workl(irz) , workl(bounds),
  482. & workl)
  483. c
  484. c if (msglvl .gt. 2) then
  485. c call dvout ( ncv, workl(irz), ndigit,
  486. c & '_seupd: Ritz values after calling _SGETS.')
  487. c call dvout ( ncv, workl(bounds), ndigit,
  488. c & '_seupd: Ritz value indices after calling _SGETS.')
  489. c end if
  490. c
  491. c %-----------------------------------------------------%
  492. c | Record indices of the converged wanted Ritz values |
  493. c | Mark the select array for possible reordering |
  494. c %-----------------------------------------------------%
  495. c
  496. numcnv = 0
  497. do 11 j = 1,ncv
  498. temp1 = max(eps23, abs(workl(irz+ncv-j)) )
  499. jj = int(workl(bounds + ncv - j))
  500. if (numcnv .lt. nconv .and.
  501. & workl(ibd+jj-1) .le. tol*temp1) then
  502. select(jj) = .true.
  503. numcnv = numcnv + 1
  504. if (jj .gt. nconv) reord = .true.
  505. endif
  506. 11 continue
  507. c
  508. c %-----------------------------------------------------------%
  509. c | Check the count (numcnv) of converged Ritz values with |
  510. c | the number (nconv) reported by _saupd. If these two |
  511. c | are different then there has probably been an error |
  512. c | caused by incorrect passing of the _saupd data. |
  513. c %-----------------------------------------------------------%
  514. c
  515. c if (msglvl .gt. 2) then
  516. c call ivout ( 1, numcnv, ndigit,
  517. c & '_seupd: Number of specified eigenvalues')
  518. c call ivout ( 1, nconv, ndigit,
  519. c & '_seupd: Number of "converged" eigenvalues')
  520. c end if
  521. c
  522. if (numcnv .ne. nconv) then
  523. info = -17
  524. go to 9000
  525. end if
  526. c
  527. c %-----------------------------------------------------------%
  528. c | Call LAPACK routine _steqr to compute the eigenvalues and |
  529. c | eigenvectors of the final symmetric tridiagonal matrix H. |
  530. c | Initialize the eigenvector matrix Q to the identity. |
  531. c %-----------------------------------------------------------%
  532. c
  533. call dcopy (ncv-1, workl(ih+1), 1, workl(ihb), 1)
  534. call dcopy (ncv, workl(ih+ldh), 1, workl(ihd), 1)
  535. c
  536. call dsteqr ('Identity', ncv, workl(ihd), workl(ihb),
  537. & workl(iq) , ldq, workl(iw), ierr)
  538. c
  539. if (ierr .ne. 0) then
  540. info = -8
  541. go to 9000
  542. end if
  543. c
  544. c if (msglvl .gt. 1) then
  545. c call dcopy (ncv, workl(iq+ncv-1), ldq, workl(iw), 1)
  546. c call dvout ( ncv, workl(ihd), ndigit,
  547. c & '_seupd: NCV Ritz values of the final H matrix')
  548. c call dvout ( ncv, workl(iw), ndigit,
  549. c & '_seupd: last row of the eigenvector matrix for H')
  550. c end if
  551. c
  552. if (reord) then
  553. c
  554. c %---------------------------------------------%
  555. c | Reordered the eigenvalues and eigenvectors |
  556. c | computed by _steqr so that the "converged" |
  557. c | eigenvalues appear in the first NCONV |
  558. c | positions of workl(ihd), and the associated |
  559. c | eigenvectors appear in the first NCONV |
  560. c | columns. |
  561. c %---------------------------------------------%
  562. c
  563. leftptr = 1
  564. rghtptr = ncv
  565. c
  566. if (ncv .eq. 1) go to 30
  567. c
  568. 20 if (select(leftptr)) then
  569. c
  570. c %-------------------------------------------%
  571. c | Search, from the left, for the first Ritz |
  572. c | value that has not converged. |
  573. c %-------------------------------------------%
  574. c
  575. leftptr = leftptr + 1
  576. c
  577. else if ( .not. select(rghtptr)) then
  578. c
  579. c %----------------------------------------------%
  580. c | Search, from the right, the first Ritz value |
  581. c | that has converged. |
  582. c %----------------------------------------------%
  583. c
  584. rghtptr = rghtptr - 1
  585. c
  586. else
  587. c
  588. c %----------------------------------------------%
  589. c | Swap the Ritz value on the left that has not |
  590. c | converged with the Ritz value on the right |
  591. c | that has converged. Swap the associated |
  592. c | eigenvector of the tridiagonal matrix H as |
  593. c | well. |
  594. c %----------------------------------------------%
  595. c
  596. temp = workl(ihd+leftptr-1)
  597. workl(ihd+leftptr-1) = workl(ihd+rghtptr-1)
  598. workl(ihd+rghtptr-1) = temp
  599. call dcopy (ncv, workl(iq+ncv*(leftptr-1)), 1,
  600. & workl(iw), 1)
  601. call dcopy (ncv, workl(iq+ncv*(rghtptr-1)), 1,
  602. & workl(iq+ncv*(leftptr-1)), 1)
  603. call dcopy (ncv, workl(iw), 1,
  604. & workl(iq+ncv*(rghtptr-1)), 1)
  605. leftptr = leftptr + 1
  606. rghtptr = rghtptr - 1
  607. c
  608. end if
  609. c
  610. if (leftptr .lt. rghtptr) go to 20
  611. c
  612. end if
  613. c
  614. 30 continue
  615. c if (msglvl .gt. 2) then
  616. c call dvout ( ncv, workl(ihd), ndigit,
  617. c & '_seupd: The eigenvalues of H--reordered')
  618. c end if
  619. c
  620. c %----------------------------------------%
  621. c | Load the converged Ritz values into D. |
  622. c %----------------------------------------%
  623. c
  624. call dcopy (nconv, workl(ihd), 1, d, 1)
  625. c
  626. else
  627. c
  628. c %-----------------------------------------------------%
  629. c | Ritz vectors not required. Load Ritz values into D. |
  630. c %-----------------------------------------------------%
  631. c
  632. call dcopy (nconv, workl(ritz), 1, d, 1)
  633. call dcopy (ncv, workl(ritz), 1, workl(ihd), 1)
  634. c
  635. end if
  636. c
  637. c %------------------------------------------------------------------%
  638. c | Transform the Ritz values and possibly vectors and corresponding |
  639. c | Ritz estimates of OP to those of A*x=lambda*B*x. The Ritz values |
  640. c | (and corresponding data) are returned in ascending order. |
  641. c %------------------------------------------------------------------%
  642. c
  643. if (type .eq. 'REGULR') then
  644. c
  645. c %---------------------------------------------------------%
  646. c | Ascending sort of wanted Ritz values, vectors and error |
  647. c | bounds. Not necessary if only Ritz values are desired. |
  648. c %---------------------------------------------------------%
  649. c
  650. if (rvec) then
  651. call dsesrt ('LA', rvec , nconv, d, ncv, workl(iq), ldq)
  652. else
  653. call dcopy (ncv, workl(bounds), 1, workl(ihb), 1)
  654. end if
  655. c
  656. else
  657. c
  658. c %-------------------------------------------------------------%
  659. c | * Make a copy of all the Ritz values. |
  660. c | * Transform the Ritz values back to the original system. |
  661. c | For TYPE = 'SHIFTI' the transformation is |
  662. c | lambda = 1/theta + sigma |
  663. c | For TYPE = 'BUCKLE' the transformation is |
  664. c | lambda = sigma * theta / ( theta - 1 ) |
  665. c | For TYPE = 'CAYLEY' the transformation is |
  666. c | lambda = sigma * (theta + 1) / (theta - 1 ) |
  667. c | where the theta are the Ritz values returned by dsaupd . |
  668. c | NOTES: |
  669. c | *The Ritz vectors are not affected by the transformation. |
  670. c | They are only reordered. |
  671. c %-------------------------------------------------------------%
  672. c
  673. call dcopy (ncv, workl(ihd), 1, workl(iw), 1)
  674. if (type .eq. 'SHIFTI') then
  675. do 40 k=1, ncv
  676. workl(ihd+k-1) = one / workl(ihd+k-1) + sigma
  677. 40 continue
  678. else if (type .eq. 'BUCKLE') then
  679. do 50 k=1, ncv
  680. workl(ihd+k-1) = sigma * workl(ihd+k-1) /
  681. & (workl(ihd+k-1) - one)
  682. 50 continue
  683. else if (type .eq. 'CAYLEY') then
  684. do 60 k=1, ncv
  685. workl(ihd+k-1) = sigma * (workl(ihd+k-1) + one) /
  686. & (workl(ihd+k-1) - one)
  687. 60 continue
  688. end if
  689. c
  690. c %-------------------------------------------------------------%
  691. c | * Store the wanted NCONV lambda values into D. |
  692. c | * Sort the NCONV wanted lambda in WORKL(IHD:IHD+NCONV-1) |
  693. c | into ascending order and apply sort to the NCONV theta |
  694. c | values in the transformed system. We will need this to |
  695. c | compute Ritz estimates in the original system. |
  696. c | * Finally sort the lambda`s into ascending order and apply |
  697. c | to Ritz vectors if wanted. Else just sort lambda`s into |
  698. c | ascending order. |
  699. c | NOTES: |
  700. c | *workl(iw:iw+ncv-1) contain the theta ordered so that they |
  701. c | match the ordering of the lambda. We`ll use them again for |
  702. c | Ritz vector purification. |
  703. c %-------------------------------------------------------------%
  704. c
  705. call dcopy (nconv, workl(ihd), 1, d, 1)
  706. call dsortr ('LA', .true., nconv, workl(ihd), workl(iw))
  707. if (rvec) then
  708. call dsesrt ('LA', rvec , nconv, d, ncv, workl(iq), ldq)
  709. else
  710. call dcopy (ncv, workl(bounds), 1, workl(ihb), 1)
  711. call dscal (ncv, bnorm2/rnorm, workl(ihb), 1)
  712. call dsortr ('LA', .true., nconv, d, workl(ihb))
  713. end if
  714. c
  715. end if
  716. c
  717. c %------------------------------------------------%
  718. c | Compute the Ritz vectors. Transform the wanted |
  719. c | eigenvectors of the symmetric tridiagonal H by |
  720. c | the Lanczos basis matrix V. |
  721. c %------------------------------------------------%
  722. c
  723. if (rvec .and. howmny .eq. 'A') then
  724. c
  725. c %----------------------------------------------------------%
  726. c | Compute the QR factorization of the matrix representing |
  727. c | the wanted invariant subspace located in the first NCONV |
  728. c | columns of workl(iq,ldq). |
  729. c %----------------------------------------------------------%
  730. c
  731. call dgeqr2 (ncv, nconv , workl(iq) ,
  732. & ldq, workl(iw+ncv), workl(ihb),
  733. & ierr)
  734. c
  735. c %--------------------------------------------------------%
  736. c | * Postmultiply V by Q. |
  737. c | * Copy the first NCONV columns of VQ into Z. |
  738. c | The N by NCONV matrix Z is now a matrix representation |
  739. c | of the approximate invariant subspace associated with |
  740. c | the Ritz values in workl(ihd). |
  741. c %--------------------------------------------------------%
  742. c
  743. call dorm2r ('Right', 'Notranspose', n ,
  744. & ncv , nconv , workl(iq),
  745. & ldq , workl(iw+ncv), v ,
  746. & ldv , workd(n+1) , ierr)
  747. call dlacpy ('All', n, nconv, v, ldv, z, ldz)
  748. c
  749. c %-----------------------------------------------------%
  750. c | In order to compute the Ritz estimates for the Ritz |
  751. c | values in both systems, need the last row of the |
  752. c | eigenvector matrix. Remember, it`s in factored form |
  753. c %-----------------------------------------------------%
  754. c
  755. do 65 j = 1, ncv-1
  756. workl(ihb+j-1) = zero
  757. 65 continue
  758. workl(ihb+ncv-1) = one
  759. call dorm2r ('Left', 'Transpose' , ncv ,
  760. & 1 , nconv , workl(iq) ,
  761. & ldq , workl(iw+ncv), workl(ihb),
  762. & ncv , temp , ierr)
  763. c
  764. c %-----------------------------------------------------%
  765. c | Make a copy of the last row into |
  766. c | workl(iw+ncv:iw+2*ncv), as it is needed again in |
  767. c | the Ritz vector purification step below |
  768. c %-----------------------------------------------------%
  769. c
  770. do 67 j = 1, nconv
  771. workl(iw+ncv+j-1) = workl(ihb+j-1)
  772. 67 continue
  773.  
  774. else if (rvec .and. howmny .eq. 'S') then
  775. c
  776. c Not yet implemented. See remark 2 above.
  777. c
  778. end if
  779. c
  780. if (type .eq. 'REGULR' .and. rvec) then
  781. c
  782. do 70 j=1, ncv
  783. workl(ihb+j-1) = rnorm * abs( workl(ihb+j-1) )
  784. 70 continue
  785. c
  786. else if (type .ne. 'REGULR' .and. rvec) then
  787. c
  788. c %-------------------------------------------------%
  789. c | * Determine Ritz estimates of the theta. |
  790. c | If RVEC = .true. then compute Ritz estimates |
  791. c | of the theta. |
  792. c | If RVEC = .false. then copy Ritz estimates |
  793. c | as computed by dsaupd . |
  794. c | * Determine Ritz estimates of the lambda. |
  795. c %-------------------------------------------------%
  796. c
  797. call dscal (ncv, bnorm2, workl(ihb), 1)
  798. if (type .eq. 'SHIFTI') then
  799. c
  800. do 80 k=1, ncv
  801. workl(ihb+k-1) = abs( workl(ihb+k-1) )
  802. & / workl(iw+k-1)**2
  803. 80 continue
  804. c
  805. else if (type .eq. 'BUCKLE') then
  806. c
  807. do 90 k=1, ncv
  808. workl(ihb+k-1) = sigma * abs( workl(ihb+k-1) )
  809. & / (workl(iw+k-1)-one )**2
  810. 90 continue
  811. c
  812. else if (type .eq. 'CAYLEY') then
  813. c
  814. do 100 k=1, ncv
  815. workl(ihb+k-1) = abs( workl(ihb+k-1)
  816. & / workl(iw+k-1)*(workl(iw+k-1)-one) )
  817. 100 continue
  818. c
  819. end if
  820. c
  821. end if
  822. c
  823. c if (type .ne. 'REGULR' .and. msglvl .gt. 1) then
  824. c call dvout ( nconv, d, ndigit,
  825. c & '_seupd: Untransformed converged Ritz values')
  826. c call dvout ( nconv, workl(ihb), ndigit,
  827. c & '_seupd: Ritz estimates of the untransformed Ritz values')
  828. c else if (msglvl .gt. 1) then
  829. c call dvout ( nconv, d, ndigit,
  830. c & '_seupd: Converged Ritz values')
  831. c call dvout ( nconv, workl(ihb), ndigit,
  832. c & '_seupd: Associated Ritz estimates')
  833. c end if
  834. c
  835. c %-------------------------------------------------%
  836. c | Ritz vector purification step. Formally perform |
  837. c | one of inverse subspace iteration. Only used |
  838. c | for MODE = 3,4,5. See reference 7 |
  839. c %-------------------------------------------------%
  840. c
  841. if (rvec .and. (type .eq. 'SHIFTI' .or. type .eq. 'CAYLEY')) then
  842. c
  843. do 110 k=0, nconv-1
  844. workl(iw+k) = workl(iw+ncv+k)
  845. & / workl(iw+k)
  846. 110 continue
  847. c
  848. else if (rvec .and. type .eq. 'BUCKLE') then
  849. c
  850. do 120 k=0, nconv-1
  851. workl(iw+k) = workl(iw+ncv+k)
  852. & / (workl(iw+k)-one)
  853. 120 continue
  854. c
  855. end if
  856. c
  857. if (rvec .and. type .ne. 'REGULR')
  858. & call dger (n, nconv, one, resid, 1, workl(iw), 1, z, ldz)
  859. c
  860. 9000 continue
  861. c
  862. c return
  863. c
  864. c %---------------%
  865. c | End of dseupd |
  866. c %---------------%
  867. c
  868. end
  869.  
  870.  
  871.  

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