Numérotation des lignes :

1. C DNAPPS SOURCE BP208322 20/02/06 21:15:19 10512
2. c-----------------------------------------------------------------------
3. c\BeginDoc
4. c
5. c\Name: dnapps
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 implicit shifts 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 which is the product of rotations
17. c and reflections resulting from the NP bulge chage sweeps.
18. c The updated Arnoldi 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 dnapps
24. c ( N, KEV, NP, SHIFTR, SHIFTI, V, LDV, H, LDH, RESID, Q, LDQ,
25. c WORKL, WORKD )
26. c
27. c\Arguments
28. c N Integer. (INPUT)
29. c Problem size, i.e. size of matrix A.
30. c
31. c KEV Integer. (INPUT/OUTPUT)
32. c KEV+NP is the size of the input matrix H.
33. c KEV is the size of the updated matrix HNEW. KEV is only
34. c updated on ouput when fewer than NP shifts are applied in
35. c order to keep the conjugate pair together.
36. c
37. c NP Integer. (INPUT)
38. c Number of implicit shifts to be applied.
39. c
40. c SHIFTR, Double precision array of length NP. (INPUT)
41. c SHIFTI Real and imaginary part of the shifts to be applied.
42. c Upon, entry to dnapps, the shifts must be sorted so that the
43. c conjugate pairs are in consecutive locations.
44. c
45. c V REAL*8 N by (KEV+NP) array. (INPUT/OUTPUT)
46. c On INPUT, V contains the current KEV+NP Arnoldi vectors.
47. c On OUTPUT, V contains the updated KEV Arnoldi vectors
48. c in the first KEV columns of V.
49. c
50. c LDV Integer. (INPUT)
51. c Leading dimension of V exactly as declared in the calling
52. c program.
53. c
54. c H REAL*8 (KEV+NP) by (KEV+NP) array. (INPUT/OUTPUT)
55. c On INPUT, H contains the current KEV+NP by KEV+NP upper
56. c Hessenber matrix of the Arnoldi factorization.
57. c On OUTPUT, H contains the updated KEV by KEV upper Hessenberg
58. c matrix in the KEV leading submatrix.
59. c
60. c LDH Integer. (INPUT)
61. c Leading dimension of H exactly as declared in the calling
62. c program.
63. c
64. c RESID Double precision array of length N. (INPUT/OUTPUT)
65. c On INPUT, RESID contains the the residual vector r_{k+p}.
66. c On OUTPUT, RESID is the update residual vector rnew_{k}
67. c in the first KEV locations.
68. c
69. c Q REAL*8 KEV+NP by KEV+NP work array. (WORKSPACE)
70. c Work array used to accumulate the rotations and reflections
71. c during the bulge chase sweep.
72. c
73. c LDQ Integer. (INPUT)
74. c Leading dimension of Q exactly as declared in the calling
75. c program.
76. c
77. c WORKL Double precision work array of length (KEV+NP). (WORKSPACE)
78. c Private (replicated) array on each PE or array allocated on
79. c the front end.
80. c
81. c WORKD Double precision work array of length 2*N. (WORKSPACE)
82. c Distributed array used in the application of the accumulated
83. c orthogonal matrix Q.
84. c
85. c\EndDoc
86. c
87. c-----------------------------------------------------------------------
88. c
89. c\BeginLib
90. c
91. c\Local variables:
92. c xxxxxx real
93. c
94. c\References:
95. c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
96. c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
97. c pp 357-385.
98. c
99. c\Routines called:
100. c ivout ARPACK utility routine that prints integers.
101. c arscnd ARPACK utility routine for timing. -> deleted by BP in 2020
102. c dmout ARPACK utility routine that prints matrices.
103. c dvout ARPACK utility routine that prints vectors.
104. c dlabad LAPACK routine that computes machine constants.
105. c dlacpy LAPACK matrix copy routine.
106. c dlamch LAPACK routine that determines machine constants.
107. c dlanhs LAPACK routine that computes various norms of a matrix.
108. c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
109. c dlarf LAPACK routine that applies Householder reflection to
110. c a matrix.
111. c dlarfg LAPACK Householder reflection construction routine.
112. c dlartg LAPACK Givens rotation construction routine.
113. c dlaset LAPACK matrix initialization routine.
114. c dgemv Level 2 BLAS routine for matrix vector multiplication.
115. c daxpy Level 1 BLAS that computes a vector triad.
116. c dcopy Level 1 BLAS that copies one vector to another .
117. c dscal Level 1 BLAS that scales a vector.
118. c
119. c\Author
120. c Danny Sorensen Phuong Vu
121. c Richard Lehoucq CRPC / Rice University
122. c Dept. of Computational & Houston, Texas
123. c Applied Mathematics
124. c Rice University
125. c Houston, Texas
126. c
127. c\Revision history:
128. c xx/xx/92: Version ' 2.4'
129. c
130. c\SCCS Information: @(#)
131. c FILE: napps.F SID: 2.4 DATE OF SID: 3/28/97 RELEASE: 2
132. c
133. c\Remarks
134. c 1. In this version, each shift is applied to all the sublocks of
135. c the Hessenberg matrix H and not just to the submatrix that it
136. c comes from. Deflation as in LAPACK routine dlahqr (QR algorithm
137. c for upper Hessenberg matrices ) is used.
138. c The subdiagonals of H are enforced to be non-negative.
139. c
140. c\EndLib
141. c
142. c-----------------------------------------------------------------------
143. c
144. subroutine dnapps
145. & ( n, kev, np, shiftr, shifti, v, ldv, h, ldh, resid, q, ldq,
146. & workl, workd )
147. c
148. c %----------------------------------------------------%
149. c | Include files for debugging and timing information |
150. c -INC TARTRAK
151. c %----------------------------------------------------%
152. c
153. c
154. c %------------------%
155. c | Scalar Arguments |
156. c %------------------%
157. c
158. integer kev, ldh, ldq, ldv, n, np
159. c
160. c %-----------------%
161. c | Array Arguments |
162. c %-----------------%
163. c
164. REAL*8
165. & h(ldh,kev+np), resid(n), shifti(np), shiftr(np),
166. & v(ldv,kev+np), q(ldq,kev+np), workd(2*n), workl(kev+np)
167. c
168. c %------------%
169. c | Parameters |
170. c %------------%
171. c
172. REAL*8
173. & one, zero
174. parameter (one = 1.0D+0, zero = 0.0D+0)
175. c
176. c %------------------------%
177. c | Local Scalars & Arrays |
178. c %------------------------%
179. c
180. integer i, iend, ir, istart, j, jj, kplusp, msglvl, nr
181. logical cconj, first
182. REAL*8
183. & c, f, g, h11, h12, h21, h22, h32, ovfl, r, s, sigmai,
184. & sigmar, smlnum, ulp, unfl, u(3), t, tau, tst1
185. save first, ovfl, smlnum, ulp, unfl
186. parameter (msglvl=0)
187. c
188. c %----------------------%
189. c | External Subroutines |
190. c %----------------------%
191. c
192. c
193. c %--------------------%
194. c | External Functions |
195. c %--------------------%
196. c
197. REAL*8
198. external dlamch, dlanhs, dlapy2
199. c
200. c %----------------------%
201. **c | Intrinsics Functions |
202. **c %----------------------%
203. **c
204. ** intrinsic abs, max, min
205. **c
206. **c %----------------%
207. **c | Data statments |
208. **c %----------------%
209. **c
210. data first / .true. /
211. **c
212. **c %-----------------------%
213. **c | Executable Statements |
214. c %-----------------------%
215. c
216. if (first) then
217. c
218. c %-----------------------------------------------%
219. c | Set machine-dependent constants for the |
220. c | stopping criterion. If norm(H) &lt;= sqrt(OVFL), |
221. c | overflow should not occur. |
222. c | REFERENCE: LAPACK subroutine dlahqr |
223. c %-----------------------------------------------%
224. c
225. unfl = dlamch( 'safe minimum' )
226. ovfl = one / unfl
227. call dlabad( unfl, ovfl )
228. ulp = dlamch( 'precision' )
229. smlnum = unfl*( n / ulp )
230. first = .false.
231. end if
232. c
233. c %-------------------------------%
234. c | Initialize timing statistics |
235. c | & message level for debugging |
236. c %-------------------------------%
237. c
238. * call arscnd (t0)
239. c msglvl = mnapps
240. kplusp = kev + np
241. c
242. c %--------------------------------------------%
243. c | Initialize Q to the identity to accumulate |
244. c | the rotations and reflections |
245. c %--------------------------------------------%
246. c
247. call dlaset ('All', kplusp, kplusp, zero, one, q, ldq)
248. c
249. c %----------------------------------------------%
250. c | Quick return if there are no shifts to apply |
251. c %----------------------------------------------%
252. c
253. if (np .eq. 0) go to 9000
254. c
255. c %----------------------------------------------%
256. c | Chase the bulge with the application of each |
257. c | implicit shift. Each shift is applied to the |
258. c | whole matrix including each block. |
259. c %----------------------------------------------%
260. c
261. cconj = .false.
262. do 110 jj = 1, np
263. sigmar = shiftr(jj)
264. sigmai = shifti(jj)
265. c
266. c if (msglvl .gt. 2 ) then
267. c call ivout ( 1, jj, ndigit,
268. c & '_napps: shift number.')
269. c call dvout ( 1, sigmar, ndigit,
270. c & '_napps: The real part of the shift ')
271. c call dvout ( 1, sigmai, ndigit,
272. c & '_napps: The imaginary part of the shift ')
273. c end if
274. c
275. c %-------------------------------------------------%
276. c | The following set of conditionals is necessary |
277. c | in order that complex conjugate pairs of shifts |
278. c | are applied together or not at all. |
279. c %-------------------------------------------------%
280. c
281. if ( cconj ) then
282. c
283. c %-----------------------------------------%
284. c | cconj = .true. means the previous shift |
285. c | had non-zero imaginary part. |
286. c %-----------------------------------------%
287. c
288. cconj = .false.
289. go to 110
290. else if ( jj .lt. np .and. abs( sigmai ) .gt. zero ) then
291. c
292. c %------------------------------------%
293. c | Start of a complex conjugate pair. |
294. c %------------------------------------%
295. c
296. cconj = .true.
297. else if ( jj .eq. np .and. abs( sigmai ) .gt. zero ) then
298. c
299. c %----------------------------------------------%
300. c | The last shift has a nonzero imaginary part. |
301. c | Don't apply it; thus the order of the |
302. c | compressed H is order KEV+1 since only np-1 |
303. c | were applied. |
304. c %----------------------------------------------%
305. c
306. kev = kev + 1
307. go to 110
308. end if
309. istart = 1
310. 20 continue
311. c
312. c %--------------------------------------------------%
313. c | if sigmai = 0 then |
314. c | Apply the jj-th shift ... |
315. c | else |
316. c | Apply the jj-th and (jj+1)-th together ... |
317. c | (Note that jj &lt; np at this point in the code) |
318. c | end |
319. c | to the current block of H. The next do loop |
320. c | determines the current block ; |
321. c %--------------------------------------------------%
322. c
323. do 30 i = istart, kplusp-1
324. c
325. c %----------------------------------------%
326. c | Check for splitting and deflation. Use |
327. c | a standard test as in the QR algorithm |
328. c | REFERENCE: LAPACK subroutine dlahqr |
329. c %----------------------------------------%
330. c
331. tst1 = abs( h( i, i ) ) + abs( h( i+1, i+1 ) )
332. if( tst1.eq.zero )
333. & tst1 = dlanhs( '1', kplusp-jj+1, h, ldh, workl )
334. if( abs( h( i+1,i ) ).le.max( ulp*tst1, smlnum ) ) then
335. if (msglvl .gt. 0) then
336. call ivout ( 1, i, ndigit,
337. & '_napps: matrix splitting at row/column no.')
338. call ivout ( 1, jj, ndigit,
339. & '_napps: matrix splitting with shift number.')
340. call dvout ( 1, h(i+1,i), ndigit,
341. & '_napps: off diagonal element.')
342. end if
343. iend = i
344. h(i+1,i) = zero
345. go to 40
346. end if
347. 30 continue
348. iend = kplusp
349. 40 continue
350. c
351. c if (msglvl .gt. 2) then
352. c call ivout ( 1, istart, ndigit,
353. c & '_napps: Start of current block ')
354. c call ivout ( 1, iend, ndigit,
355. c & '_napps: End of current block ')
356. c end if
357. c
358. c %------------------------------------------------%
359. c | No reason to apply a shift to block of order 1 |
360. c %------------------------------------------------%
361. c
362. if ( istart .eq. iend ) go to 100
363. c
364. c %------------------------------------------------------%
365. c | If istart + 1 = iend then no reason to apply a |
366. c | complex conjugate pair of shifts on a 2 by 2 matrix. |
367. c %------------------------------------------------------%
368. c
369. if ( istart + 1 .eq. iend .and. abs( sigmai ) .gt. zero )
370. & go to 100
371. c
372. h11 = h(istart,istart)
373. h21 = h(istart+1,istart)
374. if ( abs( sigmai ) .le. zero ) then
375. c
376. c %---------------------------------------------%
377. c | Real-valued shift ==> apply single shift QR |
378. c %---------------------------------------------%
379. c
380. f = h11 - sigmar
381. g = h21
382. c
383. do 80 i = istart, iend-1
384. c
385. c %-----------------------------------------------------%
386. c | Contruct the plane rotation G to zero out the bulge |
387. c %-----------------------------------------------------%
388. c
389. call dlartg (f, g, c, s, r)
390. if (i .gt. istart) then
391. c
392. c %-------------------------------------------%
393. c | The following ensures that h(1:iend-1,1), |
394. c | the first iend-2 off diagonal of elements |
395. c | H, remain non negative. |
396. c %-------------------------------------------%
397. c
398. if (r .lt. zero) then
399. r = -r
400. c = -c
401. s = -s
402. end if
403. h(i,i-1) = r
404. h(i+1,i-1) = zero
405. end if
406. c
407. c %---------------------------------------------%
408. c | Apply rotation to the left of H; H &lt;- G'*H |
409. c %---------------------------------------------%
410. c
411. do 50 j = i, kplusp
412. t = c*h(i,j) + s*h(i+1,j)
413. h(i+1,j) = -s*h(i,j) + c*h(i+1,j)
414. h(i,j) = t
415. 50 continue
416. c
417. c %---------------------------------------------%
418. c | Apply rotation to the right of H; H &lt;- H*G |
419. c %---------------------------------------------%
420. c
421. do 60 j = 1, min(i+2,iend)
422. t = c*h(j,i) + s*h(j,i+1)
423. h(j,i+1) = -s*h(j,i) + c*h(j,i+1)
424. h(j,i) = t
425. 60 continue
426. c
427. c %----------------------------------------------------%
428. c | Accumulate the rotation in the matrix Q; Q &lt;- Q*G |
429. c %----------------------------------------------------%
430. c
431. do 70 j = 1, min( i+jj, kplusp )
432. t = c*q(j,i) + s*q(j,i+1)
433. q(j,i+1) = - s*q(j,i) + c*q(j,i+1)
434. q(j,i) = t
435. 70 continue
436. c
437. c %---------------------------%
438. c | Prepare for next rotation |
439. c %---------------------------%
440. c
441. if (i .lt. iend-1) then
442. f = h(i+1,i)
443. g = h(i+2,i)
444. end if
445. 80 continue
446. c
447. c %-----------------------------------%
448. c | Finished applying the real shift. |
449. c %-----------------------------------%
450. c
451. else
452. c
453. c %----------------------------------------------------%
454. c | Complex conjugate shifts ==> apply double shift QR |
455. c %----------------------------------------------------%
456. c
457. h12 = h(istart,istart+1)
458. h22 = h(istart+1,istart+1)
459. h32 = h(istart+2,istart+1)
460. c
461. c %---------------------------------------------------------%
462. c | Compute 1st column of (H - shift*I)*(H - conj(shift)*I) |
463. c %---------------------------------------------------------%
464. c
465. s = 2.0*sigmar
466. t = dlapy2 ( sigmar, sigmai )
467. u(1) = ( h11 * (h11 - s) + t * t ) / h21 + h12
468. u(2) = h11 + h22 - s
469. u(3) = h32
470. c
471. do 90 i = istart, iend-1
472. c
473. nr = min ( 3, iend-i+1 )
474. c
475. c %-----------------------------------------------------%
476. c | Construct Householder reflector G to zero out u(1). |
477. c | G is of the form I - tau*( 1 u )' * ( 1 u' ). |
478. c %-----------------------------------------------------%
479. c
480. call dlarfg ( nr, u(1), u(2), 1, tau )
481. c
482. if (i .gt. istart) then
483. h(i,i-1) = u(1)
484. h(i+1,i-1) = zero
485. if (i .lt. iend-1) h(i+2,i-1) = zero
486. end if
487. u(1) = one
488. c
489. c %--------------------------------------%
490. c | Apply the reflector to the left of H |
491. c %--------------------------------------%
492. c
493. call dlarf ('Left', nr, kplusp-i+1, u, 1, tau,
494. & h(i,i), ldh, workl)
495. c
496. c %---------------------------------------%
497. c | Apply the reflector to the right of H |
498. c %---------------------------------------%
499. c
500. ir = min ( i+3, iend )
501. call dlarf ('Right', ir, nr, u, 1, tau,
502. & h(1,i), ldh, workl)
503. c
504. c %-----------------------------------------------------%
505. c | Accumulate the reflector in the matrix Q; Q &lt;- Q*G |
506. c %-----------------------------------------------------%
507. c
508. call dlarf ('Right', kplusp, nr, u, 1, tau,
509. & q(1,i), ldq, workl)
510. c
511. c %----------------------------%
512. c | Prepare for next reflector |
513. c %----------------------------%
514. c
515. if (i .lt. iend-1) then
516. u(1) = h(i+1,i)
517. u(2) = h(i+2,i)
518. if (i .lt. iend-2) u(3) = h(i+3,i)
519. end if
520. c
521. 90 continue
522. c
523. c %--------------------------------------------%
524. c | Finished applying a complex pair of shifts |
525. c | to the current block |
526. c %--------------------------------------------%
527. c
528. end if
529. c
530. 100 continue
531. c
532. c %---------------------------------------------------------%
533. c | Apply the same shift to the next block if there is any. |
534. c %---------------------------------------------------------%
535. c
536. istart = iend + 1
537. if (iend .lt. kplusp) go to 20
538. c
539. c %---------------------------------------------%
540. c | Loop back to the top to get the next shift. |
541. c %---------------------------------------------%
542. c
543. 110 continue
544. c
545. c %--------------------------------------------------%
546. c | Perform a similarity transformation that makes |
547. c | sure that H will have non negative sub diagonals |
548. c %--------------------------------------------------%
549. c
550. do 120 j=1,kev
551. if ( h(j+1,j) .lt. zero ) then
552. call dscal( kplusp-j+1, -one, h(j+1,j), ldh )
553. call dscal( min(j+2, kplusp), -one, h(1,j+1), 1 )
554. call dscal( min(j+np+1,kplusp), -one, q(1,j+1), 1 )
555. end if
556. 120 continue
557. c
558. do 130 i = 1, kev
559. c
560. c %--------------------------------------------%
561. c | Final check for splitting and deflation. |
562. c | Use a standard test as in the QR algorithm |
563. c | REFERENCE: LAPACK subroutine dlahqr |
564. c %--------------------------------------------%
565. c
566. tst1 = abs( h( i, i ) ) + abs( h( i+1, i+1 ) )
567. if( tst1.eq.zero )
568. & tst1 = dlanhs( '1', kev, h, ldh, workl )
569. if( h( i+1,i ) .le. max( ulp*tst1, smlnum ) )
570. & h(i+1,i) = zero
571. 130 continue
572. c
573. c %-------------------------------------------------%
574. c | Compute the (kev+1)-st column of (V*Q) and |
575. c | temporarily store the result in WORKD(N+1:2*N). |
576. c | This is needed in the residual update since we |
577. c | cannot GUARANTEE that the corresponding entry |
578. c | of H would be zero as in exact arithmetic. |
579. c %-------------------------------------------------%
580. c
581. if (h(kev+1,kev) .gt. zero)
582. & call dgemv ('N', n, kplusp, one, v, ldv, q(1,kev+1), 1, zero,
583. & workd(n+1), 1)
584. c
585. c %----------------------------------------------------------%
586. c | Compute column 1 to kev of (V*Q) in backward order |
587. c | taking advantage of the upper Hessenberg structure of Q. |
588. c %----------------------------------------------------------%
589. c
590. do 140 i = 1, kev
591. call dgemv ('N', n, kplusp-i+1, one, v, ldv,
592. & q(1,kev-i+1), 1, zero, workd, 1)
593. call dcopy (n, workd, 1, v(1,kplusp-i+1), 1)
594. 140 continue
595. c
596. c %-------------------------------------------------%
597. c | Move v(:,kplusp-kev+1:kplusp) into v(:,1:kev). |
598. c %-------------------------------------------------%
599. c
600. call dlacpy ('A', n, kev, v(1,kplusp-kev+1), ldv, v, ldv)
601. c
602. c %--------------------------------------------------------------%
603. c | Copy the (kev+1)-st column of (V*Q) in the appropriate place |
604. c %--------------------------------------------------------------%
605. c
606. if (h(kev+1,kev) .gt. zero)
607. & call dcopy (n, workd(n+1), 1, v(1,kev+1), 1)
608. c
609. c %-------------------------------------%
610. c | Update the residual vector: |
611. c | r &lt;- sigmak*r + betak*v(:,kev+1) |
612. c | where |
613. c | sigmak = (e_{kplusp}'*Q)*e_{kev} |
614. c | betak = e_{kev+1}'*H*e_{kev} |
615. c %-------------------------------------%
616. c
617. call dscal (n, q(kplusp,kev), resid, 1)
618. if (h(kev+1,kev) .gt. zero)
619. & call daxpy (n, h(kev+1,kev), v(1,kev+1), 1, resid, 1)
620. c
621. c if (msglvl .gt. 1) then
622. c call dvout ( 1, q(kplusp,kev), ndigit,
623. c & '_napps: sigmak = (e_{kev+p}^T*Q)*e_{kev}')
624. c call dvout ( 1, h(kev+1,kev), ndigit,
625. c & '_napps: betak = e_{kev+1}^T*H*e_{kev}')
626. c call ivout ( 1, kev, ndigit,
627. c & '_napps: Order of the final Hessenberg matrix ')
628. c if (msglvl .gt. 2) then
629. c call dmout (logfil, kev, kev, h, ldh, ndigit,
630. c & '_napps: updated Hessenberg matrix H for next iteration')
631. c end if
632. c c
633. c end if
634. c
635. 9000 continue
636. * call arscnd (t1)
637. c tnapps = tnapps + (t1 - t0)
638. c
639. return
640. c
641. c %---------------%
642. c | End of dnapps |
643. c %---------------%
644. c
645. end
646.
647.
648.
649.

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