Télécharger dsapps.eso

Retour à la liste

Numérotation des lignes :

  1. C DSAPPS SOURCE BP208322 15/12/17 21:15:10 8750
  2. c-----------------------------------------------------------------------
  3. c\BeginDoc
  4. c
  5. c\Name: dsapps
  6. c
  7. c\Description:
  8. c Given the Arnoldi factorization
  9. c
  10. c A*V_{k} - V_{k}*H_{k} = r_{k+p}*e_{k+p}^T,
  11. c
  12. c apply NP shifts implicitly resulting in
  13. c
  14. c A*(V_{k}*Q) - (V_{k}*Q)*(Q^T* H_{k}*Q) = r_{k+p}*e_{k+p}^T * Q
  15. c
  16. c where Q is an orthogonal matrix of order KEV+NP. Q is the product of
  17. c rotations resulting from the NP bulge chasing sweeps. The updated Arnoldi
  18. c factorization becomes:
  19. c
  20. c A*VNEW_{k} - VNEW_{k}*HNEW_{k} = rnew_{k}*e_{k}^T.
  21. c
  22. c\Usage:
  23. c call dsapps
  24. c ( N, KEV, NP, SHIFT, V, LDV, H, LDH, RESID, Q, LDQ, WORKD )
  25. c
  26. c\Arguments
  27. c N Integer. (INPUT)
  28. c Problem size, i.e. dimension of matrix A.
  29. c
  30. c KEV Integer. (INPUT)
  31. c INPUT: KEV+NP is the size of the input matrix H.
  32. c OUTPUT: KEV is the size of the updated matrix HNEW.
  33. c
  34. c NP Integer. (INPUT)
  35. c Number of implicit shifts to be applied.
  36. c
  37. c SHIFT Double precision array of length NP. (INPUT)
  38. c The shifts to be applied.
  39. c
  40. c V REAL*8 N by (KEV+NP) array. (INPUT/OUTPUT)
  41. c INPUT: V contains the current KEV+NP Arnoldi vectors.
  42. c OUTPUT: VNEW = V(1:n,1:KEV); the updated Arnoldi vectors
  43. c are in the first KEV columns of V.
  44. c
  45. c LDV Integer. (INPUT)
  46. c Leading dimension of V exactly as declared in the calling
  47. c program.
  48. c
  49. c H REAL*8 (KEV+NP) by 2 array. (INPUT/OUTPUT)
  50. c INPUT: H contains the symmetric tridiagonal matrix of the
  51. c Arnoldi factorization with the subdiagonal in the 1st column
  52. c starting at H(2,1) and the main diagonal in the 2nd column.
  53. c OUTPUT: H contains the updated tridiagonal matrix in the
  54. c KEV leading submatrix.
  55. c
  56. c LDH Integer. (INPUT)
  57. c Leading dimension of H exactly as declared in the calling
  58. c program.
  59. c
  60. c RESID Double precision array of length (N). (INPUT/OUTPUT)
  61. c INPUT: RESID contains the the residual vector r_{k+p}.
  62. c OUTPUT: RESID is the updated residual vector rnew_{k}.
  63. c
  64. c Q REAL*8 KEV+NP by KEV+NP work array. (WORKSPACE)
  65. c Work array used to accumulate the rotations during the bulge
  66. c chase sweep.
  67. c
  68. c LDQ Integer. (INPUT)
  69. c Leading dimension of Q exactly as declared in the calling
  70. c program.
  71. c
  72. c WORKD Double precision work array of length 2*N. (WORKSPACE)
  73. c Distributed array used in the application of the accumulated
  74. c orthogonal matrix Q.
  75. c
  76. c\EndDoc
  77. c
  78. c-----------------------------------------------------------------------
  79. c
  80. c\BeginLib
  81. c
  82. c\Local variables:
  83. c xxxxxx real
  84. c
  85. c\References:
  86. c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
  87. c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
  88. c pp 357-385.
  89. c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
  90. c Restarted Arnoldi Iteration", Rice University Technical Report
  91. c TR95-13, Department of Computational and Applied Mathematics.
  92. c
  93. c\Routines called:
  94. c ivout ARPACK utility routine that prints integers.
  95. c arscnd ARPACK utility routine for timing.
  96. c dvout ARPACK utility routine that prints vectors.
  97. c dlamch LAPACK routine that determines machine constants.
  98. c dlartg LAPACK Givens rotation construction routine.
  99. c dlacpy LAPACK matrix copy routine.
  100. c dlaset LAPACK matrix initialization routine.
  101. c dgemv Level 2 BLAS routine for matrix vector multiplication.
  102. c daxpy Level 1 BLAS that computes a vector triad.
  103. c dcopy Level 1 BLAS that copies one vector to another.
  104. c dscal Level 1 BLAS that scales a vector.
  105. c
  106. c\Author
  107. c Danny Sorensen Phuong Vu
  108. c Richard Lehoucq CRPC / Rice University
  109. c Dept. of Computational & Houston, Texas
  110. c Applied Mathematics
  111. c Rice University
  112. c Houston, Texas
  113. c
  114. c\Revision history:
  115. c 12/16/93: Version ' 2.4'
  116. c
  117. c\SCCS Information: @(#)
  118. c FILE: sapps.F SID: 2.6 DATE OF SID: 3/28/97 RELEASE: 2
  119. c
  120. c\Remarks
  121. c 1. In this version, each shift is applied to all the subblocks of
  122. c the tridiagonal matrix H and not just to the submatrix that it
  123. c comes from. This routine assumes that the subdiagonal elements
  124. c of H that are stored in h(1:kev+np,1) are nonegative upon input
  125. c and enforce this condition upon output. This version incorporates
  126. c deflation. See code for documentation.
  127. c
  128. c\EndLib
  129. c
  130. c-----------------------------------------------------------------------
  131. c
  132. subroutine dsapps
  133. & ( n, kev, np, shift, v, ldv, h, ldh, resid, q, ldq, workd )
  134. c
  135. c %----------------------------------------------------%
  136. c | Include files for debugging and timing information |
  137. -INC TARTRAK
  138. c %----------------------------------------------------%
  139. c
  140. c
  141. c %------------------%
  142. c | Scalar Arguments |
  143. c %------------------%
  144. c
  145. integer kev, ldh, ldq, ldv, n, np
  146. c
  147. c %-----------------%
  148. c | Array Arguments |
  149. c %-----------------%
  150. c
  151. REAL*8
  152. & h(ldh,2), q(ldq,kev+np), resid(n), shift(np),
  153. & v(ldv,kev+np), workd(2*n)
  154. c
  155. c %------------%
  156. c | Parameters |
  157. c %------------%
  158. c
  159. REAL*8
  160. & one, zero
  161. parameter (one = 1.0D+0, zero = 0.0D+0)
  162. c
  163. c %---------------%
  164. c | Local Scalars |
  165. c %---------------%
  166. c
  167. integer i, iend, istart, itop, j, jj, kplusp, msglvl
  168. logical first
  169. REAL*8
  170. & a1, a2, a3, a4, big, c, epsmch, f, g, r, s
  171. save epsmch, first
  172. c
  173. c
  174. c %----------------------%
  175. c | External Subroutines |
  176. c %----------------------%
  177. c
  178. & ivout, arscnd, dgemv
  179. c
  180. c %--------------------%
  181. c | External Functions |
  182. c %--------------------%
  183. c
  184. REAL*8
  185. external dlamch
  186. c
  187. c %----------------------%
  188. **c | Intrinsics Functions |
  189. **c %----------------------%
  190. **c
  191. ** intrinsic abs
  192. **c
  193. **c %----------------%
  194. **c | Data statments |
  195. **c %----------------%
  196. **c
  197. data first / .true. /
  198. **c
  199. **c %-----------------------%
  200. **c | Executable Statements |
  201. c %-----------------------%
  202. c
  203. if (first) then
  204. epsmch = dlamch('Epsilon-Machine')
  205. first = .false.
  206. end if
  207. itop = 1
  208. c
  209. c %-------------------------------%
  210. c | Initialize timing statistics |
  211. c | & message level for debugging |
  212. c %-------------------------------%
  213. c
  214. * call arscnd (t0)
  215. msglvl = msapps
  216. c
  217. kplusp = kev + np
  218. c
  219. c %----------------------------------------------%
  220. c | Initialize Q to the identity matrix of order |
  221. c | kplusp used to accumulate the rotations. |
  222. c %----------------------------------------------%
  223. c
  224. call dlaset ('All', kplusp, kplusp, zero, one, q, ldq)
  225. c
  226. c %----------------------------------------------%
  227. c | Quick return if there are no shifts to apply |
  228. c %----------------------------------------------%
  229. c
  230. if (np .eq. 0) go to 9000
  231. c
  232. c %----------------------------------------------------------%
  233. c | Apply the np shifts implicitly. Apply each shift to the |
  234. c | whole matrix and not just to the submatrix from which it |
  235. c | comes. |
  236. c %----------------------------------------------------------%
  237. c
  238. do 90 jj = 1, np
  239. c
  240. istart = itop
  241. c
  242. c %----------------------------------------------------------%
  243. c | Check for splitting and deflation. Currently we consider |
  244. c | an off-diagonal element h(i+1,1) negligible if |
  245. c | h(i+1,1) .le. epsmch*( |h(i,2)| + |h(i+1,2)| ) |
  246. c | for i=1:KEV+NP-1. |
  247. c | If above condition tests true then we set h(i+1,1) = 0. |
  248. c | Note that h(1:KEV+NP,1) are assumed to be non negative. |
  249. c %----------------------------------------------------------%
  250. c
  251. 20 continue
  252. c
  253. c %------------------------------------------------%
  254. c | The following loop exits early if we encounter |
  255. c | a negligible off diagonal element. |
  256. c %------------------------------------------------%
  257. c
  258. do 30 i = istart, kplusp-1
  259. big = abs(h(i,2)) + abs(h(i+1,2))
  260. if (h(i+1,1) .le. epsmch*big) then
  261. if (msglvl .gt. 0) then
  262. call ivout (logfil, 1, i, ndigit,
  263. & '_sapps: deflation at row/column no.')
  264. call ivout (logfil, 1, jj, ndigit,
  265. & '_sapps: occured before shift number.')
  266. call dvout (logfil, 1, h(i+1,1), ndigit,
  267. & '_sapps: the corresponding off diagonal element')
  268. end if
  269. h(i+1,1) = zero
  270. iend = i
  271. go to 40
  272. end if
  273. 30 continue
  274. iend = kplusp
  275. 40 continue
  276. c
  277. if (istart .lt. iend) then
  278. c
  279. c %--------------------------------------------------------%
  280. c | Construct the plane rotation G'(istart,istart+1,theta) |
  281. c | that attempts to drive h(istart+1,1) to zero. |
  282. c %--------------------------------------------------------%
  283. c
  284. f = h(istart,2) - shift(jj)
  285. g = h(istart+1,1)
  286. call dlartg (f, g, c, s, r)
  287. c
  288. c %-------------------------------------------------------%
  289. c | Apply rotation to the left and right of H; |
  290. c | H <- G' * H * G, where G = G(istart,istart+1,theta). |
  291. c | This will create a "bulge". |
  292. c %-------------------------------------------------------%
  293. c
  294. a1 = c*h(istart,2) + s*h(istart+1,1)
  295. a2 = c*h(istart+1,1) + s*h(istart+1,2)
  296. a4 = c*h(istart+1,2) - s*h(istart+1,1)
  297. a3 = c*h(istart+1,1) - s*h(istart,2)
  298. h(istart,2) = c*a1 + s*a2
  299. h(istart+1,2) = c*a4 - s*a3
  300. h(istart+1,1) = c*a3 + s*a4
  301. c
  302. c %----------------------------------------------------%
  303. c | Accumulate the rotation in the matrix Q; Q <- Q*G |
  304. c %----------------------------------------------------%
  305. c
  306. do 60 j = 1, min(istart+jj,kplusp)
  307. a1 = c*q(j,istart) + s*q(j,istart+1)
  308. q(j,istart+1) = - s*q(j,istart) + c*q(j,istart+1)
  309. q(j,istart) = a1
  310. 60 continue
  311. c
  312. c
  313. c %----------------------------------------------%
  314. c | The following loop chases the bulge created. |
  315. c | Note that the previous rotation may also be |
  316. c | done within the following loop. But it is |
  317. c | kept separate to make the distinction among |
  318. c | the bulge chasing sweeps and the first plane |
  319. c | rotation designed to drive h(istart+1,1) to |
  320. c | zero. |
  321. c %----------------------------------------------%
  322. c
  323. do 70 i = istart+1, iend-1
  324. c
  325. c %----------------------------------------------%
  326. c | Construct the plane rotation G'(i,i+1,theta) |
  327. c | that zeros the i-th bulge that was created |
  328. c | by G(i-1,i,theta). g represents the bulge. |
  329. c %----------------------------------------------%
  330. c
  331. f = h(i,1)
  332. g = s*h(i+1,1)
  333. c
  334. c %----------------------------------%
  335. c | Final update with G(i-1,i,theta) |
  336. c %----------------------------------%
  337. c
  338. h(i+1,1) = c*h(i+1,1)
  339. call dlartg (f, g, c, s, r)
  340. c
  341. c %-------------------------------------------%
  342. c | The following ensures that h(1:iend-1,1), |
  343. c | the first iend-2 off diagonal of elements |
  344. c | H, remain non negative. |
  345. c %-------------------------------------------%
  346. c
  347. if (r .lt. zero) then
  348. r = -r
  349. c = -c
  350. s = -s
  351. end if
  352. c
  353. c %--------------------------------------------%
  354. c | Apply rotation to the left and right of H; |
  355. c | H <- G * H * G', where G = G(i,i+1,theta) |
  356. c %--------------------------------------------%
  357. c
  358. h(i,1) = r
  359. c
  360. a1 = c*h(i,2) + s*h(i+1,1)
  361. a2 = c*h(i+1,1) + s*h(i+1,2)
  362. a3 = c*h(i+1,1) - s*h(i,2)
  363. a4 = c*h(i+1,2) - s*h(i+1,1)
  364. c
  365. h(i,2) = c*a1 + s*a2
  366. h(i+1,2) = c*a4 - s*a3
  367. h(i+1,1) = c*a3 + s*a4
  368. c
  369. c %----------------------------------------------------%
  370. c | Accumulate the rotation in the matrix Q; Q <- Q*G |
  371. c %----------------------------------------------------%
  372. c
  373. do 50 j = 1, min( i+jj, kplusp )
  374. a1 = c*q(j,i) + s*q(j,i+1)
  375. q(j,i+1) = - s*q(j,i) + c*q(j,i+1)
  376. q(j,i) = a1
  377. 50 continue
  378. c
  379. 70 continue
  380. c
  381. end if
  382. c
  383. c %--------------------------%
  384. c | Update the block pointer |
  385. c %--------------------------%
  386. c
  387. istart = iend + 1
  388. c
  389. c %------------------------------------------%
  390. c | Make sure that h(iend,1) is non-negative |
  391. c | If not then set h(iend,1) <-- -h(iend,1) |
  392. c | and negate the last column of Q. |
  393. c | We have effectively carried out a |
  394. c | similarity on transformation H |
  395. c %------------------------------------------%
  396. c
  397. if (h(iend,1) .lt. zero) then
  398. h(iend,1) = -h(iend,1)
  399. call dscal(kplusp, -one, q(1,iend), 1)
  400. end if
  401. c
  402. c %--------------------------------------------------------%
  403. c | Apply the same shift to the next block if there is any |
  404. c %--------------------------------------------------------%
  405. c
  406. if (iend .lt. kplusp) go to 20
  407. c
  408. c %-----------------------------------------------------%
  409. c | Check if we can increase the the start of the block |
  410. c %-----------------------------------------------------%
  411. c
  412. do 80 i = itop, kplusp-1
  413. if (h(i+1,1) .gt. zero) go to 90
  414. itop = itop + 1
  415. 80 continue
  416. c
  417. c %-----------------------------------%
  418. c | Finished applying the jj-th shift |
  419. c %-----------------------------------%
  420. c
  421. 90 continue
  422. c
  423. c %------------------------------------------%
  424. c | All shifts have been applied. Check for |
  425. c | more possible deflation that might occur |
  426. c | after the last shift is applied. |
  427. c %------------------------------------------%
  428. c
  429. do 100 i = itop, kplusp-1
  430. big = abs(h(i,2)) + abs(h(i+1,2))
  431. if (h(i+1,1) .le. epsmch*big) then
  432. if (msglvl .gt. 0) then
  433. call ivout (logfil, 1, i, ndigit,
  434. & '_sapps: deflation at row/column no.')
  435. call dvout (logfil, 1, h(i+1,1), ndigit,
  436. & '_sapps: the corresponding off diagonal element')
  437. end if
  438. h(i+1,1) = zero
  439. end if
  440. 100 continue
  441. c
  442. c %-------------------------------------------------%
  443. c | Compute the (kev+1)-st column of (V*Q) and |
  444. c | temporarily store the result in WORKD(N+1:2*N). |
  445. c | This is not necessary if h(kev+1,1) = 0. |
  446. c %-------------------------------------------------%
  447. c
  448. if ( h(kev+1,1) .gt. zero )
  449. & call dgemv ('N', n, kplusp, one, v, ldv,
  450. & q(1,kev+1), 1, zero, workd(n+1), 1)
  451. c
  452. c %-------------------------------------------------------%
  453. c | Compute column 1 to kev of (V*Q) in backward order |
  454. c | taking advantage that Q is an upper triangular matrix |
  455. c | with lower bandwidth np. |
  456. c | Place results in v(:,kplusp-kev:kplusp) temporarily. |
  457. c %-------------------------------------------------------%
  458. c
  459. do 130 i = 1, kev
  460. call dgemv ('N', n, kplusp-i+1, one, v, ldv,
  461. & q(1,kev-i+1), 1, zero, workd, 1)
  462. call dcopy (n, workd, 1, v(1,kplusp-i+1), 1)
  463. 130 continue
  464. c
  465. c %-------------------------------------------------%
  466. c | Move v(:,kplusp-kev+1:kplusp) into v(:,1:kev). |
  467. c %-------------------------------------------------%
  468. c
  469. call dlacpy ('All', n, kev, v(1,np+1), ldv, v, ldv)
  470. c
  471. c %--------------------------------------------%
  472. c | Copy the (kev+1)-st column of (V*Q) in the |
  473. c | appropriate place if h(kev+1,1) .ne. zero. |
  474. c %--------------------------------------------%
  475. c
  476. if ( h(kev+1,1) .gt. zero )
  477. & call dcopy (n, workd(n+1), 1, v(1,kev+1), 1)
  478. c
  479. c %-------------------------------------%
  480. c | Update the residual vector: |
  481. c | r <- sigmak*r + betak*v(:,kev+1) |
  482. c | where |
  483. c | sigmak = (e_{kev+p}'*Q)*e_{kev} |
  484. c | betak = e_{kev+1}'*H*e_{kev} |
  485. c %-------------------------------------%
  486. c
  487. call dscal (n, q(kplusp,kev), resid, 1)
  488. if (h(kev+1,1) .gt. zero)
  489. & call daxpy (n, h(kev+1,1), v(1,kev+1), 1, resid, 1)
  490. c
  491. if (msglvl .gt. 1) then
  492. call dvout (logfil, 1, q(kplusp,kev), ndigit,
  493. & '_sapps: sigmak of the updated residual vector')
  494. call dvout (logfil, 1, h(kev+1,1), ndigit,
  495. & '_sapps: betak of the updated residual vector')
  496. call dvout (logfil, kev, h(1,2), ndigit,
  497. & '_sapps: updated main diagonal of H for next iteration')
  498. if (kev .gt. 1) then
  499. call dvout (logfil, kev-1, h(2,1), ndigit,
  500. & '_sapps: updated sub diagonal of H for next iteration')
  501. end if
  502. end if
  503. c
  504. * call arscnd (t1)
  505. tsapps = tsapps + (t1 - t0)
  506. c
  507. 9000 continue
  508. return
  509. c
  510. c %---------------%
  511. c | End of dsapps |
  512. c %---------------%
  513. c
  514. end
  515.  
  516.  
  517.  

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