Télécharger dsaitr.eso

Retour à la liste

Numérotation des lignes :

  1. C DSAITR SOURCE BP208322 20/02/06 21:15:26 10512
  2. c-----------------------------------------------------------------------
  3. c\BeginDoc
  4. c
  5. c\Name: dsaitr
  6. c
  7. c\Description:
  8. c Reverse communication interface for applying NP additional steps to
  9. c a K step symmetric Arnoldi factorization.
  10. c
  11. c Input: OP*V_{k} - V_{k}*H = r_{k}*e_{k}^T
  12. c
  13. c with (V_{k}^T)*B*V_{k} = I, (V_{k}^T)*B*r_{k} = 0.
  14. c
  15. c Output: OP*V_{k+p} - V_{k+p}*H = r_{k+p}*e_{k+p}^T
  16. c
  17. c with (V_{k+p}^T)*B*V_{k+p} = I, (V_{k+p}^T)*B*r_{k+p} = 0.
  18. c
  19. c where OP and B are as in dsaupd. The B-norm of r_{k+p} is also
  20. c computed and returned.
  21. c
  22. c\Usage:
  23. c call dsaitr
  24. c ( IDO, BMAT, N, K, NP, MODE, RESID, RNORM, V, LDV, H, LDH,
  25. c IPNTR, WORKD, INFO )
  26. c
  27. c\Arguments
  28. c IDO Integer. (INPUT/OUTPUT)
  29. c Reverse communication flag.
  30. c -------------------------------------------------------------
  31. c IDO = 0: first call to the reverse communication interface
  32. c IDO = -1: compute Y = OP * X where
  33. c IPNTR(1) is the pointer into WORK for X,
  34. c IPNTR(2) is the pointer into WORK for Y.
  35. c This is for the restart phase to force the new
  36. c starting vector into the range of OP.
  37. c IDO = 1: compute Y = OP * X where
  38. c IPNTR(1) is the pointer into WORK for X,
  39. c IPNTR(2) is the pointer into WORK for Y,
  40. c IPNTR(3) is the pointer into WORK for B * X.
  41. c IDO = 2: compute Y = B * X where
  42. c IPNTR(1) is the pointer into WORK for X,
  43. c IPNTR(2) is the pointer into WORK for Y.
  44. c IDO = 99: done
  45. c -------------------------------------------------------------
  46. c When the routine is used in the "shift-and-invert" mode, the
  47. c vector B * Q is already available and does not need to be
  48. c recomputed in forming OP * Q.
  49. c
  50. c BMAT Character*1. (INPUT)
  51. c BMAT specifies the type of matrix B that defines the
  52. c semi-inner product for the operator OP. See dsaupd.
  53. c B = 'I' -> standard eigenvalue problem A*x = lambda*x
  54. c B = 'G' -> generalized eigenvalue problem A*x = lambda*M*x
  55. c
  56. c N Integer. (INPUT)
  57. c Dimension of the eigenproblem.
  58. c
  59. c K Integer. (INPUT)
  60. c Current order of H and the number of columns of V.
  61. c
  62. c NP Integer. (INPUT)
  63. c Number of additional Arnoldi steps to take.
  64. c
  65. c MODE Integer. (INPUT)
  66. c Signifies which form for "OP". If MODE=2 then
  67. c a reduction in the number of B matrix vector multiplies
  68. c is possible since the B-norm of OP*x is equivalent to
  69. c the inv(B)-norm of A*x.
  70. c
  71. c RESID Double precision array of length N. (INPUT/OUTPUT)
  72. c On INPUT: RESID contains the residual vector r_{k}.
  73. c On OUTPUT: RESID contains the residual vector r_{k+p}.
  74. c
  75. c RNORM Double precision scalar. (INPUT/OUTPUT)
  76. c On INPUT the B-norm of r_{k}.
  77. c On OUTPUT the B-norm of the updated residual r_{k+p}.
  78. c
  79. c V REAL*8 N by K+NP array. (INPUT/OUTPUT)
  80. c On INPUT: V contains the Arnoldi vectors in the first K
  81. c columns.
  82. c On OUTPUT: V contains the new NP Arnoldi vectors in the next
  83. c NP columns. The first K columns are unchanged.
  84. c
  85. c LDV Integer. (INPUT)
  86. c Leading dimension of V exactly as declared in the calling
  87. c program.
  88. c
  89. c H REAL*8 (K+NP) by 2 array. (INPUT/OUTPUT)
  90. c H is used to store the generated symmetric tridiagonal matrix
  91. c with the subdiagonal in the first column starting at H(2,1)
  92. c and the main diagonal in the second column.
  93. c
  94. c LDH Integer. (INPUT)
  95. c Leading dimension of H exactly as declared in the calling
  96. c program.
  97. c
  98. c IPNTR Integer array of length 3. (OUTPUT)
  99. c Pointer to mark the starting locations in the WORK for
  100. c vectors used by the Arnoldi iteration.
  101. c -------------------------------------------------------------
  102. c IPNTR(1): pointer to the current operand vector X.
  103. c IPNTR(2): pointer to the current result vector Y.
  104. c IPNTR(3): pointer to the vector B * X when used in the
  105. c shift-and-invert mode. X is the current operand.
  106. c -------------------------------------------------------------
  107. c
  108. c WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION)
  109. c Distributed array to be used in the basic Arnoldi iteration
  110. c for reverse communication. The calling program should not
  111. c use WORKD as temporary workspace during the iteration !!!!!!
  112. c On INPUT, WORKD(1:N) = B*RESID where RESID is associated
  113. c with the K step Arnoldi factorization. Used to save some
  114. c computation at the first step.
  115. c On OUTPUT, WORKD(1:N) = B*RESID where RESID is associated
  116. c with the K+NP step Arnoldi factorization.
  117. c
  118. c INFO Integer. (OUTPUT)
  119. c = 0: Normal exit.
  120. c > 0: Size of an invariant subspace of OP is found that is
  121. c less than K + NP.
  122. c
  123. c\EndDoc
  124. c
  125. c-----------------------------------------------------------------------
  126. c
  127. c\BeginLib
  128. c
  129. c\Local variables:
  130. c xxxxxx real
  131. c
  132. c\Routines called:
  133. c dgetv0 ARPACK routine to generate the initial vector.
  134. c ivout ARPACK utility routine that prints integers.
  135. c dmout ARPACK utility routine that prints matrices.
  136. c dvout ARPACK utility routine that prints vectors.
  137. c dlamch LAPACK routine that determines machine constants.
  138. c dlascl LAPACK routine for careful scaling of a matrix.
  139. c dgemv Level 2 BLAS routine for matrix vector multiplication.
  140. c daxpy Level 1 BLAS that computes a vector triad.
  141. c dscal Level 1 BLAS that scales a vector.
  142. c dcopy Level 1 BLAS that copies one vector to another .
  143. c ddot Level 1 BLAS that computes the scalar product of two vectors.
  144. c dnrm2 Level 1 BLAS that computes the norm of a vector.
  145. c
  146. c\Author
  147. c Danny Sorensen Phuong Vu
  148. c Richard Lehoucq CRPC / Rice University
  149. c Dept. of Computational & Houston, Texas
  150. c Applied Mathematics
  151. c Rice University
  152. c Houston, Texas
  153. c
  154. c\Revision history:
  155. c xx/xx/93: Version ' 2.4'
  156. c
  157. c\SCCS Information: @(#)
  158. c FILE: saitr.F SID: 2.6 DATE OF SID: 8/28/96 RELEASE: 2
  159. c
  160. c\Remarks
  161. c The algorithm implemented is:
  162. c
  163. c restart = .false.
  164. c Given V_{k} = [v_{1}, ..., v_{k}], r_{k};
  165. c r_{k} contains the initial residual vector even for k = 0;
  166. c Also assume that rnorm = || B*r_{k} || and B*r_{k} are already
  167. c computed by the calling program.
  168. c
  169. c betaj = rnorm ; p_{k+1} = B*r_{k} ;
  170. c For j = k+1, ..., k+np Do
  171. c 1) if ( betaj < tol ) stop or restart depending on j.
  172. c if ( restart ) generate a new starting vector.
  173. c 2) v_{j} = r(j-1)/betaj; V_{j} = [V_{j-1}, v_{j}];
  174. c p_{j} = p_{j}/betaj
  175. c 3) r_{j} = OP*v_{j} where OP is defined as in dsaupd
  176. c For shift-invert mode p_{j} = B*v_{j} is already available.
  177. c wnorm = || OP*v_{j} ||
  178. c 4) Compute the j-th step residual vector.
  179. c w_{j} = V_{j}^T * B * OP * v_{j}
  180. c r_{j} = OP*v_{j} - V_{j} * w_{j}
  181. c alphaj <- j-th component of w_{j}
  182. c rnorm = || r_{j} ||
  183. c betaj+1 = rnorm
  184. c If (rnorm > 0.717*wnorm) accept step and go back to 1)
  185. c 5) Re-orthogonalization step:
  186. c s = V_{j}'*B*r_{j}
  187. c r_{j} = r_{j} - V_{j}*s; rnorm1 = || r_{j} ||
  188. c alphaj = alphaj + s_{j};
  189. c 6) Iterative refinement step:
  190. c If (rnorm1 > 0.717*rnorm) then
  191. c rnorm = rnorm1
  192. c accept step and go back to 1)
  193. c Else
  194. c rnorm = rnorm1
  195. c If this is the first time in step 6), go to 5)
  196. c Else r_{j} lies in the span of V_{j} numerically.
  197. c Set r_{j} = 0 and rnorm = 0; go to 1)
  198. c EndIf
  199. c End Do
  200. c
  201. c\EndLib
  202. c
  203. c-----------------------------------------------------------------------
  204. c
  205. subroutine dsaitr
  206. & (ido, bmat, n, k, np, mode, resid, rnorm, v, ldv, h, ldh,
  207. & ipntr, workd, ITRAK, info)
  208. c
  209. c %----------------------------------------------------%
  210. c | Include files for debugging and timing information |
  211. c -INC TARTRAK
  212. c %----------------------------------------------------%
  213. c
  214. c
  215. c %------------------%
  216. c | Scalar Arguments |
  217. c %------------------%
  218. c
  219. character bmat*1
  220. integer ido, info, k, ldh, ldv, n, mode, np
  221. REAL*8
  222. & rnorm
  223. real*8 T0,T1,T2,T3,T4,T5
  224. c
  225. c %-----------------%
  226. c | Array Arguments |
  227. c %-----------------%
  228. c
  229. integer ipntr(3)
  230. INTEGER ITRAK(5)
  231. REAL*8
  232. & h(ldh,2), resid(n), v(ldv,k+np), workd(3*n)
  233. c
  234. c %------------%
  235. c | Parameters |
  236. c %------------%
  237. c
  238. REAL*8
  239. & one, zero
  240. parameter (one = 1.0D+0, zero = 0.0D+0)
  241. c
  242. c %---------------%
  243. c | Local Scalars |
  244. c %---------------%
  245. c
  246. logical first, orth1, orth2, rstart, step3, step4
  247. integer i, ierr, ipj, irj, ivj, iter, itry, j, msglvl,
  248. & infol, jj
  249. REAL*8
  250. & rnorm1, wnorm, safmin, temp1
  251. save orth1, orth2, rstart, step3, step4,
  252. c & ierr, ipj, irj, ivj, iter, itry, j, msglvl,
  253. & ierr, ipj, irj, ivj, iter, itry, j,
  254. & rnorm1, safmin, wnorm
  255. parameter (msglvl=0)
  256. c
  257. c %-----------------------%
  258. c | Local Array Arguments |
  259. c %-----------------------%
  260. c
  261. REAL*8
  262. & xtemp(2)
  263. c
  264. c %----------------------%
  265. c | External Subroutines |
  266. c %----------------------%
  267. c
  268. c
  269. c %--------------------%
  270. c | External Functions |
  271. c %--------------------%
  272. c
  273. REAL*8
  274. external ddot, dnrm2, dlamch
  275. c
  276. c %-----------------%
  277. c | Data statements |
  278. c %-----------------%
  279. c
  280. data first / .true. /
  281. c
  282. c %-----------------------%
  283. c | Executable Statements |
  284. c %-----------------------%
  285. c T0=0.D0
  286. c T1=0.D0
  287. c T2=0.D0
  288. c T3=0.D0
  289. c T4=0.D0
  290. c T5=0.D0
  291. NOPX =ITRAK(1)
  292. NBX =ITRAK(2)
  293. NRORTH=ITRAK(3)
  294. NITREF=ITRAK(4)
  295. NRSTRT=ITRAK(5)
  296. c
  297. if (first) then
  298. first = .false.
  299. c
  300. c %--------------------------------%
  301. c | safmin = safe minimum is such |
  302. c | that 1/sfmin does not overflow |
  303. c %--------------------------------%
  304. c
  305. safmin = dlamch('safmin')
  306. end if
  307. c
  308. if (ido .eq. 0) then
  309. c
  310. c %-------------------------------%
  311. c | Initialize timing statistics |
  312. c | & message level for debugging |
  313. c %-------------------------------%
  314. c
  315. * call arscnd (t0)
  316. c msglvl = msaitr
  317. c
  318. c %------------------------------%
  319. c | Initial call to this routine |
  320. c %------------------------------%
  321. c
  322. info = 0
  323. step3 = .false.
  324. step4 = .false.
  325. rstart = .false.
  326. orth1 = .false.
  327. orth2 = .false.
  328. c
  329. c %--------------------------------%
  330. c | Pointer to the current step of |
  331. c | the factorization to build |
  332. c %--------------------------------%
  333. c
  334. j = k + 1
  335. c
  336. c %------------------------------------------%
  337. c | Pointers used for reverse communication |
  338. c | when using WORKD. |
  339. c %------------------------------------------%
  340. c
  341. ipj = 1
  342. irj = ipj + n
  343. ivj = irj + n
  344. end if
  345. c
  346. c %-------------------------------------------------%
  347. c | When in reverse communication mode one of: |
  348. c | STEP3, STEP4, ORTH1, ORTH2, RSTART |
  349. c | will be .true. |
  350. c | STEP3: return from computing OP*v_{j}. |
  351. c | STEP4: return from computing B-norm of OP*v_{j} |
  352. c | ORTH1: return from computing B-norm of r_{j+1} |
  353. c | ORTH2: return from computing B-norm of |
  354. c | correction to the residual vector. |
  355. c | RSTART: return from OP computations needed by |
  356. c | dgetv0. |
  357. c %-------------------------------------------------%
  358. c
  359. if (step3) go to 50
  360. if (step4) go to 60
  361. if (orth1) go to 70
  362. if (orth2) go to 90
  363. if (rstart) go to 30
  364. c
  365. c %------------------------------%
  366. c | Else this is the first step. |
  367. c %------------------------------%
  368. c
  369. c %--------------------------------------------------------------%
  370. c | |
  371. c | A R N O L D I I T E R A T I O N L O O P |
  372. c | |
  373. c | Note: B*r_{j-1} is already in WORKD(1:N)=WORKD(IPJ:IPJ+N-1) |
  374. c %--------------------------------------------------------------%
  375. c
  376. 1000 continue
  377. c
  378. c if (msglvl .gt. 2) then
  379. c call ivout ( 1, j, ndigit,
  380. c & '_saitr: generating Arnoldi vector no.')
  381. c call dvout ( 1, rnorm, ndigit,
  382. c & '_saitr: B-norm of the current residual =')
  383. c end if
  384. c
  385. c %---------------------------------------------------------%
  386. c | Check for exact zero. Equivalent to determing whether a |
  387. c | j-step Arnoldi factorization is present. |
  388. c %---------------------------------------------------------%
  389. c
  390. if (rnorm .gt. zero) go to 40
  391. c
  392. c %---------------------------------------------------%
  393. c | Invariant subspace found, generate a new starting |
  394. c | vector which is orthogonal to the current Arnoldi |
  395. c | basis and continue the iteration. |
  396. c %---------------------------------------------------%
  397. c
  398. if (msglvl .gt. 0) then
  399. call ivout ( 1, j, ndigit,
  400. & '_saitr: ****** restart at step ******')
  401. end if
  402. c
  403. c %---------------------------------------------%
  404. c | ITRY is the loop variable that controls the |
  405. c | maximum amount of times that a restart is |
  406. c %---------------------------------------------%
  407. c
  408. nrstrt = nrstrt + 1
  409. itry = 1
  410. 20 continue
  411. rstart = .true.
  412. ido = 0
  413. 30 continue
  414. c
  415. c %--------------------------------------%
  416. c | If in reverse communication mode and |
  417. c | RSTART = .true. flow returns here. |
  418. c %--------------------------------------%
  419. c
  420. call dgetv0 (ido, bmat, itry, .false., n, j, v, ldv,
  421. & resid, rnorm, ipntr, workd, nbx,nopx, ierr)
  422. if (ido .ne. 99) go to 9000
  423. if (ierr .lt. 0) then
  424. itry = itry + 1
  425. if (itry .le. 3) go to 20
  426. c
  427. c %------------------------------------------------%
  428. c | Give up after several restart attempts. |
  429. c | Set INFO to the size of the invariant subspace |
  430. c | which spans OP and exit. |
  431. c %------------------------------------------------%
  432. c
  433. info = j - 1
  434. * call arscnd (t1)
  435. c tsaitr = tsaitr + (t1 - t0)
  436. ido = 99
  437. go to 9000
  438. end if
  439. c
  440. 40 continue
  441. c
  442. c %---------------------------------------------------------%
  443. c | STEP 2: v_{j} = r_{j-1}/rnorm and p_{j} = p_{j}/rnorm |
  444. c | Note that p_{j} = B*r_{j-1}. In order to avoid overflow |
  445. c | when reciprocating a small RNORM, test against lower |
  446. c | machine bound. |
  447. c %---------------------------------------------------------%
  448. c
  449. call dcopy (n, resid, 1, v(1,j), 1)
  450. if (rnorm .ge. safmin) then
  451. temp1 = one / rnorm
  452. call dscal (n, temp1, v(1,j), 1)
  453. call dscal (n, temp1, workd(ipj), 1)
  454. else
  455. c
  456. c %-----------------------------------------%
  457. c | To scale both v_{j} and p_{j} carefully |
  458. c | use LAPACK routine SLASCL |
  459. c %-----------------------------------------%
  460. c
  461. call dlascl ('General', i, i, rnorm, one, n, 1,
  462. & v(1,j), n, infol)
  463. call dlascl ('General', i, i, rnorm, one, n, 1,
  464. & workd(ipj), n, infol)
  465. end if
  466. c
  467. c %------------------------------------------------------%
  468. c | STEP 3: r_{j} = OP*v_{j}; Note that p_{j} = B*v_{j} |
  469. c | Note that this is not quite yet r_{j}. See STEP 4 |
  470. c %------------------------------------------------------%
  471. c
  472. step3 = .true.
  473. nopx = nopx + 1
  474. * call arscnd (t2)
  475. call dcopy (n, v(1,j), 1, workd(ivj), 1)
  476. ipntr(1) = ivj
  477. ipntr(2) = irj
  478. ipntr(3) = ipj
  479. ido = 1
  480. c
  481. c %-----------------------------------%
  482. c | Exit in order to compute OP*v_{j} |
  483. c %-----------------------------------%
  484. c
  485. go to 9000
  486. 50 continue
  487. c
  488. c %-----------------------------------%
  489. c | Back from reverse communication; |
  490. c | WORKD(IRJ:IRJ+N-1) := OP*v_{j}. |
  491. c %-----------------------------------%
  492. c
  493. * call arscnd (t3)
  494. c tmvopx = tmvopx + (t3 - t2)
  495. c
  496. step3 = .false.
  497. c
  498. c %------------------------------------------%
  499. c | Put another copy of OP*v_{j} into RESID. |
  500. c %------------------------------------------%
  501. c
  502. call dcopy (n, workd(irj), 1, resid, 1)
  503. c
  504. c %-------------------------------------------%
  505. c | STEP 4: Finish extending the symmetric |
  506. c | Arnoldi to length j. If MODE = 2 |
  507. c | then B*OP = B*inv(B)*A = A and |
  508. c | we don't need to compute B*OP. |
  509. c | NOTE: If MODE = 2 WORKD(IVJ:IVJ+N-1) is |
  510. c | assumed to have A*v_{j}. |
  511. c %-------------------------------------------%
  512. c
  513. if (mode .eq. 2) go to 65
  514. * call arscnd (t2)
  515. if (bmat .eq. 'G') then
  516. nbx = nbx + 1
  517. step4 = .true.
  518. ipntr(1) = irj
  519. ipntr(2) = ipj
  520. ido = 2
  521. c
  522. c %-------------------------------------%
  523. c | Exit in order to compute B*OP*v_{j} |
  524. c %-------------------------------------%
  525. c
  526. go to 9000
  527. else if (bmat .eq. 'I') then
  528. call dcopy(n, resid, 1 , workd(ipj), 1)
  529. end if
  530. 60 continue
  531. c
  532. c %-----------------------------------%
  533. c | Back from reverse communication; |
  534. c | WORKD(IPJ:IPJ+N-1) := B*OP*v_{j}. |
  535. c %-----------------------------------%
  536. c
  537. c if (bmat .eq. 'G') then
  538. * call arscnd (t3)
  539. c tmvbx = tmvbx + (t3 - t2)
  540. c end if
  541. c
  542. step4 = .false.
  543. c
  544. c %-------------------------------------%
  545. c | The following is needed for STEP 5. |
  546. c | Compute the B-norm of OP*v_{j}. |
  547. c %-------------------------------------%
  548. c
  549. 65 continue
  550. if (mode .eq. 2) then
  551. c
  552. c %----------------------------------%
  553. c | Note that the B-norm of OP*v_{j} |
  554. c | is the inv(B)-norm of A*v_{j}. |
  555. c %----------------------------------%
  556. c
  557. wnorm = ddot (n, resid, 1, workd(ivj), 1)
  558. wnorm = sqrt(abs(wnorm))
  559. else if (bmat .eq. 'G') then
  560. wnorm = ddot (n, resid, 1, workd(ipj), 1)
  561. wnorm = sqrt(abs(wnorm))
  562. else if (bmat .eq. 'I') then
  563. wnorm = dnrm2(n, resid, 1)
  564. end if
  565. c
  566. c %-----------------------------------------%
  567. c | Compute the j-th residual corresponding |
  568. c | to the j step factorization. |
  569. c | Use Classical Gram Schmidt and compute: |
  570. c | w_{j} <- V_{j}^T * B * OP * v_{j} |
  571. c | r_{j} <- OP*v_{j} - V_{j} * w_{j} |
  572. c %-----------------------------------------%
  573. c
  574. c
  575. c %------------------------------------------%
  576. c | Compute the j Fourier coefficients w_{j} |
  577. c | WORKD(IPJ:IPJ+N-1) contains B*OP*v_{j}. |
  578. c %------------------------------------------%
  579. c
  580. if (mode .ne. 2 ) then
  581. call dgemv('T', n, j, one, v, ldv, workd(ipj), 1, zero,
  582. & workd(irj), 1)
  583. else if (mode .eq. 2) then
  584. call dgemv('T', n, j, one, v, ldv, workd(ivj), 1, zero,
  585. & workd(irj), 1)
  586. end if
  587. c
  588. c %--------------------------------------%
  589. c | Orthgonalize r_{j} against V_{j}. |
  590. c | RESID contains OP*v_{j}. See STEP 3. |
  591. c %--------------------------------------%
  592. c
  593. call dgemv('N', n, j, -one, v, ldv, workd(irj), 1, one,
  594. & resid, 1)
  595. c
  596. c %--------------------------------------%
  597. c | Extend H to have j rows and columns. |
  598. c %--------------------------------------%
  599. c
  600. h(j,2) = workd(irj + j - 1)
  601. if (j .eq. 1 .or. rstart) then
  602. h(j,1) = zero
  603. else
  604. h(j,1) = rnorm
  605. end if
  606. * call arscnd (t4)
  607. c
  608. orth1 = .true.
  609. iter = 0
  610. c
  611. * call arscnd (t2)
  612. if (bmat .eq. 'G') then
  613. nbx = nbx + 1
  614. call dcopy (n, resid, 1, workd(irj), 1)
  615. ipntr(1) = irj
  616. ipntr(2) = ipj
  617. ido = 2
  618. c
  619. c %----------------------------------%
  620. c | Exit in order to compute B*r_{j} |
  621. c %----------------------------------%
  622. c
  623. go to 9000
  624. else if (bmat .eq. 'I') then
  625. call dcopy (n, resid, 1, workd(ipj), 1)
  626. end if
  627. 70 continue
  628. c
  629. c %---------------------------------------------------%
  630. c | Back from reverse communication if ORTH1 = .true. |
  631. c | WORKD(IPJ:IPJ+N-1) := B*r_{j}. |
  632. c %---------------------------------------------------%
  633. c
  634. c if (bmat .eq. 'G') then
  635. * call arscnd (t3)
  636. c tmvbx = tmvbx + (t3 - t2)
  637. c end if
  638. c
  639. orth1 = .false.
  640. c
  641. c %------------------------------%
  642. c | Compute the B-norm of r_{j}. |
  643. c %------------------------------%
  644. c
  645. if (bmat .eq. 'G') then
  646. rnorm = ddot (n, resid, 1, workd(ipj), 1)
  647. rnorm = sqrt(abs(rnorm))
  648. else if (bmat .eq. 'I') then
  649. rnorm = dnrm2(n, resid, 1)
  650. end if
  651. c
  652. c %-----------------------------------------------------------%
  653. c | STEP 5: Re-orthogonalization / Iterative refinement phase |
  654. c | Maximum NITER_ITREF tries. |
  655. c | |
  656. c | s = V_{j}^T * B * r_{j} |
  657. c | r_{j} = r_{j} - V_{j}*s |
  658. c | alphaj = alphaj + s_{j} |
  659. c | |
  660. c | The stopping criteria used for iterative refinement is |
  661. c | discussed in Parlett's book SEP, page 107 and in Gragg & |
  662. c | Reichel ACM TOMS paper; Algorithm 686, Dec. 1990. |
  663. c | Determine if we need to correct the residual. The goal is |
  664. c | to enforce ||v(:,1:j)^T * r_{j}|| .le. eps * || r_{j} || |
  665. c %-----------------------------------------------------------%
  666. c
  667. if (rnorm .gt. 0.717*wnorm) go to 100
  668. nrorth = nrorth + 1
  669. c
  670. c %---------------------------------------------------%
  671. c | Enter the Iterative refinement phase. If further |
  672. c | refinement is necessary, loop back here. The loop |
  673. c | variable is ITER. Perform a step of Classical |
  674. c | Gram-Schmidt using all the Arnoldi vectors V_{j} |
  675. c %---------------------------------------------------%
  676. c
  677. 80 continue
  678. c
  679. c if (msglvl .gt. 2) then
  680. c xtemp(1) = wnorm
  681. c xtemp(2) = rnorm
  682. c call dvout ( 2, xtemp, ndigit,
  683. c & '_saitr: re-orthonalization ; wnorm and rnorm are')
  684. c end if
  685. c
  686. c %----------------------------------------------------%
  687. c | Compute V_{j}^T * B * r_{j}. |
  688. c | WORKD(IRJ:IRJ+J-1) = v(:,1:J)'*WORKD(IPJ:IPJ+N-1). |
  689. c %----------------------------------------------------%
  690. c
  691. call dgemv ('T', n, j, one, v, ldv, workd(ipj), 1,
  692. & zero, workd(irj), 1)
  693. c
  694. c %----------------------------------------------%
  695. c | Compute the correction to the residual: |
  696. c | r_{j} = r_{j} - V_{j} * WORKD(IRJ:IRJ+J-1). |
  697. c | The correction to H is v(:,1:J)*H(1:J,1:J) + |
  698. c | v(:,1:J)*WORKD(IRJ:IRJ+J-1)*e'_j, but only |
  699. c | H(j,j) is updated. |
  700. c %----------------------------------------------%
  701. c
  702. call dgemv ('N', n, j, -one, v, ldv, workd(irj), 1,
  703. & one, resid, 1)
  704. c
  705. if (j .eq. 1 .or. rstart) h(j,1) = zero
  706. h(j,2) = h(j,2) + workd(irj + j - 1)
  707. c
  708. orth2 = .true.
  709. * call arscnd (t2)
  710. if (bmat .eq. 'G') then
  711. nbx = nbx + 1
  712. call dcopy (n, resid, 1, workd(irj), 1)
  713. ipntr(1) = irj
  714. ipntr(2) = ipj
  715. ido = 2
  716. c
  717. c %-----------------------------------%
  718. c | Exit in order to compute B*r_{j}. |
  719. c | r_{j} is the corrected residual. |
  720. c %-----------------------------------%
  721. c
  722. go to 9000
  723. else if (bmat .eq. 'I') then
  724. call dcopy (n, resid, 1, workd(ipj), 1)
  725. end if
  726. 90 continue
  727. c
  728. c %---------------------------------------------------%
  729. c | Back from reverse communication if ORTH2 = .true. |
  730. c %---------------------------------------------------%
  731. c
  732. c if (bmat .eq. 'G') then
  733. * call arscnd (t3)
  734. c tmvbx = tmvbx + (t3 - t2)
  735. c end if
  736. c
  737. c %-----------------------------------------------------%
  738. c | Compute the B-norm of the corrected residual r_{j}. |
  739. c %-----------------------------------------------------%
  740. c
  741. if (bmat .eq. 'G') then
  742. rnorm1 = ddot (n, resid, 1, workd(ipj), 1)
  743. rnorm1 = sqrt(abs(rnorm1))
  744. else if (bmat .eq. 'I') then
  745. rnorm1 = dnrm2(n, resid, 1)
  746. end if
  747. c
  748. if (msglvl .gt. 0 .and. iter .gt. 0) then
  749. call ivout ( 1, j, ndigit,
  750. & '_saitr: Iterative refinement for Arnoldi residual')
  751. c if (msglvl .gt. 2) then
  752. c xtemp(1) = rnorm
  753. c xtemp(2) = rnorm1
  754. c call dvout ( 2, xtemp, ndigit,
  755. c & '_saitr: iterative refinement ; rnorm and rnorm1 are')
  756. c end if
  757. end if
  758. c
  759. c %-----------------------------------------%
  760. c | Determine if we need to perform another |
  761. c | step of re-orthogonalization. |
  762. c %-----------------------------------------%
  763. c
  764. if (rnorm1 .gt. 0.717*rnorm) then
  765. c
  766. c %--------------------------------%
  767. c | No need for further refinement |
  768. c %--------------------------------%
  769. c
  770. rnorm = rnorm1
  771. c
  772. else
  773. c
  774. c %-------------------------------------------%
  775. c | Another step of iterative refinement step |
  776. c %-------------------------------------------%
  777. c
  778. nitref = nitref + 1
  779. rnorm = rnorm1
  780. iter = iter + 1
  781. if (iter .le. 1) go to 80
  782. c
  783. c %-------------------------------------------------%
  784. c | Otherwise RESID is numerically in the span of V |
  785. c %-------------------------------------------------%
  786. c
  787. do 95 jj = 1, n
  788. resid(jj) = zero
  789. 95 continue
  790. rnorm = zero
  791. end if
  792. c
  793. c %----------------------------------------------%
  794. c | Branch here directly if iterative refinement |
  795. c | wasn't necessary or after at most NITER_REF |
  796. c | steps of iterative refinement. |
  797. c %----------------------------------------------%
  798. c
  799. 100 continue
  800. c
  801. rstart = .false.
  802. orth2 = .false.
  803. c
  804. * call arscnd (t5)
  805. c titref = titref + (t5 - t4)
  806. c
  807. c %----------------------------------------------------------%
  808. c | Make sure the last off-diagonal element is non negative |
  809. c | If not perform a similarity transformation on H(1:j,1:j) |
  810. c | and scale v(:,j) by -1. |
  811. c %----------------------------------------------------------%
  812. c
  813. if (h(j,1) .lt. zero) then
  814. h(j,1) = -h(j,1)
  815. if ( j .lt. k+np) then
  816. call dscal(n, -one, v(1,j+1), 1)
  817. else
  818. call dscal(n, -one, resid, 1)
  819. end if
  820. end if
  821. c
  822. c %------------------------------------%
  823. c | STEP 6: Update j = j+1; Continue |
  824. c %------------------------------------%
  825. c
  826. j = j + 1
  827. if (j .gt. k+np) then
  828. * call arscnd (t1)
  829. c tsaitr = tsaitr + (t1 - t0)
  830. ido = 99
  831. c
  832. c if (msglvl .gt. 1) then
  833. c call dvout ( k+np, h(1,2), ndigit,
  834. c & '_saitr: main diagonal of matrix H of step K+NP.')
  835. c if (k+np .gt. 1) then
  836. c call dvout ( k+np-1, h(2,1), ndigit,
  837. c & '_saitr: sub diagonal of matrix H of step K+NP.')
  838. c end if
  839. c end if
  840. c
  841. go to 9000
  842. end if
  843. c
  844. c %--------------------------------------------------------%
  845. c | Loop back to extend the factorization by another step. |
  846. c %--------------------------------------------------------%
  847. c
  848. go to 1000
  849. c
  850. c %---------------------------------------------------------------%
  851. c | |
  852. c | E N D O F M A I N I T E R A T I O N L O O P |
  853. c | |
  854. c %---------------------------------------------------------------%
  855. c
  856. 9000 continue
  857. ITRAK(1)=NOPX
  858. ITRAK(2)=NBX
  859. ITRAK(3)=NRORTH
  860. ITRAK(4)=NITREF
  861. ITRAK(5)=NRSTRT
  862. c
  863. c %---------------%
  864. c | End of dsaitr |
  865. c %---------------%
  866. c
  867. end
  868.  
  869.  
  870.  
  871.  

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