Télécharger dseupd.eso

Retour à la liste

Numérotation des lignes :

  1. C DSEUPD SOURCE BP208322 15/10/13 21:15:53 8670
  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. -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. integer iparam(7), ipntr(11)
  247. logical select(ncv)
  248. REAL*8
  249. & d(nev) , resid(n) , v(ldv,ncv),
  250. & z(ldz, nev), workd(2*n), workl(lworkl)
  251. c
  252. c %------------%
  253. c | Parameters |
  254. c %------------%
  255. c
  256. REAL*8
  257. & one, zero
  258. parameter (one = 1.0D+0 , zero = 0.0D+0 )
  259. c
  260. c %---------------%
  261. c | Local Scalars |
  262. c %---------------%
  263. c
  264. character type*6
  265. integer bounds , ierr , ih , ihb , ihd ,
  266. & iq , iw , j , k , ldh ,
  267. & ldq , mode , msglvl, nconv , next ,
  268. & ritz , irz , ibd , np , ishift,
  269. & leftptr, rghtptr, numcnv, jj
  270. REAL*8
  271. & bnorm2 , rnorm, temp, temp1, eps23
  272. logical reord
  273. c
  274. c %----------------------%
  275. c | External Subroutines |
  276. c %----------------------%
  277. c
  278. external dcopy , dger , dgeqr2 , dlacpy , dorm2r , dscal ,
  279. c
  280. c %--------------------%
  281. c | External Functions |
  282. c %--------------------%
  283. c
  284. REAL*8
  285. external dnrm2 , dlamch
  286. c
  287. c %---------------------%
  288. **c | Intrinsic Functions |
  289. **c %---------------------%
  290. **c
  291. ** intrinsic min
  292. **c
  293. **c %-----------------------%
  294. **c | Executable Statements |
  295. c %-----------------------%
  296. c
  297. c %------------------------%
  298. c | Set default parameters |
  299. c %------------------------%
  300. c
  301. msglvl = mseupd
  302. mode = iparam(7)
  303. nconv = iparam(5)
  304. info = 0
  305. c
  306. c %--------------%
  307. c | Quick return |
  308. c %--------------%
  309. c
  310. if (nconv .eq. 0) go to 9000
  311. ierr = 0
  312. c
  313. if (nconv .le. 0) ierr = -14
  314. if (n .le. 0) ierr = -1
  315. if (nev .le. 0) ierr = -2
  316. if (ncv .le. nev .or. ncv .gt. n) ierr = -3
  317. if (which .ne. 'LM' .and.
  318. & which .ne. 'SM' .and.
  319. & which .ne. 'LA' .and.
  320. & which .ne. 'SA' .and.
  321. & which .ne. 'BE') ierr = -5
  322. if (bmat .ne. 'I' .and. bmat .ne. 'G') ierr = -6
  323. if ( (howmny .ne. 'A' .and.
  324. & howmny .ne. 'P' .and.
  325. & howmny .ne. 'S') .and. rvec )
  326. & ierr = -15
  327. if (rvec .and. howmny .eq. 'S') ierr = -16
  328. c
  329. if (rvec .and. lworkl .lt. ncv**2+8*ncv) ierr = -7
  330. c
  331. if (mode .eq. 1 .or. mode .eq. 2) then
  332. type = 'REGULR'
  333. else if (mode .eq. 3 ) then
  334. type = 'SHIFTI'
  335. else if (mode .eq. 4 ) then
  336. type = 'BUCKLE'
  337. else if (mode .eq. 5 ) then
  338. type = 'CAYLEY'
  339. else
  340. ierr = -10
  341. end if
  342. if (mode .eq. 1 .and. bmat .eq. 'G') ierr = -11
  343. if (nev .eq. 1 .and. which .eq. 'BE') ierr = -12
  344. c
  345. c %------------%
  346. c | Error Exit |
  347. c %------------%
  348. c
  349. if (ierr .ne. 0) then
  350. info = ierr
  351. go to 9000
  352. end if
  353. c
  354. c %-------------------------------------------------------%
  355. c | Pointer into WORKL for address of H, RITZ, BOUNDS, Q |
  356. c | etc... and the remaining workspace. |
  357. c | Also update pointer to be used on output. |
  358. c | Memory is laid out as follows: |
  359. c | workl(1:2*ncv) := generated tridiagonal matrix H |
  360. c | The subdiagonal is stored in workl(2:ncv). |
  361. c | The dead spot is workl(1) but upon exiting |
  362. c | dsaupd stores the B-norm of the last residual |
  363. c | vector in workl(1). We use this !!! |
  364. c | workl(2*ncv+1:2*ncv+ncv) := ritz values |
  365. c | The wanted values are in the first NCONV spots. |
  366. c | workl(3*ncv+1:3*ncv+ncv) := computed Ritz estimates |
  367. c | The wanted values are in the first NCONV spots. |
  368. c | NOTE: workl(1:4*ncv) is set by dsaupd and is not |
  369. c | modified by dseupd . |
  370. c %-------------------------------------------------------%
  371. c
  372. c %-------------------------------------------------------%
  373. c | The following is used and set by dseupd . |
  374. c | workl(4*ncv+1:4*ncv+ncv) := used as workspace during |
  375. c | computation of the eigenvectors of H. Stores |
  376. c | the diagonal of H. Upon EXIT contains the NCV |
  377. c | Ritz values of the original system. The first |
  378. c | NCONV spots have the wanted values. If MODE = |
  379. c | 1 or 2 then will equal workl(2*ncv+1:3*ncv). |
  380. c | workl(5*ncv+1:5*ncv+ncv) := used as workspace during |
  381. c | computation of the eigenvectors of H. Stores |
  382. c | the subdiagonal of H. Upon EXIT contains the |
  383. c | NCV corresponding Ritz estimates of the |
  384. c | original system. The first NCONV spots have the |
  385. c | wanted values. If MODE = 1,2 then will equal |
  386. c | workl(3*ncv+1:4*ncv). |
  387. c | workl(6*ncv+1:6*ncv+ncv*ncv) := orthogonal Q that is |
  388. c | the eigenvector matrix for H as returned by |
  389. c | dsteqr . Not referenced if RVEC = .False. |
  390. c | Ordering follows that of workl(4*ncv+1:5*ncv) |
  391. c | workl(6*ncv+ncv*ncv+1:6*ncv+ncv*ncv+2*ncv) := |
  392. c | Workspace. Needed by dsteqr and by dseupd . |
  393. c | GRAND total of NCV*(NCV+8) locations. |
  394. c %-------------------------------------------------------%
  395. c
  396. c
  397. ih = ipntr(5)
  398. ritz = ipntr(6)
  399. bounds = ipntr(7)
  400. ldh = ncv
  401. ldq = ncv
  402. ihd = bounds + ldh
  403. ihb = ihd + ldh
  404. iq = ihb + ldh
  405. iw = iq + ldh*ncv
  406. next = iw + 2*ncv
  407. ipntr(4) = next
  408. ipntr(8) = ihd
  409. ipntr(9) = ihb
  410. ipntr(10) = iq
  411. c
  412. c %----------------------------------------%
  413. c | irz points to the Ritz values computed |
  414. c | by _seigt before exiting _saup2. |
  415. c | ibd points to the Ritz estimates |
  416. c | computed by _seigt before exiting |
  417. c | _saup2. |
  418. c %----------------------------------------%
  419. c
  420. irz = ipntr(11)+ncv
  421. ibd = irz+ncv
  422. c
  423. c
  424. c %---------------------------------%
  425. c | Set machine dependent constant. |
  426. c %---------------------------------%
  427. c
  428. eps23 = dlamch ('Epsilon-Machine')
  429. eps23 = eps23**(2.0D+0 / 3.0D+0 )
  430. c
  431. c %---------------------------------------%
  432. c | RNORM is B-norm of the RESID(1:N). |
  433. c | BNORM2 is the 2 norm of B*RESID(1:N). |
  434. c | Upon exit of dsaupd WORKD(1:N) has |
  435. c | B*RESID(1:N). |
  436. c %---------------------------------------%
  437. c
  438. rnorm = workl(ih)
  439. if (bmat .eq. 'I') then
  440. bnorm2 = rnorm
  441. else if (bmat .eq. 'G') then
  442. bnorm2 = dnrm2 (n, workd, 1)
  443. end if
  444. c
  445. if (msglvl .gt. 2) then
  446. call dvout (logfil, ncv, workl(irz), ndigit,
  447. & '_seupd: Ritz values passed in from _SAUPD.')
  448. call dvout (logfil, ncv, workl(ibd), ndigit,
  449. & '_seupd: Ritz estimates passed in from _SAUPD.')
  450. end if
  451. c
  452. if (rvec) then
  453. c
  454. reord = .false.
  455. c
  456. c %---------------------------------------------------%
  457. c | Use the temporary bounds array to store indices |
  458. c | These will be used to mark the select array later |
  459. c %---------------------------------------------------%
  460. c
  461. do 10 j = 1,ncv
  462. workl(bounds+j-1) = j
  463. select(j) = .false.
  464. 10 continue
  465. c
  466. c %-------------------------------------%
  467. c | Select the wanted Ritz values. |
  468. c | Sort the Ritz values so that the |
  469. c | wanted ones appear at the tailing |
  470. c | NEV positions of workl(irr) and |
  471. c | workl(iri). Move the corresponding |
  472. c | error estimates in workl(bound) |
  473. c | accordingly. |
  474. c %-------------------------------------%
  475. c
  476. np = ncv - nev
  477. ishift = 0
  478. call dsgets (ishift, which , nev ,
  479. & np , workl(irz) , workl(bounds),
  480. & workl)
  481. c
  482. if (msglvl .gt. 2) then
  483. call dvout (logfil, ncv, workl(irz), ndigit,
  484. & '_seupd: Ritz values after calling _SGETS.')
  485. call dvout (logfil, ncv, workl(bounds), ndigit,
  486. & '_seupd: Ritz value indices after calling _SGETS.')
  487. end if
  488. c
  489. c %-----------------------------------------------------%
  490. c | Record indices of the converged wanted Ritz values |
  491. c | Mark the select array for possible reordering |
  492. c %-----------------------------------------------------%
  493. c
  494. numcnv = 0
  495. do 11 j = 1,ncv
  496. temp1 = max(eps23, abs(workl(irz+ncv-j)) )
  497. jj = int(workl(bounds + ncv - j))
  498. if (numcnv .lt. nconv .and.
  499. & workl(ibd+jj-1) .le. tol*temp1) then
  500. select(jj) = .true.
  501. numcnv = numcnv + 1
  502. if (jj .gt. nconv) reord = .true.
  503. endif
  504. 11 continue
  505. c
  506. c %-----------------------------------------------------------%
  507. c | Check the count (numcnv) of converged Ritz values with |
  508. c | the number (nconv) reported by _saupd. If these two |
  509. c | are different then there has probably been an error |
  510. c | caused by incorrect passing of the _saupd data. |
  511. c %-----------------------------------------------------------%
  512. c
  513. if (msglvl .gt. 2) then
  514. call ivout(logfil, 1, numcnv, ndigit,
  515. & '_seupd: Number of specified eigenvalues')
  516. call ivout(logfil, 1, nconv, ndigit,
  517. & '_seupd: Number of "converged" eigenvalues')
  518. end if
  519. c
  520. if (numcnv .ne. nconv) then
  521. info = -17
  522. go to 9000
  523. end if
  524. c
  525. c %-----------------------------------------------------------%
  526. c | Call LAPACK routine _steqr to compute the eigenvalues and |
  527. c | eigenvectors of the final symmetric tridiagonal matrix H. |
  528. c | Initialize the eigenvector matrix Q to the identity. |
  529. c %-----------------------------------------------------------%
  530. c
  531. call dcopy (ncv-1, workl(ih+1), 1, workl(ihb), 1)
  532. call dcopy (ncv, workl(ih+ldh), 1, workl(ihd), 1)
  533. c
  534. call dsteqr ('Identity', ncv, workl(ihd), workl(ihb),
  535. & workl(iq) , ldq, workl(iw), ierr)
  536. c
  537. if (ierr .ne. 0) then
  538. info = -8
  539. go to 9000
  540. end if
  541. c
  542. if (msglvl .gt. 1) then
  543. call dcopy (ncv, workl(iq+ncv-1), ldq, workl(iw), 1)
  544. call dvout (logfil, ncv, workl(ihd), ndigit,
  545. & '_seupd: NCV Ritz values of the final H matrix')
  546. call dvout (logfil, ncv, workl(iw), ndigit,
  547. & '_seupd: last row of the eigenvector matrix for H')
  548. end if
  549. c
  550. if (reord) then
  551. c
  552. c %---------------------------------------------%
  553. c | Reordered the eigenvalues and eigenvectors |
  554. c | computed by _steqr so that the "converged" |
  555. c | eigenvalues appear in the first NCONV |
  556. c | positions of workl(ihd), and the associated |
  557. c | eigenvectors appear in the first NCONV |
  558. c | columns. |
  559. c %---------------------------------------------%
  560. c
  561. leftptr = 1
  562. rghtptr = ncv
  563. c
  564. if (ncv .eq. 1) go to 30
  565. c
  566. 20 if (select(leftptr)) then
  567. c
  568. c %-------------------------------------------%
  569. c | Search, from the left, for the first Ritz |
  570. c | value that has not converged. |
  571. c %-------------------------------------------%
  572. c
  573. leftptr = leftptr + 1
  574. c
  575. else if ( .not. select(rghtptr)) then
  576. c
  577. c %----------------------------------------------%
  578. c | Search, from the right, the first Ritz value |
  579. c | that has converged. |
  580. c %----------------------------------------------%
  581. c
  582. rghtptr = rghtptr - 1
  583. c
  584. else
  585. c
  586. c %----------------------------------------------%
  587. c | Swap the Ritz value on the left that has not |
  588. c | converged with the Ritz value on the right |
  589. c | that has converged. Swap the associated |
  590. c | eigenvector of the tridiagonal matrix H as |
  591. c | well. |
  592. c %----------------------------------------------%
  593. c
  594. temp = workl(ihd+leftptr-1)
  595. workl(ihd+leftptr-1) = workl(ihd+rghtptr-1)
  596. workl(ihd+rghtptr-1) = temp
  597. call dcopy (ncv, workl(iq+ncv*(leftptr-1)), 1,
  598. & workl(iw), 1)
  599. call dcopy (ncv, workl(iq+ncv*(rghtptr-1)), 1,
  600. & workl(iq+ncv*(leftptr-1)), 1)
  601. call dcopy (ncv, workl(iw), 1,
  602. & workl(iq+ncv*(rghtptr-1)), 1)
  603. leftptr = leftptr + 1
  604. rghtptr = rghtptr - 1
  605. c
  606. end if
  607. c
  608. if (leftptr .lt. rghtptr) go to 20
  609. c
  610. end if
  611. c
  612. 30 if (msglvl .gt. 2) then
  613. call dvout (logfil, ncv, workl(ihd), ndigit,
  614. & '_seupd: The eigenvalues of H--reordered')
  615. end if
  616. c
  617. c %----------------------------------------%
  618. c | Load the converged Ritz values into D. |
  619. c %----------------------------------------%
  620. c
  621. call dcopy (nconv, workl(ihd), 1, d, 1)
  622. c
  623. else
  624. c
  625. c %-----------------------------------------------------%
  626. c | Ritz vectors not required. Load Ritz values into D. |
  627. c %-----------------------------------------------------%
  628. c
  629. call dcopy (nconv, workl(ritz), 1, d, 1)
  630. call dcopy (ncv, workl(ritz), 1, workl(ihd), 1)
  631. c
  632. end if
  633. c
  634. c %------------------------------------------------------------------%
  635. c | Transform the Ritz values and possibly vectors and corresponding |
  636. c | Ritz estimates of OP to those of A*x=lambda*B*x. The Ritz values |
  637. c | (and corresponding data) are returned in ascending order. |
  638. c %------------------------------------------------------------------%
  639. c
  640. if (type .eq. 'REGULR') then
  641. c
  642. c %---------------------------------------------------------%
  643. c | Ascending sort of wanted Ritz values, vectors and error |
  644. c | bounds. Not necessary if only Ritz values are desired. |
  645. c %---------------------------------------------------------%
  646. c
  647. if (rvec) then
  648. call dsesrt ('LA', rvec , nconv, d, ncv, workl(iq), ldq)
  649. else
  650. call dcopy (ncv, workl(bounds), 1, workl(ihb), 1)
  651. end if
  652. c
  653. else
  654. c
  655. c %-------------------------------------------------------------%
  656. c | * Make a copy of all the Ritz values. |
  657. c | * Transform the Ritz values back to the original system. |
  658. c | For TYPE = 'SHIFTI' the transformation is |
  659. c | lambda = 1/theta + sigma |
  660. c | For TYPE = 'BUCKLE' the transformation is |
  661. c | lambda = sigma * theta / ( theta - 1 ) |
  662. c | For TYPE = 'CAYLEY' the transformation is |
  663. c | lambda = sigma * (theta + 1) / (theta - 1 ) |
  664. c | where the theta are the Ritz values returned by dsaupd . |
  665. c | NOTES: |
  666. c | *The Ritz vectors are not affected by the transformation. |
  667. c | They are only reordered. |
  668. c %-------------------------------------------------------------%
  669. c
  670. call dcopy (ncv, workl(ihd), 1, workl(iw), 1)
  671. if (type .eq. 'SHIFTI') then
  672. do 40 k=1, ncv
  673. workl(ihd+k-1) = one / workl(ihd+k-1) + sigma
  674. 40 continue
  675. else if (type .eq. 'BUCKLE') then
  676. do 50 k=1, ncv
  677. workl(ihd+k-1) = sigma * workl(ihd+k-1) /
  678. & (workl(ihd+k-1) - one)
  679. 50 continue
  680. else if (type .eq. 'CAYLEY') then
  681. do 60 k=1, ncv
  682. workl(ihd+k-1) = sigma * (workl(ihd+k-1) + one) /
  683. & (workl(ihd+k-1) - one)
  684. 60 continue
  685. end if
  686. c
  687. c %-------------------------------------------------------------%
  688. c | * Store the wanted NCONV lambda values into D. |
  689. c | * Sort the NCONV wanted lambda in WORKL(IHD:IHD+NCONV-1) |
  690. c | into ascending order and apply sort to the NCONV theta |
  691. c | values in the transformed system. We will need this to |
  692. c | compute Ritz estimates in the original system. |
  693. c | * Finally sort the lambda`s into ascending order and apply |
  694. c | to Ritz vectors if wanted. Else just sort lambda`s into |
  695. c | ascending order. |
  696. c | NOTES: |
  697. c | *workl(iw:iw+ncv-1) contain the theta ordered so that they |
  698. c | match the ordering of the lambda. We`ll use them again for |
  699. c | Ritz vector purification. |
  700. c %-------------------------------------------------------------%
  701. c
  702. call dcopy (nconv, workl(ihd), 1, d, 1)
  703. call dsortr ('LA', .true., nconv, workl(ihd), workl(iw))
  704. if (rvec) then
  705. call dsesrt ('LA', rvec , nconv, d, ncv, workl(iq), ldq)
  706. else
  707. call dcopy (ncv, workl(bounds), 1, workl(ihb), 1)
  708. call dscal (ncv, bnorm2/rnorm, workl(ihb), 1)
  709. call dsortr ('LA', .true., nconv, d, workl(ihb))
  710. end if
  711. c
  712. end if
  713. c
  714. c %------------------------------------------------%
  715. c | Compute the Ritz vectors. Transform the wanted |
  716. c | eigenvectors of the symmetric tridiagonal H by |
  717. c | the Lanczos basis matrix V. |
  718. c %------------------------------------------------%
  719. c
  720. if (rvec .and. howmny .eq. 'A') then
  721. c
  722. c %----------------------------------------------------------%
  723. c | Compute the QR factorization of the matrix representing |
  724. c | the wanted invariant subspace located in the first NCONV |
  725. c | columns of workl(iq,ldq). |
  726. c %----------------------------------------------------------%
  727. c
  728. call dgeqr2 (ncv, nconv , workl(iq) ,
  729. & ldq, workl(iw+ncv), workl(ihb),
  730. & ierr)
  731. c
  732. c %--------------------------------------------------------%
  733. c | * Postmultiply V by Q. |
  734. c | * Copy the first NCONV columns of VQ into Z. |
  735. c | The N by NCONV matrix Z is now a matrix representation |
  736. c | of the approximate invariant subspace associated with |
  737. c | the Ritz values in workl(ihd). |
  738. c %--------------------------------------------------------%
  739. c
  740. call dorm2r ('Right', 'Notranspose', n ,
  741. & ncv , nconv , workl(iq),
  742. & ldq , workl(iw+ncv), v ,
  743. & ldv , workd(n+1) , ierr)
  744. call dlacpy ('All', n, nconv, v, ldv, z, ldz)
  745. c
  746. c %-----------------------------------------------------%
  747. c | In order to compute the Ritz estimates for the Ritz |
  748. c | values in both systems, need the last row of the |
  749. c | eigenvector matrix. Remember, it`s in factored form |
  750. c %-----------------------------------------------------%
  751. c
  752. do 65 j = 1, ncv-1
  753. workl(ihb+j-1) = zero
  754. 65 continue
  755. workl(ihb+ncv-1) = one
  756. call dorm2r ('Left', 'Transpose' , ncv ,
  757. & 1 , nconv , workl(iq) ,
  758. & ldq , workl(iw+ncv), workl(ihb),
  759. & ncv , temp , ierr)
  760. c
  761. c %-----------------------------------------------------%
  762. c | Make a copy of the last row into |
  763. c | workl(iw+ncv:iw+2*ncv), as it is needed again in |
  764. c | the Ritz vector purification step below |
  765. c %-----------------------------------------------------%
  766. c
  767. do 67 j = 1, nconv
  768. workl(iw+ncv+j-1) = workl(ihb+j-1)
  769. 67 continue
  770.  
  771. else if (rvec .and. howmny .eq. 'S') then
  772. c
  773. c Not yet implemented. See remark 2 above.
  774. c
  775. end if
  776. c
  777. if (type .eq. 'REGULR' .and. rvec) then
  778. c
  779. do 70 j=1, ncv
  780. workl(ihb+j-1) = rnorm * abs( workl(ihb+j-1) )
  781. 70 continue
  782. c
  783. else if (type .ne. 'REGULR' .and. rvec) then
  784. c
  785. c %-------------------------------------------------%
  786. c | * Determine Ritz estimates of the theta. |
  787. c | If RVEC = .true. then compute Ritz estimates |
  788. c | of the theta. |
  789. c | If RVEC = .false. then copy Ritz estimates |
  790. c | as computed by dsaupd . |
  791. c | * Determine Ritz estimates of the lambda. |
  792. c %-------------------------------------------------%
  793. c
  794. call dscal (ncv, bnorm2, workl(ihb), 1)
  795. if (type .eq. 'SHIFTI') then
  796. c
  797. do 80 k=1, ncv
  798. workl(ihb+k-1) = abs( workl(ihb+k-1) )
  799. & / workl(iw+k-1)**2
  800. 80 continue
  801. c
  802. else if (type .eq. 'BUCKLE') then
  803. c
  804. do 90 k=1, ncv
  805. workl(ihb+k-1) = sigma * abs( workl(ihb+k-1) )
  806. & / (workl(iw+k-1)-one )**2
  807. 90 continue
  808. c
  809. else if (type .eq. 'CAYLEY') then
  810. c
  811. do 100 k=1, ncv
  812. workl(ihb+k-1) = abs( workl(ihb+k-1)
  813. & / workl(iw+k-1)*(workl(iw+k-1)-one) )
  814. 100 continue
  815. c
  816. end if
  817. c
  818. end if
  819. c
  820. if (type .ne. 'REGULR' .and. msglvl .gt. 1) then
  821. call dvout (logfil, nconv, d, ndigit,
  822. & '_seupd: Untransformed converged Ritz values')
  823. call dvout (logfil, nconv, workl(ihb), ndigit,
  824. & '_seupd: Ritz estimates of the untransformed Ritz values')
  825. else if (msglvl .gt. 1) then
  826. call dvout (logfil, nconv, d, ndigit,
  827. & '_seupd: Converged Ritz values')
  828. call dvout (logfil, nconv, workl(ihb), ndigit,
  829. & '_seupd: Associated Ritz estimates')
  830. end if
  831. c
  832. c %-------------------------------------------------%
  833. c | Ritz vector purification step. Formally perform |
  834. c | one of inverse subspace iteration. Only used |
  835. c | for MODE = 3,4,5. See reference 7 |
  836. c %-------------------------------------------------%
  837. c
  838. if (rvec .and. (type .eq. 'SHIFTI' .or. type .eq. 'CAYLEY')) then
  839. c
  840. do 110 k=0, nconv-1
  841. workl(iw+k) = workl(iw+ncv+k)
  842. & / workl(iw+k)
  843. 110 continue
  844. c
  845. else if (rvec .and. type .eq. 'BUCKLE') then
  846. c
  847. do 120 k=0, nconv-1
  848. workl(iw+k) = workl(iw+ncv+k)
  849. & / (workl(iw+k)-one)
  850. 120 continue
  851. c
  852. end if
  853. c
  854. if (rvec .and. type .ne. 'REGULR')
  855. & call dger (n, nconv, one, resid, 1, workl(iw), 1, z, ldz)
  856. c
  857. 9000 continue
  858. c
  859. return
  860. c
  861. c %---------------%
  862. c | End of dseupd |
  863. c %---------------%
  864. c
  865. end
  866.  
  867.  

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