Télécharger dnapps.eso

Retour à la liste

Numérotation des lignes :

  1. C DNAPPS SOURCE BP208322 15/12/17 21:15:09 8750
  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.
  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. -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. c
  187. c %----------------------%
  188. c | External Subroutines |
  189. c %----------------------%
  190. c
  191. & dlaset, dlabad, arscnd, dlartg
  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) <= 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. 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. if (msglvl .gt. 2 ) then
  267. call ivout (logfil, 1, jj, ndigit,
  268. & '_napps: shift number.')
  269. call dvout (logfil, 1, sigmar, ndigit,
  270. & '_napps: The real part of the shift ')
  271. call dvout (logfil, 1, sigmai, ndigit,
  272. & '_napps: The imaginary part of the shift ')
  273. 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 < 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 (logfil, 1, i, ndigit,
  337. & '_napps: matrix splitting at row/column no.')
  338. call ivout (logfil, 1, jj, ndigit,
  339. & '_napps: matrix splitting with shift number.')
  340. call dvout (logfil, 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. if (msglvl .gt. 2) then
  352. call ivout (logfil, 1, istart, ndigit,
  353. & '_napps: Start of current block ')
  354. call ivout (logfil, 1, iend, ndigit,
  355. & '_napps: End of current block ')
  356. 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 <- 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 <- 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 <- 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 <- 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 <- 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. if (msglvl .gt. 1) then
  622. call dvout (logfil, 1, q(kplusp,kev), ndigit,
  623. & '_napps: sigmak = (e_{kev+p}^T*Q)*e_{kev}')
  624. call dvout (logfil, 1, h(kev+1,kev), ndigit,
  625. & '_napps: betak = e_{kev+1}^T*H*e_{kev}')
  626. call ivout (logfil, 1, kev, ndigit,
  627. & '_napps: Order of the final Hessenberg matrix ')
  628. if (msglvl .gt. 2) then
  629. call dmout (logfil, kev, kev, h, ldh, ndigit,
  630. & '_napps: updated Hessenberg matrix H for next iteration')
  631. end if
  632. c
  633. end if
  634. c
  635. 9000 continue
  636. * call arscnd (t1)
  637. 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.  

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