Télécharger dnaitr.eso

Retour à la liste

Numérotation des lignes :

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

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