Télécharger dnaupd.eso

Retour à la liste

Numérotation des lignes :

dnaupd
  1. C DNAUPD SOURCE FANDEUR 22/05/02 21:15:11 11359
  2. c\BeginDoc
  3. c
  4. c\Name: dnaupd
  5. c
  6. c\Description:
  7. c Reverse communication interface for the Implicitly Restarted Arnoldi
  8. c iteration. This subroutine computes approximations to a few eigenpairs
  9. c of a linear operator "OP" with respect to a semi-inner product defined by
  10. c a symmetric positive semi-definite real matrix B. B may be the identity
  11. c matrix. NOTE: If the linear operator "OP" is real and symmetric
  12. c with respect to the real positive semi-definite symmetric matrix B,
  13. c i.e. B*OP = (OP`)*B, then subroutine dsaupd should be used instead.
  14. c
  15. c The computed approximate eigenvalues are called Ritz values and
  16. c the corresponding approximate eigenvectors are called Ritz vectors.
  17. c
  18. c dnaupd is usually called iteratively to solve one of the
  19. c following problems:
  20. c
  21. c Mode 1: A*x = lambda*x.
  22. c ===> OP = A and B = I.
  23. c
  24. c Mode 2: A*x = lambda*M*x, M symmetric positive definite
  25. c ===> OP = inv[M]*A and B = M.
  26. c ===> (If M can be factored see remark 3 below)
  27. c
  28. c Mode 3: A*x = lambda*M*x, M symmetric semi-definite
  29. c ===> OP = Real_Part{ inv[A - sigma*M]*M } and B = M.
  30. c ===> shift-and-invert mode (in real arithmetic)
  31. c If OP*x = amu*x, then
  32. c amu = 1/2 * [ 1/(lambda-sigma) + 1/(lambda-conjg(sigma)) ].
  33. c Note: If sigma is real, i.e. imaginary part of sigma is zero;
  34. c Real_Part{ inv[A - sigma*M]*M } == inv[A - sigma*M]*M
  35. c amu == 1/(lambda-sigma).
  36. c
  37. c Mode 4: A*x = lambda*M*x, M symmetric semi-definite
  38. c ===> OP = Imaginary_Part{ inv[A - sigma*M]*M } and B = M.
  39. c ===> shift-and-invert mode (in real arithmetic)
  40. c If OP*x = amu*x, then
  41. c amu = 1/2i * [ 1/(lambda-sigma) - 1/(lambda-conjg(sigma)) ].
  42. c
  43. c Both mode 3 and 4 give the same enhancement to eigenvalues close to
  44. c the (complex) shift sigma. However, as lambda goes to infinity,
  45. c the operator OP in mode 4 dampens the eigenvalues more strongly than
  46. c does OP defined in mode 3.
  47. c
  48. c NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v
  49. c should be accomplished either by a direct method
  50. c using a sparse matrix factorization and solving
  51. c
  52. c [A - sigma*M]*w = v or M*w = v,
  53. c
  54. c or through an iterative method for solving these
  55. c systems. If an iterative method is used, the
  56. c convergence test must be more stringent than
  57. c the accuracy requirements for the eigenvalue
  58. c approximations.
  59. c
  60. c\Usage:
  61. c call dnaupd
  62. c ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
  63. c IPNTR, WORKD, WORKL, LWORKL, INFO )
  64. c
  65. c\Arguments
  66. c IDO Integer. (INPUT/OUTPUT)
  67. c Reverse communication flag. IDO must be zero on the first
  68. c call to dnaupd . IDO will be set internally to
  69. c indicate the type of operation to be performed. Control is
  70. c then given back to the calling routine which has the
  71. c responsibility to carry out the requested operation and call
  72. c dnaupd with the result. The operand is given in
  73. c WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)).
  74. c -------------------------------------------------------------
  75. c IDO = 0: first call to the reverse communication interface
  76. c IDO = -1: compute Y = OP * X where
  77. c IPNTR(1) is the pointer into WORKD for X,
  78. c IPNTR(2) is the pointer into WORKD for Y.
  79. c This is for the initialization phase to force the
  80. c starting vector into the range of OP.
  81. c IDO = 1: compute Y = OP * X where
  82. c IPNTR(1) is the pointer into WORKD for X,
  83. c IPNTR(2) is the pointer into WORKD for Y.
  84. c In mode 3 and 4, the vector B * X is already
  85. c available in WORKD(ipntr(3)). It does not
  86. c need to be recomputed in forming OP * X.
  87. c IDO = 2: compute Y = B * X where
  88. c IPNTR(1) is the pointer into WORKD for X,
  89. c IPNTR(2) is the pointer into WORKD for Y.
  90. c IDO = 3: compute the IPARAM(8) real and imaginary parts
  91. c of the shifts where INPTR(14) is the pointer
  92. c into WORKL for placing the shifts. See Remark
  93. c 5 below.
  94. c IDO = 99: done
  95. c -------------------------------------------------------------
  96. c
  97. c BMAT Character*1. (INPUT)
  98. c BMAT specifies the type of the matrix B that defines the
  99. c semi-inner product for the operator OP.
  100. c BMAT = 'I' -> standard eigenvalue problem A*x = lambda*x
  101. c BMAT = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
  102. c
  103. c N Integer. (INPUT)
  104. c Dimension of the eigenproblem.
  105. c
  106. c WHICH Character*2. (INPUT)
  107. c 'LM' -> want the NEV eigenvalues of largest magnitude.
  108. c 'SM' -> want the NEV eigenvalues of smallest magnitude.
  109. c 'LR' -> want the NEV eigenvalues of largest real part.
  110. c 'SR' -> want the NEV eigenvalues of smallest real part.
  111. c 'LI' -> want the NEV eigenvalues of largest imaginary part.
  112. c 'SI' -> want the NEV eigenvalues of smallest imaginary part.
  113. c
  114. c NEV Integer. (INPUT)
  115. c Number of eigenvalues of OP to be computed. 0 < NEV < N-1.
  116. c
  117. c TOL Double precision scalar. (INPUT/OUTPUT)
  118. c Stopping criterion: the relative accuracy of the Ritz value
  119. c is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I))
  120. c where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex.
  121. c DEFAULT = DLAMCH ('EPS') (machine precision as computed
  122. c by the LAPACK auxiliary subroutine DLAMCH ).
  123. c
  124. c RESID Double precision array of length N. (INPUT/OUTPUT)
  125. c On INPUT:
  126. c If INFO .EQ. 0, a random initial residual vector is used.
  127. c If INFO .NE. 0, RESID contains the initial residual vector,
  128. c possibly from a previous run.
  129. c On OUTPUT:
  130. c RESID contains the final residual vector.
  131. c
  132. c NCV Integer. (INPUT)
  133. c Number of columns of the matrix V. NCV must satisfy the two
  134. c inequalities 2 <= NCV-NEV and NCV <= N.
  135. c This will indicate how many Arnoldi vectors are generated
  136. c at each iteration. After the startup phase in which NEV
  137. c Arnoldi vectors are generated, the algorithm generates
  138. c approximately NCV-NEV Arnoldi vectors at each subsequent update
  139. c iteration. Most of the cost in generating each Arnoldi vector is
  140. c in the matrix-vector operation OP*x.
  141. c NOTE: 2 <= NCV-NEV in order that complex conjugate pairs of Ritz
  142. c values are kept together. (See remark 4 below)
  143. c
  144. c V REAL*8 array N by NCV. (OUTPUT)
  145. c Contains the final set of Arnoldi basis vectors.
  146. c
  147. c LDV Integer. (INPUT)
  148. c Leading dimension of V exactly as declared in the calling program.
  149. c
  150. c IPARAM Integer array of length 11. (INPUT/OUTPUT)
  151. c IPARAM(1) = ISHIFT: method for selecting the implicit shifts.
  152. c The shifts selected at each iteration are used to restart
  153. c the Arnoldi iteration in an implicit fashion.
  154. c -------------------------------------------------------------
  155. c ISHIFT = 0: the shifts are provided by the user via
  156. c reverse communication. The real and imaginary
  157. c parts of the NCV eigenvalues of the Hessenberg
  158. c matrix H are returned in the part of the WORKL
  159. c array corresponding to RITZR and RITZI. See remark
  160. c 5 below.
  161. c ISHIFT = 1: exact shifts with respect to the current
  162. c Hessenberg matrix H. This is equivalent to
  163. c restarting the iteration with a starting vector
  164. c that is a linear combination of approximate Schur
  165. c vectors associated with the "wanted" Ritz values.
  166. c -------------------------------------------------------------
  167. c
  168. c IPARAM(2) = No longer referenced.
  169. c
  170. c IPARAM(3) = MXITER
  171. c On INPUT: maximum number of Arnoldi update iterations allowed.
  172. c On OUTPUT: actual number of Arnoldi update iterations taken.
  173. c
  174. c IPARAM(4) = NB: blocksize to be used in the recurrence.
  175. c The code currently works only for NB = 1.
  176. c
  177. c IPARAM(5) = NCONV: number of "converged" Ritz values.
  178. c This represents the number of Ritz values that satisfy
  179. c the convergence criterion.
  180. c
  181. c IPARAM(6) = IUPD
  182. c No longer referenced. Implicit restarting is ALWAYS used.
  183. c
  184. c IPARAM(7) = MODE
  185. c On INPUT determines what type of eigenproblem is being solved.
  186. c Must be 1,2,3,4; See under \Description of dnaupd for the
  187. c four modes available.
  188. c
  189. c IPARAM(8) = NP
  190. c When ido = 3 and the user provides shifts through reverse
  191. c communication (IPARAM(1)=0), dnaupd returns NP, the number
  192. c of shifts the user is to provide. 0 < NP <=NCV-NEV. See Remark
  193. c 5 below.
  194. c
  195. c IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,
  196. c OUTPUT: NUMOP = total number of OP*x operations,
  197. c NUMOPB = total number of B*x operations if BMAT='G',
  198. c NUMREO = total number of steps of re-orthogonalization.
  199. c
  200. c IPNTR Integer array of length 14. (OUTPUT)
  201. c Pointer to mark the starting locations in the WORKD and WORKL
  202. c arrays for matrices/vectors used by the Arnoldi iteration.
  203. c -------------------------------------------------------------
  204. c IPNTR(1): pointer to the current operand vector X in WORKD.
  205. c IPNTR(2): pointer to the current result vector Y in WORKD.
  206. c IPNTR(3): pointer to the vector B * X in WORKD when used in
  207. c the shift-and-invert mode.
  208. c IPNTR(4): pointer to the next available location in WORKL
  209. c that is untouched by the program.
  210. c IPNTR(5): pointer to the NCV by NCV upper Hessenberg matrix
  211. c H in WORKL.
  212. c IPNTR(6): pointer to the real part of the ritz value array
  213. c RITZR in WORKL.
  214. c IPNTR(7): pointer to the imaginary part of the ritz value array
  215. c RITZI in WORKL.
  216. c IPNTR(8): pointer to the Ritz estimates in array WORKL associated
  217. c with the Ritz values located in RITZR and RITZI in WORKL.
  218. c
  219. c IPNTR(14): pointer to the NP shifts in WORKL. See Remark 5 below.
  220. c
  221. c Note: IPNTR(9:13) is only referenced by dneupd . See Remark 2 below.
  222. c
  223. c IPNTR(9): pointer to the real part of the NCV RITZ values of the
  224. c original system.
  225. c IPNTR(10): pointer to the imaginary part of the NCV RITZ values of
  226. c the original system.
  227. c IPNTR(11): pointer to the NCV corresponding error bounds.
  228. c IPNTR(12): pointer to the NCV by NCV upper quasi-triangular
  229. c Schur matrix for H.
  230. c IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors
  231. c of the upper Hessenberg matrix H. Only referenced by
  232. c dneupd if RVEC = .TRUE. See Remark 2 below.
  233. c -------------------------------------------------------------
  234. c
  235. c WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION)
  236. c Distributed array to be used in the basic Arnoldi iteration
  237. c for reverse communication. The user should not use WORKD
  238. c as temporary workspace during the iteration. Upon termination
  239. c WORKD(1:N) contains B*RESID(1:N). If an invariant subspace
  240. c associated with the converged Ritz values is desired, see remark
  241. c 2 below, subroutine dneupd uses this output.
  242. c See Data Distribution Note below.
  243. c
  244. c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)
  245. c Private (replicated) array on each PE or array allocated on
  246. c the front end. See Data Distribution Note below.
  247. c
  248. c LWORKL Integer. (INPUT)
  249. c LWORKL must be at least 3*NCV**2 + 6*NCV.
  250. c
  251. c INFO Integer. (INPUT/OUTPUT)
  252. c If INFO .EQ. 0, a randomly initial residual vector is used.
  253. c If INFO .NE. 0, RESID contains the initial residual vector,
  254. c possibly from a previous run.
  255. c Error flag on output.
  256. c = 0: Normal exit.
  257. c = 1: Maximum number of iterations taken.
  258. c All possible eigenvalues of OP has been found. IPARAM(5)
  259. c returns the number of wanted converged Ritz values.
  260. c = 2: No longer an informational error. Deprecated starting
  261. c with release 2 of ARPACK.
  262. c = 3: No shifts could be applied during a cycle of the
  263. c Implicitly restarted Arnoldi iteration. One possibility
  264. c is to increase the size of NCV relative to NEV.
  265. c See remark 4 below.
  266. c = -1: N must be positive.
  267. c = -2: NEV must be positive.
  268. c = -3: NCV-NEV >= 2 and less than or equal to N.
  269. c = -4: The maximum number of Arnoldi update iteration
  270. c must be greater than zero.
  271. c = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
  272. c = -6: BMAT must be one of 'I' or 'G'.
  273. c = -7: Length of private work array is not sufficient.
  274. c = -8: Error return from LAPACK eigenvalue calculation;
  275. c = -9: Starting vector is zero.
  276. c = -10: IPARAM(7) must be 1,2,3,4.
  277. c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
  278. c = -12: IPARAM(1) must be equal to 0 or 1.
  279. c = -9999: Could not build an Arnoldi factorization.
  280. c IPARAM(5) returns the size of the current Arnoldi
  281. c factorization.
  282. c
  283. c\Remarks
  284. c 1. The computed Ritz values are approximate eigenvalues of OP. The
  285. c selection of WHICH should be made with this in mind when
  286. c Mode = 3 and 4. After convergence, approximate eigenvalues of the
  287. c original problem may be obtained with the ARPACK subroutine dneupd .
  288. c
  289. c 2. If a basis for the invariant subspace corresponding to the converged Ritz
  290. c values is needed, the user must call dneupd immediately following
  291. c completion of dnaupd . This is new starting with release 2 of ARPACK.
  292. c
  293. c 3. If M can be factored into a Cholesky factorization M = LL`
  294. c then Mode = 2 should not be selected. Instead one should use
  295. c Mode = 1 with OP = inv(L)*A*inv(L`). Appropriate triangular
  296. c linear systems should be solved with L and L` rather
  297. c than computing inverses. After convergence, an approximate
  298. c eigenvector z of the original problem is recovered by solving
  299. c L`z = x where x is a Ritz vector of OP.
  300. c
  301. c 4. At present there is no a-priori analysis to guide the selection
  302. c of NCV relative to NEV. The only formal requrement is that NCV > NEV + 2.
  303. c However, it is recommended that NCV .ge. 2*NEV+1. If many problems of
  304. c the same type are to be solved, one should experiment with increasing
  305. c NCV while keeping NEV fixed for a given test problem. This will
  306. c usually decrease the required number of OP*x operations but it
  307. c also increases the work and storage required to maintain the orthogonal
  308. c basis vectors. The optimal "cross-over" with respect to CPU time
  309. c is problem dependent and must be determined empirically.
  310. c See Chapter 8 of Reference 2 for further information.
  311. c
  312. c 5. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the
  313. c NP = IPARAM(8) real and imaginary parts of the shifts in locations
  314. c real part imaginary part
  315. c ----------------------- --------------
  316. c 1 WORKL(IPNTR(14)) WORKL(IPNTR(14)+NP)
  317. c 2 WORKL(IPNTR(14)+1) WORKL(IPNTR(14)+NP+1)
  318. c . .
  319. c . .
  320. c . .
  321. c NP WORKL(IPNTR(14)+NP-1) WORKL(IPNTR(14)+2*NP-1).
  322. c
  323. c Only complex conjugate pairs of shifts may be applied and the pairs
  324. c must be placed in consecutive locations. The real part of the
  325. c eigenvalues of the current upper Hessenberg matrix are located in
  326. c WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1) and the imaginary part
  327. c in WORKL(IPNTR(7)) through WORKL(IPNTR(7)+NCV-1). They are ordered
  328. c according to the order defined by WHICH. The complex conjugate
  329. c pairs are kept together and the associated Ritz estimates are located in
  330. c WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1).
  331. c
  332. c-----------------------------------------------------------------------
  333. c
  334. c\Data Distribution Note:
  335. c
  336. c Fortran-D syntax:
  337. c ================
  338. c Double precision resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
  339. c decompose d1(n), d2(n,ncv)
  340. c align resid(i) with d1(i)
  341. c align v(i,j) with d2(i,j)
  342. c align workd(i) with d1(i) range (1:n)
  343. c align workd(i) with d1(i-n) range (n+1:2*n)
  344. c align workd(i) with d1(i-2*n) range (2*n+1:3*n)
  345. c distribute d1(block), d2(block,:)
  346. c replicated workl(lworkl)
  347. c
  348. c Cray MPP syntax:
  349. c ===============
  350. c Double precision resid(n), v(ldv,ncv), workd(n,3), workl(lworkl)
  351. c shared resid(block), v(block,:), workd(block,:)
  352. c replicated workl(lworkl)
  353. c
  354. c CM2/CM5 syntax:
  355. c ==============
  356. c
  357. c-----------------------------------------------------------------------
  358. c
  359. c include 'ex-nonsym.doc'
  360. c
  361. c-----------------------------------------------------------------------
  362. c
  363. c\BeginLib
  364. c
  365. c\Local variables:
  366. c xxxxxx real
  367. c
  368. c\References:
  369. c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
  370. c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
  371. c pp 357-385.
  372. c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
  373. c Restarted Arnoldi Iteration", Rice University Technical Report
  374. c TR95-13, Department of Computational and Applied Mathematics.
  375. c 3. B.N. Parlett & Y. Saad, "Complex Shift and Invert Strategies for
  376. c Real Matrices", Linear Algebra and its Applications, vol 88/89,
  377. c pp 575-595, (1987).
  378. c
  379. c\Routines called:
  380. c dnaup2 ARPACK routine that implements the Implicitly Restarted
  381. c Arnoldi Iteration.
  382. c ivout ARPACK utility routine that prints integers.
  383. c arscnd ARPACK utility routine for timing. -> deleted by BP in 2020
  384. c dvout ARPACK utility routine that prints vectors.
  385. c dlamch LAPACK routine that determines machine constants.
  386. c
  387. c\Author
  388. c Danny Sorensen Phuong Vu
  389. c Richard Lehoucq CRPC / Rice University
  390. c Dept. of Computational & Houston, Texas
  391. c Applied Mathematics
  392. c Rice University
  393. c Houston, Texas
  394. c
  395. c\Revision history:
  396. c 12/16/93: Version '1.1'
  397. c
  398. c\SCCS Information: @(#)
  399. c FILE: naupd.F SID: 2.8 DATE OF SID: 04/10/01 RELEASE: 2
  400. c
  401. c\Remarks
  402. c
  403. c\EndLib
  404. c
  405. c-----------------------------------------------------------------------
  406. c
  407. subroutine dnaupd
  408. & ( ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam,
  409. & ipntr, workd, workl, lworkl, ITRAK, info )
  410. c
  411. c %----------------------------------------------------%
  412. c | Include files for debugging and timing information |
  413. c -INC TARTRAK
  414. c %----------------------------------------------------%
  415. c
  416. c
  417. c %------------------%
  418. c | Scalar Arguments |
  419. c %------------------%
  420. c
  421. character bmat*1, which*2
  422. integer ido, info, ldv, lworkl, n, ncv, nev
  423. REAL*8
  424. & tol
  425. c
  426. c %-----------------%
  427. c | Array Arguments |
  428. c %-----------------%
  429. c
  430. integer iparam(11), ipntr(14)
  431. INTEGER ITRAK(5)
  432. REAL*8
  433. & resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
  434. c
  435. c %------------%
  436. c | Parameters |
  437. c %------------%
  438. c
  439. REAL*8
  440. & one, zero
  441. parameter (one = 1.0D+0 , zero = 0.0D+0 )
  442. c
  443. c %---------------%
  444. c | Local Scalars |
  445. c %---------------%
  446. c
  447. integer bounds, ierr, ih, iq, ishift, iupd, iw,
  448. & ldh, ldq, levec, mode, msglvl, mxiter, nb,
  449. & nev0, next, np, ritzi, ritzr, j
  450. save bounds, ih, iq, ishift, iupd, iw, ldh, ldq,
  451. c & levec, mode, msglvl, mxiter, nb, nev0, next,
  452. & levec, mode, mxiter, nb, nev0, next,
  453. & np, ritzi, ritzr
  454. parameter (msglvl=0)
  455. c
  456. c %----------------------%
  457. c | External Subroutines |
  458. c %----------------------%
  459. c
  460. external dnaup2 , dvout , ivout, dstatn
  461. c
  462. c %--------------------%
  463. c | External Functions |
  464. c %--------------------%
  465. c
  466. REAL*8
  467. external dlamch
  468. c
  469. c %-----------------------%
  470. c | Executable Statements |
  471. c %-----------------------%
  472. nopx =ITRAK(1)
  473. nbx =ITRAK(2)
  474. nrorth=ITRAK(3)
  475. nitref=ITRAK(4)
  476. nrstrt=ITRAK(5)
  477. c
  478. if (ido .eq. 0) then
  479. c
  480. c %-------------------------------%
  481. c | Initialize timing statistics |
  482. c | & message level for debugging |
  483. c %-------------------------------%
  484. c
  485. call dstatn(nopx,nbx,nrorth,nitref,nrstrt)
  486. * call arscnd (t0)
  487. c msglvl = mnaupd
  488. c
  489. c %----------------%
  490. c | Error checking |
  491. c %----------------%
  492. c
  493. ierr = 0
  494. ishift = iparam(1)
  495. c levec = iparam(2)
  496. mxiter = iparam(3)
  497. c nb = iparam(4)
  498. nb = 1
  499. c
  500. c %--------------------------------------------%
  501. c | Revision 2 performs only implicit restart. |
  502. c %--------------------------------------------%
  503. c
  504. iupd = 1
  505. mode = iparam(7)
  506. c
  507. if (n .le. 0) then
  508. ierr = -1
  509. else if (nev .le. 0) then
  510. ierr = -2
  511. else if (ncv .le. nev+1 .or. ncv .gt. n) then
  512. ierr = -3
  513. else if (mxiter .le. 0) then
  514. ierr = -4
  515. else if (which .ne. 'LM' .and.
  516. & which .ne. 'SM' .and.
  517. & which .ne. 'LR' .and.
  518. & which .ne. 'SR' .and.
  519. & which .ne. 'LI' .and.
  520. & which .ne. 'SI') then
  521. ierr = -5
  522. else if (bmat .ne. 'I' .and. bmat .ne. 'G') then
  523. ierr = -6
  524. else if (lworkl .lt. 3*ncv**2 + 6*ncv) then
  525. ierr = -7
  526. else if (mode .lt. 1 .or. mode .gt. 4) then
  527. ierr = -10
  528. else if (mode .eq. 1 .and. bmat .eq. 'G') then
  529. ierr = -11
  530. else if (ishift .lt. 0 .or. ishift .gt. 1) then
  531. ierr = -12
  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, 3*ncv**2 + 6*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:ncv*ncv) := generated Hessenberg matrix |
  575. c | workl(ncv*ncv+1:ncv*ncv+2*ncv) := real and imaginary |
  576. c | parts of ritz values |
  577. c | workl(ncv*ncv+2*ncv+1:ncv*ncv+3*ncv) := error bounds |
  578. c | workl(ncv*ncv+3*ncv+1:2*ncv*ncv+3*ncv) := rotation matrix Q |
  579. c | workl(2*ncv*ncv+3*ncv+1:3*ncv*ncv+6*ncv) := workspace |
  580. c | The final workspace is needed by subroutine dneigh called |
  581. c | by dnaup2 . Subroutine dneigh calls LAPACK routines for |
  582. c | calculating eigenvalues and the last row of the eigenvector |
  583. c | matrix. |
  584. c %-------------------------------------------------------------%
  585. c
  586. ldh = ncv
  587. ldq = ncv
  588. ih = 1
  589. ritzr = ih + ldh*ncv
  590. ritzi = ritzr + ncv
  591. bounds = ritzi + ncv
  592. iq = bounds + ncv
  593. iw = iq + ldq*ncv
  594. next = iw + ncv**2 + 3*ncv
  595. c
  596. ipntr(4) = next
  597. ipntr(5) = ih
  598. ipntr(6) = ritzr
  599. ipntr(7) = ritzi
  600. ipntr(8) = bounds
  601. ipntr(14) = iw
  602. c
  603. end if
  604. c
  605. c %-------------------------------------------------------%
  606. c | Carry out the Implicitly restarted Arnoldi Iteration. |
  607. c %-------------------------------------------------------%
  608. c
  609. call dnaup2
  610. & ( ido, bmat, n, which, nev0, np, tol, resid, mode, iupd,
  611. & ishift, mxiter, v, ldv, workl(ih), ldh, workl(ritzr),
  612. & workl(ritzi), workl(bounds), workl(iq), ldq, workl(iw),
  613. & ipntr, workd, ITRAK, info )
  614. c recup
  615. nopx =ITRAK(1)
  616. nbx =ITRAK(2)
  617. nrorth=ITRAK(3)
  618. nitref=ITRAK(4)
  619. nrstrt=ITRAK(5)
  620. c
  621. c %--------------------------------------------------%
  622. c | ido .ne. 99 implies use of reverse communication |
  623. c | to compute operations involving OP or shifts. |
  624. c %--------------------------------------------------%
  625. c
  626. if (ido .eq. 3) iparam(8) = np
  627. if (ido .ne. 99) go to 9000
  628. c
  629. iparam(3) = mxiter
  630. iparam(5) = np
  631. iparam(9) = nopx
  632. iparam(10) = nbx
  633. iparam(11) = nrorth
  634. c
  635. c %------------------------------------%
  636. c | Exit if there was an informational |
  637. c | error within dnaup2 . |
  638. c %------------------------------------%
  639. c
  640. if (info .lt. 0) go to 9000
  641. if (info .eq. 2) info = 3
  642. c
  643. if (msglvl .gt. 0) then
  644. call ivout ( 1, mxiter, ndigit,
  645. & '_naupd: Number of update iterations taken')
  646. call ivout ( 1, np, ndigit,
  647. & '_naupd: Number of wanted "converged" Ritz values')
  648. call dvout ( np, workl(ritzr), ndigit,
  649. & '_naupd: Real part of the final Ritz values')
  650. call dvout ( np, workl(ritzi), ndigit,
  651. & '_naupd: Imaginary part of the final Ritz values')
  652. call dvout ( np, workl(bounds), ndigit,
  653. & '_naupd: Associated Ritz estimates')
  654. end if
  655. c
  656. * call arscnd (t1)
  657. c tnaupd = t1 - t0
  658. c
  659. if (msglvl .gt. 0) then
  660. c
  661. c %--------------------------------------------------------%
  662. c | Version Number & Version Date are defined in version.h |
  663. c %--------------------------------------------------------%
  664. c
  665. write (6,1000)
  666. write (6,1100) mxiter, nopx, nbx, nrorth, nitref, nrstrt
  667. c & ,tmvopx, tmvbx, tnaupd, tnaup2, tnaitr, titref,
  668. c & tgetv0, tneigh, tngets, tnapps, tnconv, trvec
  669. 1000 format (//,
  670. & 5x, '=============================================',/
  671. & 5x, '= Nonsymmetric implicit Arnoldi update code =',/
  672. & 5x, '= Version Number: ', ' 2.4' , 21x, ' =',/
  673. & 5x, '= Version Date: ', ' 07/31/96' , 16x, ' =',/
  674. & 5x, '=============================================',//)
  675. c & 5x, '= Summary of timing statistics =',/
  676. c & 5x, '=============================================',//)
  677. 1100 format (
  678. & 5x, 'Total number update iterations = ', i5,/
  679. & 5x, 'Total number of OP*x operations = ', i5,/
  680. & 5x, 'Total number of B*x operations = ', i5,/
  681. & 5x, 'Total number of reorthogonalization steps = ', i5,/
  682. & 5x, 'Total number of iterative refinement steps = ', i5,/
  683. & 5x, 'Total number of restart steps = ', i5,/)
  684. c & 5x, 'Total time in user OP*x operation = ', f12.6,/
  685. c & 5x, 'Total time in user B*x operation = ', f12.6,/
  686. c & 5x, 'Total time in Arnoldi update routine = ', f12.6,/
  687. c & 5x, 'Total time in naup2 routine = ', f12.6,/
  688. c & 5x, 'Total time in basic Arnoldi iteration loop = ', f12.6,/
  689. c & 5x, 'Total time in reorthogonalization phase = ', f12.6,/
  690. c & 5x, 'Total time in (re)start vector generation = ', f12.6,/
  691. c & 5x, 'Total time in Hessenberg eig. subproblem = ', f12.6,/
  692. c & 5x, 'Total time in getting the shifts = ', f12.6,/
  693. c & 5x, 'Total time in applying the shifts = ', f12.6,/
  694. c & 5x, 'Total time in convergence testing = ', f12.6,/
  695. c & 5x, 'Total time in computing final Ritz vectors = ', f12.6/)
  696. end if
  697. c
  698. 9000 continue
  699. c
  700. c ITRAK(..)=... inutile car ITRAK fourni a dsaup2 par ex et pas modifie ici
  701. c
  702. c
  703. c %---------------%
  704. c | End of dnaupd |
  705. c %---------------%
  706. c
  707. end
  708.  
  709.  
  710.  
  711.  

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