Télécharger dsaupd.eso

Retour à la liste

Numérotation des lignes :

  1. C DSAUPD SOURCE BP208322 15/10/13 21:15:51 8670
  2. c-----------------------------------------------------------------------
  3. c\BeginDoc
  4. c
  5. c\Name: dsaupd
  6. c
  7. c\Description:
  8. c
  9. c Reverse communication interface for the Implicitly Restarted Arnoldi
  10. c Iteration. For symmetric problems this reduces to a variant of the Lanczos
  11. c method. This method has been designed to compute approximations to a
  12. c few eigenpairs of a linear operator OP that is real and symmetric
  13. c with respect to a real positive semi-definite symmetric matrix B,
  14. c i.e.
  15. c
  16. c B*OP = (OP`)*B.
  17. c
  18. c Another way to express this condition is
  19. c
  20. c < x,OPy > = < OPx,y > where < z,w > = z`Bw .
  21. c
  22. c In the standard eigenproblem B is the identity matrix.
  23. c ( A` denotes transpose of A)
  24. c
  25. c The computed approximate eigenvalues are called Ritz values and
  26. c the corresponding approximate eigenvectors are called Ritz vectors.
  27. c
  28. c dsaupd is usually called iteratively to solve one of the
  29. c following problems:
  30. c
  31. c Mode 1: A*x = lambda*x, A symmetric
  32. c ===> OP = A and B = I.
  33. c
  34. c Mode 2: A*x = lambda*M*x, A symmetric, M symmetric positive definite
  35. c ===> OP = inv[M]*A and B = M.
  36. c ===> (If M can be factored see remark 3 below)
  37. c
  38. c Mode 3: K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
  39. c ===> OP = (inv[K - sigma*M])*M and B = M.
  40. c ===> Shift-and-Invert mode
  41. c
  42. c Mode 4: K*x = lambda*KG*x, K symmetric positive semi-definite,
  43. c KG symmetric indefinite
  44. c ===> OP = (inv[K - sigma*KG])*K and B = K.
  45. c ===> Buckling mode
  46. c
  47. c Mode 5: A*x = lambda*M*x, A symmetric, M symmetric positive semi-definite
  48. c ===> OP = inv[A - sigma*M]*[A + sigma*M] and B = M.
  49. c ===> Cayley transformed mode
  50. c
  51. c NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v
  52. c should be accomplished either by a direct method
  53. c using a sparse matrix factorization and solving
  54. c
  55. c [A - sigma*M]*w = v or M*w = v,
  56. c
  57. c or through an iterative method for solving these
  58. c systems. If an iterative method is used, the
  59. c convergence test must be more stringent than
  60. c the accuracy requirements for the eigenvalue
  61. c approximations.
  62. c
  63. c\Usage:
  64. c call dsaupd
  65. c ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
  66. c IPNTR, WORKD, WORKL, LWORKL, INFO )
  67. c
  68. c\Arguments
  69. c IDO Integer. (INPUT/OUTPUT)
  70. c Reverse communication flag. IDO must be zero on the first
  71. c call to dsaupd . IDO will be set internally to
  72. c indicate the type of operation to be performed. Control is
  73. c then given back to the calling routine which has the
  74. c responsibility to carry out the requested operation and call
  75. c dsaupd with the result. The operand is given in
  76. c WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)).
  77. c (If Mode = 2 see remark 5 below)
  78. c -------------------------------------------------------------
  79. c IDO = 0: first call to the reverse communication interface
  80. c IDO = -1: compute Y = OP * X where
  81. c IPNTR(1) is the pointer into WORKD for X,
  82. c IPNTR(2) is the pointer into WORKD for Y.
  83. c This is for the initialization phase to force the
  84. c starting vector into the range of OP.
  85. c IDO = 1: compute Y = OP * X where
  86. c IPNTR(1) is the pointer into WORKD for X,
  87. c IPNTR(2) is the pointer into WORKD for Y.
  88. c In mode 3,4 and 5, the vector B * X is already
  89. c available in WORKD(ipntr(3)). It does not
  90. c need to be recomputed in forming OP * X.
  91. c IDO = 2: compute Y = B * X where
  92. c IPNTR(1) is the pointer into WORKD for X,
  93. c IPNTR(2) is the pointer into WORKD for Y.
  94. c IDO = 3: compute the IPARAM(8) shifts where
  95. c IPNTR(11) is the pointer into WORKL for
  96. c placing the shifts. See remark 6 below.
  97. c IDO = 99: done
  98. c -------------------------------------------------------------
  99. c
  100. c BMAT Character*1. (INPUT)
  101. c BMAT specifies the type of the matrix B that defines the
  102. c semi-inner product for the operator OP.
  103. c B = 'I' -> standard eigenvalue problem A*x = lambda*x
  104. c B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
  105. c
  106. c N Integer. (INPUT)
  107. c Dimension of the eigenproblem.
  108. c
  109. c WHICH Character*2. (INPUT)
  110. c Specify which of the Ritz values of OP to compute.
  111. c
  112. c 'LA' - compute the NEV largest (algebraic) eigenvalues.
  113. c 'SA' - compute the NEV smallest (algebraic) eigenvalues.
  114. c 'LM' - compute the NEV largest (in magnitude) eigenvalues.
  115. c 'SM' - compute the NEV smallest (in magnitude) eigenvalues.
  116. c 'BE' - compute NEV eigenvalues, half from each end of the
  117. c spectrum. When NEV is odd, compute one more from the
  118. c high end than from the low end.
  119. c (see remark 1 below)
  120. c
  121. c NEV Integer. (INPUT)
  122. c Number of eigenvalues of OP to be computed. 0 < NEV < N.
  123. c
  124. c TOL Double precision scalar. (INPUT)
  125. c Stopping criterion: the relative accuracy of the Ritz value
  126. c is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I)).
  127. c If TOL .LE. 0. is passed a default is set:
  128. c DEFAULT = DLAMCH ('EPS') (machine precision as computed
  129. c by the LAPACK auxiliary subroutine DLAMCH ).
  130. c
  131. c RESID Double precision array of length N. (INPUT/OUTPUT)
  132. c On INPUT:
  133. c If INFO .EQ. 0, a random initial residual vector is used.
  134. c If INFO .NE. 0, RESID contains the initial residual vector,
  135. c possibly from a previous run.
  136. c On OUTPUT:
  137. c RESID contains the final residual vector.
  138. c
  139. c NCV Integer. (INPUT)
  140. c Number of columns of the matrix V (less than or equal to N).
  141. c This will indicate how many Lanczos vectors are generated
  142. c at each iteration. After the startup phase in which NEV
  143. c Lanczos vectors are generated, the algorithm generates
  144. c NCV-NEV Lanczos vectors at each subsequent update iteration.
  145. c Most of the cost in generating each Lanczos vector is in the
  146. c matrix-vector product OP*x. (See remark 4 below).
  147. c
  148. c V REAL*8 N by NCV array. (OUTPUT)
  149. c The NCV columns of V contain the Lanczos basis vectors.
  150. c
  151. c LDV Integer. (INPUT)
  152. c Leading dimension of V exactly as declared in the calling
  153. c program.
  154. c
  155. c IPARAM Integer array of length 11. (INPUT/OUTPUT)
  156. c IPARAM(1) = ISHIFT: method for selecting the implicit shifts.
  157. c The shifts selected at each iteration are used to restart
  158. c the Arnoldi iteration in an implicit fashion.
  159. c -------------------------------------------------------------
  160. c ISHIFT = 0: the shifts are provided by the user via
  161. c reverse communication. The NCV eigenvalues of
  162. c the current tridiagonal matrix T are returned in
  163. c the part of WORKL array corresponding to RITZ.
  164. c See remark 6 below.
  165. c ISHIFT = 1: exact shifts with respect to the reduced
  166. c tridiagonal matrix T. This is equivalent to
  167. c restarting the iteration with a starting vector
  168. c that is a linear combination of Ritz vectors
  169. c associated with the "wanted" Ritz values.
  170. c -------------------------------------------------------------
  171. c
  172. c IPARAM(2) = LEVEC
  173. c No longer referenced. See remark 2 below.
  174. c
  175. c IPARAM(3) = MXITER
  176. c On INPUT: maximum number of Arnoldi update iterations allowed.
  177. c On OUTPUT: actual number of Arnoldi update iterations taken.
  178. c
  179. c IPARAM(4) = NB: blocksize to be used in the recurrence.
  180. c The code currently works only for NB = 1.
  181. c
  182. c IPARAM(5) = NCONV: number of "converged" Ritz values.
  183. c This represents the number of Ritz values that satisfy
  184. c the convergence criterion.
  185. c
  186. c IPARAM(6) = IUPD
  187. c No longer referenced. Implicit restarting is ALWAYS used.
  188. c
  189. c IPARAM(7) = MODE
  190. c On INPUT determines what type of eigenproblem is being solved.
  191. c Must be 1,2,3,4,5; See under \Description of dsaupd for the
  192. c five modes available.
  193. c
  194. c IPARAM(8) = NP
  195. c When ido = 3 and the user provides shifts through reverse
  196. c communication (IPARAM(1)=0), dsaupd returns NP, the number
  197. c of shifts the user is to provide. 0 < NP <=NCV-NEV. See Remark
  198. c 6 below.
  199. c
  200. c IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,
  201. c OUTPUT: NUMOP = total number of OP*x operations,
  202. c NUMOPB = total number of B*x operations if BMAT='G',
  203. c NUMREO = total number of steps of re-orthogonalization.
  204. c
  205. c IPNTR Integer array of length 11. (OUTPUT)
  206. c Pointer to mark the starting locations in the WORKD and WORKL
  207. c arrays for matrices/vectors used by the Lanczos iteration.
  208. c -------------------------------------------------------------
  209. c IPNTR(1): pointer to the current operand vector X in WORKD.
  210. c IPNTR(2): pointer to the current result vector Y in WORKD.
  211. c IPNTR(3): pointer to the vector B * X in WORKD when used in
  212. c the shift-and-invert mode.
  213. c IPNTR(4): pointer to the next available location in WORKL
  214. c that is untouched by the program.
  215. c IPNTR(5): pointer to the NCV by 2 tridiagonal matrix T in WORKL.
  216. c IPNTR(6): pointer to the NCV RITZ values array in WORKL.
  217. c IPNTR(7): pointer to the Ritz estimates in array WORKL associated
  218. c with the Ritz values located in RITZ in WORKL.
  219. c IPNTR(11): pointer to the NP shifts in WORKL. See Remark 6 below.
  220. c
  221. c Note: IPNTR(8:10) is only referenced by dseupd . See Remark 2.
  222. c IPNTR(8): pointer to the NCV RITZ values of the original system.
  223. c IPNTR(9): pointer to the NCV corresponding error bounds.
  224. c IPNTR(10): pointer to the NCV by NCV matrix of eigenvectors
  225. c of the tridiagonal matrix T. Only referenced by
  226. c dseupd if RVEC = .TRUE. See Remarks.
  227. c -------------------------------------------------------------
  228. c
  229. c WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION)
  230. c Distributed array to be used in the basic Arnoldi iteration
  231. c for reverse communication. The user should not use WORKD
  232. c as temporary workspace during the iteration. Upon termination
  233. c WORKD(1:N) contains B*RESID(1:N). If the Ritz vectors are desired
  234. c subroutine dseupd uses this output.
  235. c See Data Distribution Note below.
  236. c
  237. c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)
  238. c Private (replicated) array on each PE or array allocated on
  239. c the front end. See Data Distribution Note below.
  240. c
  241. c LWORKL Integer. (INPUT)
  242. c LWORKL must be at least NCV**2 + 8*NCV .
  243. c
  244. c INFO Integer. (INPUT/OUTPUT)
  245. c If INFO .EQ. 0, a randomly initial residual vector is used.
  246. c If INFO .NE. 0, RESID contains the initial residual vector,
  247. c possibly from a previous run.
  248. c Error flag on output.
  249. c = 0: Normal exit.
  250. c = 1: Maximum number of iterations taken.
  251. c All possible eigenvalues of OP has been found. IPARAM(5)
  252. c returns the number of wanted converged Ritz values.
  253. c = 2: No longer an informational error. Deprecated starting
  254. c with release 2 of ARPACK.
  255. c = 3: No shifts could be applied during a cycle of the
  256. c Implicitly restarted Arnoldi iteration. One possibility
  257. c is to increase the size of NCV relative to NEV.
  258. c See remark 4 below.
  259. c = -1: N must be positive.
  260. c = -2: NEV must be positive.
  261. c = -3: NCV must be greater than NEV and less than or equal to N.
  262. c = -4: The maximum number of Arnoldi update iterations allowed
  263. c must be greater than zero.
  264. c = -5: WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.
  265. c = -6: BMAT must be one of 'I' or 'G'.
  266. c = -7: Length of private work array WORKL is not sufficient.
  267. c = -8: Error return from trid. eigenvalue calculation;
  268. c Informatinal error from LAPACK routine dsteqr .
  269. c = -9: Starting vector is zero.
  270. c = -10: IPARAM(7) must be 1,2,3,4,5.
  271. c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
  272. c = -12: IPARAM(1) must be equal to 0 or 1.
  273. c = -13: NEV and WHICH = 'BE' are incompatable.
  274. c = -9999: Could not build an Arnoldi factorization.
  275. c IPARAM(5) returns the size of the current Arnoldi
  276. c factorization. The user is advised to check that
  277. c enough workspace and array storage has been allocated.
  278. c
  279. c
  280. c\Remarks
  281. c 1. The converged Ritz values are always returned in ascending
  282. c algebraic order. The computed Ritz values are approximate
  283. c eigenvalues of OP. The selection of WHICH should be made
  284. c with this in mind when Mode = 3,4,5. After convergence,
  285. c approximate eigenvalues of the original problem may be obtained
  286. c with the ARPACK subroutine dseupd .
  287. c
  288. c 2. If the Ritz vectors corresponding to the converged Ritz values
  289. c are needed, the user must call dseupd immediately following completion
  290. c of dsaupd . This is new starting with version 2.1 of ARPACK.
  291. c
  292. c 3. If M can be factored into a Cholesky factorization M = LL`
  293. c then Mode = 2 should not be selected. Instead one should use
  294. c Mode = 1 with OP = inv(L)*A*inv(L`). Appropriate triangular
  295. c linear systems should be solved with L and L` rather
  296. c than computing inverses. After convergence, an approximate
  297. c eigenvector z of the original problem is recovered by solving
  298. c L`z = x where x is a Ritz vector of OP.
  299. c
  300. c 4. At present there is no a-priori analysis to guide the selection
  301. c of NCV relative to NEV. The only formal requrement is that NCV > NEV.
  302. c However, it is recommended that NCV .ge. 2*NEV. If many problems of
  303. c the same type are to be solved, one should experiment with increasing
  304. c NCV while keeping NEV fixed for a given test problem. This will
  305. c usually decrease the required number of OP*x operations but it
  306. c also increases the work and storage required to maintain the orthogonal
  307. c basis vectors. The optimal "cross-over" with respect to CPU time
  308. c is problem dependent and must be determined empirically.
  309. c
  310. c 5. If IPARAM(7) = 2 then in the Reverse commuication interface the user
  311. c must do the following. When IDO = 1, Y = OP * X is to be computed.
  312. c When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user
  313. c must overwrite X with A*X. Y is then the solution to the linear set
  314. c of equations B*Y = A*X.
  315. c
  316. c 6. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the
  317. c NP = IPARAM(8) shifts in locations:
  318. c 1 WORKL(IPNTR(11))
  319. c 2 WORKL(IPNTR(11)+1)
  320. c .
  321. c .
  322. c .
  323. c NP WORKL(IPNTR(11)+NP-1).
  324. c
  325. c The eigenvalues of the current tridiagonal matrix are located in
  326. c WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1). They are in the
  327. c order defined by WHICH. The associated Ritz estimates are located in
  328. c WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1).
  329. c
  330. c-----------------------------------------------------------------------
  331. c
  332. c\Data Distribution Note:
  333. c
  334. c Fortran-D syntax:
  335. c ================
  336. c REAL RESID(N), V(LDV,NCV), WORKD(3*N), WORKL(LWORKL)
  337. c DECOMPOSE D1(N), D2(N,NCV)
  338. c ALIGN RESID(I) with D1(I)
  339. c ALIGN V(I,J) with D2(I,J)
  340. c ALIGN WORKD(I) with D1(I) range (1:N)
  341. c ALIGN WORKD(I) with D1(I-N) range (N+1:2*N)
  342. c ALIGN WORKD(I) with D1(I-2*N) range (2*N+1:3*N)
  343. c DISTRIBUTE D1(BLOCK), D2(BLOCK,:)
  344. c REPLICATED WORKL(LWORKL)
  345. c
  346. c Cray MPP syntax:
  347. c ===============
  348. c REAL RESID(N), V(LDV,NCV), WORKD(N,3), WORKL(LWORKL)
  349. c SHARED RESID(BLOCK), V(BLOCK,:), WORKD(BLOCK,:)
  350. c REPLICATED WORKL(LWORKL)
  351. c
  352. c
  353. c\BeginLib
  354. c
  355. c\References:
  356. c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
  357. c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
  358. c pp 357-385.
  359. c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
  360. c Restarted Arnoldi Iteration", Rice University Technical Report
  361. c TR95-13, Department of Computational and Applied Mathematics.
  362. c 3. B.N. Parlett, "The Symmetric Eigenvalue Problem". Prentice-Hall,
  363. c 1980.
  364. c 4. B.N. Parlett, B. Nour-Omid, "Towards a Black Box Lanczos Program",
  365. c Computer Physics Communications, 53 (1989), pp 169-179.
  366. c 5. B. Nour-Omid, B.N. Parlett, T. Ericson, P.S. Jensen, "How to
  367. c Implement the Spectral Transformation", Math. Comp., 48 (1987),
  368. c pp 663-673.
  369. c 6. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos
  370. c Algorithm for Solving Sparse Symmetric Generalized Eigenproblems",
  371. c SIAM J. Matr. Anal. Apps., January (1993).
  372. c 7. L. Reichel, W.B. Gragg, "Algorithm 686: FORTRAN Subroutines
  373. c for Updating the QR decomposition", ACM TOMS, December 1990,
  374. c Volume 16 Number 4, pp 369-377.
  375. c 8. R.B. Lehoucq, D.C. Sorensen, "Implementation of Some Spectral
  376. c Transformations in a k-Step Arnoldi Method". In Preparation.
  377. c
  378. c\Routines called:
  379. c dsaup2 ARPACK routine that implements the Implicitly Restarted
  380. c Arnoldi Iteration.
  381. c dstats ARPACK routine that initialize timing and other statistics
  382. c variables.
  383. c ivout ARPACK utility routine that prints integers.
  384. c arscnd ARPACK utility routine for timing.
  385. c dvout ARPACK utility routine that prints vectors.
  386. c dlamch LAPACK routine that determines machine constants.
  387. c
  388. c\Authors
  389. c Danny Sorensen Phuong Vu
  390. c Richard Lehoucq CRPC / Rice University
  391. c Dept. of Computational & Houston, Texas
  392. c Applied Mathematics
  393. c Rice University
  394. c Houston, Texas
  395. c
  396. c\Revision history:
  397. c 12/15/93: Version ' 2.4'
  398. c
  399. c\SCCS Information: @(#)
  400. c FILE: saupd.F SID: 2.8 DATE OF SID: 04/10/01 RELEASE: 2
  401. c
  402. c\Remarks
  403. c 1. None
  404. c
  405. c\EndLib
  406. c
  407. c-----------------------------------------------------------------------
  408. c
  409. subroutine dsaupd
  410. & ( ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam,
  411. & ipntr, workd, workl, lworkl, info )
  412. c
  413. c %----------------------------------------------------%
  414. c | Include files for debugging and timing information |
  415. -INC TARTRAK
  416. c %----------------------------------------------------%
  417. c
  418. c
  419. c %------------------%
  420. c | Scalar Arguments |
  421. c %------------------%
  422. c
  423. character bmat*1, which*2
  424. integer ido, info, ldv, lworkl, n, ncv, nev
  425. REAL*8
  426. & tol
  427. c
  428. c %-----------------%
  429. c | Array Arguments |
  430. c %-----------------%
  431. c
  432. integer iparam(11), ipntr(11)
  433. REAL*8
  434. & resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
  435. c
  436. c %------------%
  437. c | Parameters |
  438. c %------------%
  439. c
  440. REAL*8
  441. & one, zero
  442. parameter (one = 1.0D+0 , zero = 0.0D+0 )
  443. c
  444. c %---------------%
  445. c | Local Scalars |
  446. c %---------------%
  447. c
  448. integer bounds, ierr, ih, iq, ishift, iupd, iw,
  449. & ldh, ldq, msglvl, mxiter, mode, nb,
  450. & nev0, next, np, ritz, j
  451. save bounds, ierr, ih, iq, ishift, iupd, iw,
  452. & ldh, ldq, msglvl, mxiter, mode, nb,
  453. & nev0, next, np, ritz
  454. c
  455. c %----------------------%
  456. c | External Subroutines |
  457. c %----------------------%
  458. c
  459. external dsaup2 , dvout , ivout, arscnd, dstats
  460. c
  461. c %--------------------%
  462. c | External Functions |
  463. c %--------------------%
  464. c
  465. REAL*8
  466. external dlamch
  467. c
  468. c %-----------------------%
  469. c | Executable Statements |
  470. c %-----------------------%
  471. c
  472. if (ido .eq. 0) then
  473. c
  474. c %-------------------------------%
  475. c | Initialize timing statistics |
  476. c | & message level for debugging |
  477. c %-------------------------------%
  478. c
  479. call dstats
  480. * call arscnd (t0)
  481. msglvl = msaupd
  482. c
  483. ierr = 0
  484. ishift = iparam(1)
  485. mxiter = iparam(3)
  486. c nb = iparam(4)
  487. nb = 1
  488. c
  489. c %--------------------------------------------%
  490. c | Revision 2 performs only implicit restart. |
  491. c %--------------------------------------------%
  492. c
  493. iupd = 1
  494. mode = iparam(7)
  495. c
  496. c %----------------%
  497. c | Error checking |
  498. c %----------------%
  499. c
  500. if (n .le. 0) then
  501. ierr = -1
  502. else if (nev .le. 0) then
  503. ierr = -2
  504. else if (ncv .le. nev .or. ncv .gt. n) then
  505. ierr = -3
  506. end if
  507. c
  508. c %----------------------------------------------%
  509. c | NP is the number of additional steps to |
  510. c | extend the length NEV Lanczos factorization. |
  511. c %----------------------------------------------%
  512. c
  513. np = ncv - nev
  514. c
  515. if (mxiter .le. 0) ierr = -4
  516. if (which .ne. 'LM' .and.
  517. & which .ne. 'SM' .and.
  518. & which .ne. 'LA' .and.
  519. & which .ne. 'SA' .and.
  520. & which .ne. 'BE') ierr = -5
  521. if (bmat .ne. 'I' .and. bmat .ne. 'G') ierr = -6
  522. c
  523. if (lworkl .lt. ncv**2 + 8*ncv) ierr = -7
  524. if (mode .lt. 1 .or. mode .gt. 5) then
  525. ierr = -10
  526. else if (mode .eq. 1 .and. bmat .eq. 'G') then
  527. ierr = -11
  528. else if (ishift .lt. 0 .or. ishift .gt. 1) then
  529. ierr = -12
  530. else if (nev .eq. 1 .and. which .eq. 'BE') then
  531. ierr = -13
  532. end if
  533. c
  534. c %------------%
  535. c | Error Exit |
  536. c %------------%
  537. c
  538. if (ierr .ne. 0) then
  539. info = ierr
  540. ido = 99
  541. go to 9000
  542. end if
  543. c
  544. c %------------------------%
  545. c | Set default parameters |
  546. c %------------------------%
  547. c
  548. if (nb .le. 0) nb = 1
  549. if (tol .le. zero) tol = dlamch ('EpsMach')
  550. c
  551. c %----------------------------------------------%
  552. c | NP is the number of additional steps to |
  553. c | extend the length NEV Lanczos factorization. |
  554. c | NEV0 is the local variable designating the |
  555. c | size of the invariant subspace desired. |
  556. c %----------------------------------------------%
  557. c
  558. np = ncv - nev
  559. nev0 = nev
  560. c
  561. c %-----------------------------%
  562. c | Zero out internal workspace |
  563. c %-----------------------------%
  564. c
  565. do 10 j = 1, ncv**2 + 8*ncv
  566. workl(j) = zero
  567. 10 continue
  568. c
  569. c %-------------------------------------------------------%
  570. c | Pointer into WORKL for address of H, RITZ, BOUNDS, Q |
  571. c | etc... and the remaining workspace. |
  572. c | Also update pointer to be used on output. |
  573. c | Memory is laid out as follows: |
  574. c | workl(1:2*ncv) := generated tridiagonal matrix |
  575. c | workl(2*ncv+1:2*ncv+ncv) := ritz values |
  576. c | workl(3*ncv+1:3*ncv+ncv) := computed error bounds |
  577. c | workl(4*ncv+1:4*ncv+ncv*ncv) := rotation matrix Q |
  578. c | workl(4*ncv+ncv*ncv+1:7*ncv+ncv*ncv) := workspace |
  579. c %-------------------------------------------------------%
  580. c
  581. ldh = ncv
  582. ldq = ncv
  583. ih = 1
  584. ritz = ih + 2*ldh
  585. bounds = ritz + ncv
  586. iq = bounds + ncv
  587. iw = iq + ncv**2
  588. next = iw + 3*ncv
  589. c
  590. ipntr(4) = next
  591. ipntr(5) = ih
  592. ipntr(6) = ritz
  593. ipntr(7) = bounds
  594. ipntr(11) = iw
  595. end if
  596. c
  597. c %-------------------------------------------------------%
  598. c | Carry out the Implicitly restarted Lanczos Iteration. |
  599. c %-------------------------------------------------------%
  600. c
  601. call dsaup2
  602. & ( ido, bmat, n, which, nev0, np, tol, resid, mode, iupd,
  603. & ishift, mxiter, v, ldv, workl(ih), ldh, workl(ritz),
  604. & workl(bounds), workl(iq), ldq, workl(iw), ipntr, workd,
  605. & info )
  606. c
  607. c %--------------------------------------------------%
  608. c | ido .ne. 99 implies use of reverse communication |
  609. c | to compute operations involving OP or shifts. |
  610. c %--------------------------------------------------%
  611. c
  612. if (ido .eq. 3) iparam(8) = np
  613. if (ido .ne. 99) go to 9000
  614. c
  615. iparam(3) = mxiter
  616. iparam(5) = np
  617. iparam(9) = nopx
  618. iparam(10) = nbx
  619. iparam(11) = nrorth
  620. c
  621. c %------------------------------------%
  622. c | Exit if there was an informational |
  623. c | error within dsaup2 . |
  624. c %------------------------------------%
  625. c
  626. if (info .lt. 0) go to 9000
  627. if (info .eq. 2) info = 3
  628. c
  629. if (msglvl .gt. 0) then
  630. call ivout (logfil, 1, mxiter, ndigit,
  631. & '_saupd: number of update iterations taken')
  632. call ivout (logfil, 1, np, ndigit,
  633. & '_saupd: number of "converged" Ritz values')
  634. call dvout (logfil, np, workl(Ritz), ndigit,
  635. & '_saupd: final Ritz values')
  636. call dvout (logfil, np, workl(Bounds), ndigit,
  637. & '_saupd: corresponding error bounds')
  638. end if
  639. c
  640. * call arscnd (t1)
  641. tsaupd = t1 - t0
  642. c
  643. if (msglvl .gt. 0) then
  644. c
  645. c %--------------------------------------------------------%
  646. c | Version Number & Version Date are defined in version.h |
  647. c %--------------------------------------------------------%
  648. c
  649. write (6,1000)
  650. write (6,1100) mxiter, nopx, nbx, nrorth, nitref, nrstrt,
  651. & tmvopx, tmvbx, tsaupd, tsaup2, tsaitr, titref,
  652. & tgetv0, tseigt, tsgets, tsapps, tsconv
  653. 1000 format (//,
  654. & 5x, '==========================================',/
  655. & 5x, '= Symmetric implicit Arnoldi update code =',/
  656. & 5x, '= Version Number:', ' 2.4' , 19x, ' =',/
  657. & 5x, '= Version Date: ', ' 07/31/96' , 14x, ' =',/
  658. & 5x, '==========================================',/
  659. & 5x, '= Summary of timing statistics =',/
  660. & 5x, '==========================================',//)
  661. 1100 format (
  662. & 5x, 'Total number update iterations = ', i5,/
  663. & 5x, 'Total number of OP*x operations = ', i5,/
  664. & 5x, 'Total number of B*x operations = ', i5,/
  665. & 5x, 'Total number of reorthogonalization steps = ', i5,/
  666. & 5x, 'Total number of iterative refinement steps = ', i5,/
  667. & 5x, 'Total number of restart steps = ', i5,/
  668. & 5x, 'Total time in user OP*x operation = ', f12.6,/
  669. & 5x, 'Total time in user B*x operation = ', f12.6,/
  670. & 5x, 'Total time in Arnoldi update routine = ', f12.6,/
  671. & 5x, 'Total time in saup2 routine = ', f12.6,/
  672. & 5x, 'Total time in basic Arnoldi iteration loop = ', f12.6,/
  673. & 5x, 'Total time in reorthogonalization phase = ', f12.6,/
  674. & 5x, 'Total time in (re)start vector generation = ', f12.6,/
  675. & 5x, 'Total time in trid eigenvalue subproblem = ', f12.6,/
  676. & 5x, 'Total time in getting the shifts = ', f12.6,/
  677. & 5x, 'Total time in applying the shifts = ', f12.6,/
  678. & 5x, 'Total time in convergence testing = ', f12.6)
  679. end if
  680. c
  681. 9000 continue
  682. c
  683. return
  684. c
  685. c %---------------%
  686. c | End of dsaupd |
  687. c %---------------%
  688. c
  689. end
  690.  
  691.  

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