Télécharger dnaup2.eso

Retour à la liste

Numérotation des lignes :

  1. C DNAUP2 SOURCE BP208322 15/10/13 21:15:44 8670
  2. c\BeginDoc
  3. c
  4. c\Name: dnaup2
  5. c
  6. c\Description:
  7. c Intermediate level interface called by dnaupd .
  8. c
  9. c\Usage:
  10. c call dnaup2
  11. c ( IDO, BMAT, N, WHICH, NEV, NP, TOL, RESID, MODE, IUPD,
  12. c ISHIFT, MXITER, V, LDV, H, LDH, RITZR, RITZI, BOUNDS,
  13. c Q, LDQ, WORKL, IPNTR, WORKD, INFO )
  14. c
  15. c\Arguments
  16. c
  17. c IDO, BMAT, N, WHICH, NEV, TOL, RESID: same as defined in dnaupd .
  18. c MODE, ISHIFT, MXITER: see the definition of IPARAM in dnaupd .
  19. c
  20. c NP Integer. (INPUT/OUTPUT)
  21. c Contains the number of implicit shifts to apply during
  22. c each Arnoldi iteration.
  23. c If ISHIFT=1, NP is adjusted dynamically at each iteration
  24. c to accelerate convergence and prevent stagnation.
  25. c This is also roughly equal to the number of matrix-vector
  26. c products (involving the operator OP) per Arnoldi iteration.
  27. c The logic for adjusting is contained within the current
  28. c subroutine.
  29. c If ISHIFT=0, NP is the number of shifts the user needs
  30. c to provide via reverse comunication. 0 < NP < NCV-NEV.
  31. c NP may be less than NCV-NEV for two reasons. The first, is
  32. c to keep complex conjugate pairs of "wanted" Ritz values
  33. c together. The second, is that a leading block of the current
  34. c upper Hessenberg matrix has split off and contains "unwanted"
  35. c Ritz values.
  36. c Upon termination of the IRA iteration, NP contains the number
  37. c of "converged" wanted Ritz values.
  38. c
  39. c IUPD Integer. (INPUT)
  40. c IUPD .EQ. 0: use explicit restart instead implicit update.
  41. c IUPD .NE. 0: use implicit update.
  42. c
  43. c V REAL*8 N by (NEV+NP) array. (INPUT/OUTPUT)
  44. c The Arnoldi basis vectors are returned in the first NEV
  45. c columns of V.
  46. c
  47. c LDV Integer. (INPUT)
  48. c Leading dimension of V exactly as declared in the calling
  49. c program.
  50. c
  51. c H REAL*8 (NEV+NP) by (NEV+NP) array. (OUTPUT)
  52. c H is used to store the generated upper Hessenberg matrix
  53. c
  54. c LDH Integer. (INPUT)
  55. c Leading dimension of H exactly as declared in the calling
  56. c program.
  57. c
  58. c RITZR, Double precision arrays of length NEV+NP. (OUTPUT)
  59. c RITZI RITZR(1:NEV) (resp. RITZI(1:NEV)) contains the real (resp.
  60. c imaginary) part of the computed Ritz values of OP.
  61. c
  62. c BOUNDS Double precision array of length NEV+NP. (OUTPUT)
  63. c BOUNDS(1:NEV) contain the error bounds corresponding to
  64. c the computed Ritz values.
  65. c
  66. c Q REAL*8 (NEV+NP) by (NEV+NP) array. (WORKSPACE)
  67. c Private (replicated) work array used to accumulate the
  68. c rotation in the shift application step.
  69. c
  70. c LDQ Integer. (INPUT)
  71. c Leading dimension of Q exactly as declared in the calling
  72. c program.
  73. c
  74. c WORKL Double precision work array of length at least
  75. c (NEV+NP)**2 + 3*(NEV+NP). (INPUT/WORKSPACE)
  76. c Private (replicated) array on each PE or array allocated on
  77. c the front end. It is used in shifts calculation, shifts
  78. c application and convergence checking.
  79. c
  80. c On exit, the last 3*(NEV+NP) locations of WORKL contain
  81. c the Ritz values (real,imaginary) and associated Ritz
  82. c estimates of the current Hessenberg matrix. They are
  83. c listed in the same order as returned from dneigh .
  84. c
  85. c If ISHIFT .EQ. O and IDO .EQ. 3, the first 2*NP locations
  86. c of WORKL are used in reverse communication to hold the user
  87. c supplied shifts.
  88. c
  89. c IPNTR Integer array of length 3. (OUTPUT)
  90. c Pointer to mark the starting locations in the WORKD for
  91. c vectors used by the Arnoldi iteration.
  92. c -------------------------------------------------------------
  93. c IPNTR(1): pointer to the current operand vector X.
  94. c IPNTR(2): pointer to the current result vector Y.
  95. c IPNTR(3): pointer to the vector B * X when used in the
  96. c shift-and-invert mode. X is the current operand.
  97. c -------------------------------------------------------------
  98. c
  99. c WORKD Double precision work array of length 3*N. (WORKSPACE)
  100. c Distributed array to be used in the basic Arnoldi iteration
  101. c for reverse communication. The user should not use WORKD
  102. c as temporary workspace during the iteration !!!!!!!!!!
  103. c See Data Distribution Note in DNAUPD.
  104. c
  105. c INFO Integer. (INPUT/OUTPUT)
  106. c If INFO .EQ. 0, a randomly initial residual vector is used.
  107. c If INFO .NE. 0, RESID contains the initial residual vector,
  108. c possibly from a previous run.
  109. c Error flag on output.
  110. c = 0: Normal return.
  111. c = 1: Maximum number of iterations taken.
  112. c All possible eigenvalues of OP has been found.
  113. c NP returns the number of converged Ritz values.
  114. c = 2: No shifts could be applied.
  115. c = -8: Error return from LAPACK eigenvalue calculation;
  116. c This should never happen.
  117. c = -9: Starting vector is zero.
  118. c = -9999: Could not build an Arnoldi factorization.
  119. c Size that was built in returned in NP.
  120. c
  121. c\EndDoc
  122. c
  123. c-----------------------------------------------------------------------
  124. c
  125. c\BeginLib
  126. c
  127. c\Local variables:
  128. c xxxxxx real
  129. c
  130. c\References:
  131. c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
  132. c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
  133. c pp 357-385.
  134. c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
  135. c Restarted Arnoldi Iteration", Rice University Technical Report
  136. c TR95-13, Department of Computational and Applied Mathematics.
  137. c
  138. c\Routines called:
  139. c dgetv0 ARPACK initial vector generation routine.
  140. c dnaitr ARPACK Arnoldi factorization routine.
  141. c dnapps ARPACK application of implicit shifts routine.
  142. c dnconv ARPACK convergence of Ritz values routine.
  143. c dneigh ARPACK compute Ritz values and error bounds routine.
  144. c dngets ARPACK reorder Ritz values and error bounds routine.
  145. c dsortc ARPACK sorting routine.
  146. c ivout ARPACK utility routine that prints integers.
  147. c arscnd ARPACK utility routine for timing.
  148. c dmout ARPACK utility routine that prints matrices
  149. c dvout ARPACK utility routine that prints vectors.
  150. c dlamch LAPACK routine that determines machine constants.
  151. c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
  152. c dcopy Level 1 BLAS that copies one vector to another .
  153. c ddot Level 1 BLAS that computes the scalar product of two vectors.
  154. c dnrm2 Level 1 BLAS that computes the norm of a vector.
  155. c dswap Level 1 BLAS that swaps two vectors.
  156. c
  157. c\Author
  158. c Danny Sorensen Phuong Vu
  159. c Richard Lehoucq CRPC / Rice University
  160. c Dept. of Computational & Houston, Texas
  161. c Applied Mathematics
  162. c Rice University
  163. c Houston, Texas
  164. c
  165. c\SCCS Information: @(#)
  166. c FILE: naup2.F SID: 2.8 DATE OF SID: 10/17/00 RELEASE: 2
  167. c
  168. c\Remarks
  169. c 1. None
  170. c
  171. c\EndLib
  172. c
  173. c-----------------------------------------------------------------------
  174. c
  175. subroutine dnaup2
  176. & ( ido, bmat, n, which, nev, np, tol, resid, mode, iupd,
  177. & ishift, mxiter, v, ldv, h, ldh, ritzr, ritzi, bounds,
  178. & q, ldq, workl, ipntr, workd, info )
  179. c
  180. c %----------------------------------------------------%
  181. c | Include files for debugging and timing information |
  182. -INC TARTRAK
  183. c %----------------------------------------------------%
  184. c
  185. c
  186. c %------------------%
  187. c | Scalar Arguments |
  188. c %------------------%
  189. c
  190. character bmat*1, which*2
  191. integer ido, info, ishift, iupd, mode, ldh, ldq, ldv, mxiter,
  192. & n, nev, np
  193. REAL*8
  194. & tol
  195. c
  196. c %-----------------%
  197. c | Array Arguments |
  198. c %-----------------%
  199. c
  200. integer ipntr(13)
  201. REAL*8
  202. & bounds(nev+np), h(ldh,nev+np), q(ldq,nev+np), resid(n),
  203. & ritzi(nev+np), ritzr(nev+np), v(ldv,nev+np),
  204. & workd(3*n), workl( (nev+np)*(nev+np+3) )
  205. c
  206. c %------------%
  207. c | Parameters |
  208. c %------------%
  209. c
  210. REAL*8
  211. & one, zero
  212. parameter (one = 1.0D+0 , zero = 0.0D+0 )
  213. c
  214. c %---------------%
  215. c | Local Scalars |
  216. c %---------------%
  217. c
  218. character wprime*2
  219. logical cnorm , getv0, initv, update, ushift
  220. integer ierr , iter , j , kplusp, msglvl, nconv,
  221. & nevbef, nev0 , np0 , nptemp, numcnv
  222. REAL*8
  223. & rnorm , temp , eps23
  224. save cnorm , getv0, initv, update, ushift,
  225. & rnorm , iter , eps23, kplusp, msglvl, nconv ,
  226. & nevbef, nev0 , np0 , numcnv
  227. c
  228. c %-----------------------%
  229. c | Local array arguments |
  230. c %-----------------------%
  231. c
  232. integer kp(4)
  233. c
  234. c %----------------------%
  235. c | External Subroutines |
  236. c %----------------------%
  237. c
  238. external dcopy , dgetv0 , dnaitr , dnconv , dneigh ,
  239. & dngets , dnapps , dvout , ivout , arscnd
  240. c
  241. c %--------------------%
  242. c | External Functions |
  243. c %--------------------%
  244. c
  245. REAL*8
  246. external ddot , dnrm2 , dlapy2 , dlamch
  247. c
  248. c %---------------------%
  249. **c | Intrinsic Functions |
  250. **c %---------------------%
  251. **c
  252. ** intrinsic min, max, abs, sqrt
  253. **c
  254. **c %-----------------------%
  255. **c | Executable Statements |
  256. c %-----------------------%
  257. c
  258. if (ido .eq. 0) then
  259. c
  260. * call arscnd (t0)
  261. c
  262. msglvl = mnaup2
  263. c
  264. c %-------------------------------------%
  265. c | Get the machine dependent constant. |
  266. c %-------------------------------------%
  267. c
  268. eps23 = dlamch ('Epsilon-Machine')
  269. eps23 = eps23**(2.0D+0 / 3.0D+0 )
  270. c
  271. nev0 = nev
  272. np0 = np
  273. c
  274. c %-------------------------------------%
  275. c | kplusp is the bound on the largest |
  276. c | Lanczos factorization built. |
  277. c | nconv is the current number of |
  278. c | "converged" eigenvlues. |
  279. c | iter is the counter on the current |
  280. c | iteration step. |
  281. c %-------------------------------------%
  282. c
  283. kplusp = nev + np
  284. nconv = 0
  285. iter = 0
  286. c
  287. c %---------------------------------------%
  288. c | Set flags for computing the first NEV |
  289. c | steps of the Arnoldi factorization. |
  290. c %---------------------------------------%
  291. c
  292. getv0 = .true.
  293. update = .false.
  294. ushift = .false.
  295. cnorm = .false.
  296. c
  297. if (info .ne. 0) then
  298. c
  299. c %--------------------------------------------%
  300. c | User provides the initial residual vector. |
  301. c %--------------------------------------------%
  302. c
  303. initv = .true.
  304. info = 0
  305. else
  306. initv = .false.
  307. end if
  308. end if
  309. c
  310. c %---------------------------------------------%
  311. c | Get a possibly random starting vector and |
  312. c | force it into the range of the operator OP. |
  313. c %---------------------------------------------%
  314. c
  315. 10 continue
  316. c
  317. if (getv0) then
  318. call dgetv0 (ido, bmat, 1, initv, n, 1, v, ldv, resid, rnorm,
  319. & ipntr, workd, info)
  320. c
  321. if (ido .ne. 99) go to 9000
  322. c
  323. if (rnorm .eq. zero) then
  324. c
  325. c %-----------------------------------------%
  326. c | The initial vector is zero. Error exit. |
  327. c %-----------------------------------------%
  328. c
  329. info = -9
  330. go to 1100
  331. end if
  332. getv0 = .false.
  333. ido = 0
  334. end if
  335. c
  336. c %-----------------------------------%
  337. c | Back from reverse communication : |
  338. c | continue with update step |
  339. c %-----------------------------------%
  340. c
  341. if (update) go to 20
  342. c
  343. c %-------------------------------------------%
  344. c | Back from computing user specified shifts |
  345. c %-------------------------------------------%
  346. c
  347. if (ushift) go to 50
  348. c
  349. c %-------------------------------------%
  350. c | Back from computing residual norm |
  351. c | at the end of the current iteration |
  352. c %-------------------------------------%
  353. c
  354. if (cnorm) go to 100
  355. c
  356. c %----------------------------------------------------------%
  357. c | Compute the first NEV steps of the Arnoldi factorization |
  358. c %----------------------------------------------------------%
  359. c
  360. call dnaitr (ido, bmat, n, 0, nev, mode, resid, rnorm, v, ldv,
  361. & h, ldh, ipntr, workd, info)
  362. c
  363. c %---------------------------------------------------%
  364. c | ido .ne. 99 implies use of reverse communication |
  365. c | to compute operations involving OP and possibly B |
  366. c %---------------------------------------------------%
  367. c
  368. if (ido .ne. 99) go to 9000
  369. c
  370. if (info .gt. 0) then
  371. np = info
  372. mxiter = iter
  373. info = -9999
  374. go to 1200
  375. end if
  376. c
  377. c %--------------------------------------------------------------%
  378. c | |
  379. c | M A I N ARNOLDI I T E R A T I O N L O O P |
  380. c | Each iteration implicitly restarts the Arnoldi |
  381. c | factorization in place. |
  382. c | |
  383. c %--------------------------------------------------------------%
  384. c
  385. 1000 continue
  386. c
  387. iter = iter + 1
  388. c
  389. c if (msglvl .gt. 0) then
  390. c call ivout (logfil, 1, iter, ndigit,
  391. c & '_naup2: **** Start of major iteration number ****')
  392. c end if
  393. c
  394. c %-----------------------------------------------------------%
  395. c | Compute NP additional steps of the Arnoldi factorization. |
  396. c | Adjust NP since NEV might have been updated by last call |
  397. c | to the shift application routine dnapps . |
  398. c %-----------------------------------------------------------%
  399. c
  400. np = kplusp - nev
  401. c
  402. c if (msglvl .gt. 1) then
  403. c call ivout (logfil, 1, nev, ndigit,
  404. c & '_naup2: The length of the current Arnoldi factorization')
  405. c call ivout (logfil, 1, np, ndigit,
  406. c & '_naup2: Extend the Arnoldi factorization by')
  407. c end if
  408. c
  409. c %-----------------------------------------------------------%
  410. c | Compute NP additional steps of the Arnoldi factorization. |
  411. c %-----------------------------------------------------------%
  412. c
  413. ido = 0
  414. 20 continue
  415. update = .true.
  416. c
  417. call dnaitr (ido , bmat, n , nev, np , mode , resid,
  418. & rnorm, v , ldv, h , ldh, ipntr, workd,
  419. & info)
  420. c
  421. c %---------------------------------------------------%
  422. c | ido .ne. 99 implies use of reverse communication |
  423. c | to compute operations involving OP and possibly B |
  424. c %---------------------------------------------------%
  425. c
  426. if (ido .ne. 99) go to 9000
  427. c
  428. if (info .gt. 0) then
  429. np = info
  430. mxiter = iter
  431. info = -9999
  432. go to 1200
  433. end if
  434. update = .false.
  435. c
  436. c if (msglvl .gt. 1) then
  437. c call dvout (logfil, 1, rnorm, ndigit,
  438. c & '_naup2: Corresponding B-norm of the residual')
  439. c end if
  440. c
  441. c %--------------------------------------------------------%
  442. c | Compute the eigenvalues and corresponding error bounds |
  443. c | of the current upper Hessenberg matrix. |
  444. c %--------------------------------------------------------%
  445. c
  446. call dneigh (rnorm, kplusp, h, ldh, ritzr, ritzi, bounds,
  447. & q, ldq, workl, ierr)
  448. c
  449. if (ierr .ne. 0) then
  450. info = -8
  451. go to 1200
  452. end if
  453. c
  454. c %----------------------------------------------------%
  455. c | Make a copy of eigenvalues and corresponding error |
  456. c | bounds obtained from dneigh . |
  457. c %----------------------------------------------------%
  458. c
  459. call dcopy (kplusp, ritzr, 1, workl(kplusp**2+1), 1)
  460. call dcopy (kplusp, ritzi, 1, workl(kplusp**2+kplusp+1), 1)
  461. call dcopy (kplusp, bounds, 1, workl(kplusp**2+2*kplusp+1), 1)
  462. c
  463. c %---------------------------------------------------%
  464. c | Select the wanted Ritz values and their bounds |
  465. c | to be used in the convergence test. |
  466. c | The wanted part of the spectrum and corresponding |
  467. c | error bounds are in the last NEV loc. of RITZR, |
  468. c | RITZI and BOUNDS respectively. The variables NEV |
  469. c | and NP may be updated if the NEV-th wanted Ritz |
  470. c | value has a non zero imaginary part. In this case |
  471. c | NEV is increased by one and NP decreased by one. |
  472. c | NOTE: The last two arguments of dngets are no |
  473. c | longer used as of version 2.1. |
  474. c %---------------------------------------------------%
  475. c
  476. nev = nev0
  477. np = np0
  478. numcnv = nev
  479. call dngets (ishift, which, nev, np, ritzr, ritzi,
  480. & bounds, workl, workl(np+1))
  481. if (nev .eq. nev0+1) numcnv = nev0+1
  482. c
  483. c %-------------------%
  484. c | Convergence test. |
  485. c %-------------------%
  486. c
  487. call dcopy (nev, bounds(np+1), 1, workl(2*np+1), 1)
  488. call dnconv (nev, ritzr(np+1), ritzi(np+1), workl(2*np+1),
  489. & tol, nconv)
  490. c
  491. if (msglvl .gt. 2) then
  492. kp(1) = nev
  493. kp(2) = np
  494. kp(3) = numcnv
  495. kp(4) = nconv
  496. c call ivout (logfil, 4, kp, ndigit,
  497. c & '_naup2: NEV, NP, NUMCNV, NCONV are')
  498. c call dvout (logfil, kplusp, ritzr, ndigit,
  499. c & '_naup2: Real part of the eigenvalues of H')
  500. c call dvout (logfil, kplusp, ritzi, ndigit,
  501. c & '_naup2: Imaginary part of the eigenvalues of H')
  502. c call dvout (logfil, kplusp, bounds, ndigit,
  503. c & '_naup2: Ritz estimates of the current NCV Ritz values')
  504. end if
  505. c
  506. c %---------------------------------------------------------%
  507. c | Count the number of unwanted Ritz values that have zero |
  508. c | Ritz estimates. If any Ritz estimates are equal to zero |
  509. c | then a leading block of H of order equal to at least |
  510. c | the number of Ritz values with zero Ritz estimates has |
  511. c | split off. None of these Ritz values may be removed by |
  512. c | shifting. Decrease NP the number of shifts to apply. If |
  513. c | no shifts may be applied, then prepare to exit |
  514. c %---------------------------------------------------------%
  515. c
  516. nptemp = np
  517. do 30 j=1, nptemp
  518. if (bounds(j) .eq. zero) then
  519. np = np - 1
  520. nev = nev + 1
  521. end if
  522. 30 continue
  523. c
  524. if ( (nconv .ge. numcnv) .or.
  525. & (iter .gt. mxiter) .or.
  526. & (np .eq. 0) ) then
  527.  
  528. c if (msglvl .gt. 4) then
  529. c call dvout (logfil, kplusp, workl(kplusp**2+1), ndigit,
  530. c & '_naup2: Real part of the eig computed by _neigh:')
  531. c call dvout (logfil, kplusp, workl(kplusp**2+kplusp+1),
  532. c & ndigit,
  533. c & '_naup2: Imag part of the eig computed by _neigh:')
  534. c call dvout (logfil, kplusp, workl(kplusp**2+kplusp*2+1),
  535. c & ndigit,
  536. c & '_naup2: Ritz eistmates computed by _neigh:')
  537. c end if
  538. c
  539. c %------------------------------------------------%
  540. c | Prepare to exit. Put the converged Ritz values |
  541. c | and corresponding bounds in RITZ(1:NCONV) and |
  542. c | BOUNDS(1:NCONV) respectively. Then sort. Be |
  543. c | careful when NCONV > NP |
  544. c %------------------------------------------------%
  545. c
  546. c %------------------------------------------%
  547. c | Use h( 3,1 ) as storage to communicate |
  548. c | rnorm to _neupd if needed |
  549. c %------------------------------------------%
  550.  
  551. h(3,1) = rnorm
  552. c
  553. c %----------------------------------------------%
  554. c | To be consistent with dngets , we first do a |
  555. c | pre-processing sort in order to keep complex |
  556. c | conjugate pairs together. This is similar |
  557. c | to the pre-processing sort used in dngets |
  558. c | except that the sort is done in the opposite |
  559. c | order. |
  560. c %----------------------------------------------%
  561. c
  562. if (which .eq. 'LM') wprime = 'SR'
  563. if (which .eq. 'SM') wprime = 'LR'
  564. if (which .eq. 'LR') wprime = 'SM'
  565. if (which .eq. 'SR') wprime = 'LM'
  566. if (which .eq. 'LI') wprime = 'SM'
  567. if (which .eq. 'SI') wprime = 'LM'
  568. c
  569. call dsortc (wprime, .true., kplusp, ritzr, ritzi, bounds)
  570. c
  571. c %----------------------------------------------%
  572. c | Now sort Ritz values so that converged Ritz |
  573. c | values appear within the first NEV locations |
  574. c | of ritzr, ritzi and bounds, and the most |
  575. c | desired one appears at the front. |
  576. c %----------------------------------------------%
  577. c
  578. if (which .eq. 'LM') wprime = 'SM'
  579. if (which .eq. 'SM') wprime = 'LM'
  580. if (which .eq. 'LR') wprime = 'SR'
  581. if (which .eq. 'SR') wprime = 'LR'
  582. if (which .eq. 'LI') wprime = 'SI'
  583. if (which .eq. 'SI') wprime = 'LI'
  584. c
  585. call dsortc (wprime, .true., kplusp, ritzr, ritzi, bounds)
  586. c
  587. c %--------------------------------------------------%
  588. c | Scale the Ritz estimate of each Ritz value |
  589. c | by 1 / max(eps23,magnitude of the Ritz value). |
  590. c %--------------------------------------------------%
  591. c
  592. do 35 j = 1, numcnv
  593. temp = max(eps23,dlapy2 (ritzr(j),
  594. & ritzi(j)))
  595. bounds(j) = bounds(j)/temp
  596. 35 continue
  597. c
  598. c %----------------------------------------------------%
  599. c | Sort the Ritz values according to the scaled Ritz |
  600. c | esitmates. This will push all the converged ones |
  601. c | towards the front of ritzr, ritzi, bounds |
  602. c | (in the case when NCONV < NEV.) |
  603. c %----------------------------------------------------%
  604. c
  605. wprime = 'LR'
  606. call dsortc (wprime, .true., numcnv, bounds, ritzr, ritzi)
  607. c
  608. c %----------------------------------------------%
  609. c | Scale the Ritz estimate back to its original |
  610. c | value. |
  611. c %----------------------------------------------%
  612. c
  613. do 40 j = 1, numcnv
  614. temp = max(eps23, dlapy2 (ritzr(j),
  615. & ritzi(j)))
  616. bounds(j) = bounds(j)*temp
  617. 40 continue
  618. c
  619. c %------------------------------------------------%
  620. c | Sort the converged Ritz values again so that |
  621. c | the "threshold" value appears at the front of |
  622. c | ritzr, ritzi and bound. |
  623. c %------------------------------------------------%
  624. c
  625. call dsortc (which, .true., nconv, ritzr, ritzi, bounds)
  626. c
  627. c if (msglvl .gt. 1) then
  628. c call dvout (logfil, kplusp, ritzr, ndigit,
  629. c & '_naup2: Sorted real part of the eigenvalues')
  630. c call dvout (logfil, kplusp, ritzi, ndigit,
  631. c & '_naup2: Sorted imaginary part of the eigenvalues')
  632. c call dvout (logfil, kplusp, bounds, ndigit,
  633. c & '_naup2: Sorted ritz estimates.')
  634. c end if
  635.  
  636. c %------------------------------------%
  637. c | Max iterations have been exceeded. |
  638. c %------------------------------------%
  639. c
  640. if (iter .gt. mxiter .and. nconv .lt. numcnv) info = 1
  641. c
  642. c %---------------------%
  643. c | No shifts to apply. |
  644. c %---------------------%
  645. c
  646. if (np .eq. 0 .and. nconv .lt. numcnv) info = 2
  647. c
  648. np = nconv
  649. go to 1100
  650. c
  651. else if ( (nconv .lt. numcnv) .and. (ishift .eq. 1) ) then
  652. c
  653. c %-------------------------------------------------%
  654. c | Do not have all the requested eigenvalues yet. |
  655. c | To prevent possible stagnation, adjust the size |
  656. c | of NEV. |
  657. c %-------------------------------------------------%
  658. c
  659. nevbef = nev
  660. nev = nev + min(nconv, np/2)
  661. if (nev .eq. 1 .and. kplusp .ge. 6) then
  662. nev = kplusp / 2
  663. else if (nev .eq. 1 .and. kplusp .gt. 3) then
  664. nev = 2
  665. end if
  666. c %---- Scipy fix ------------------------------------------------
  667. c | We must keep nev below this value, as otherwise we can get
  668. c | np == 0 (note that dngets below can bump nev by 1). If np == 0,
  669. c | the next call to `dnaitr` will write out-of-bounds.
  670. c |
  671. if (nev .gt. kplusp - 2) then
  672. nev = kplusp - 2
  673. end if
  674. c |
  675. c %---- Scipy fix end --------------------------------------------
  676. c
  677. np = kplusp - nev
  678. c
  679. c %---------------------------------------%
  680. c | If the size of NEV was just increased |
  681. c | resort the eigenvalues. |
  682. c %---------------------------------------%
  683. c
  684. if (nevbef .lt. nev)
  685. & call dngets (ishift, which, nev, np, ritzr, ritzi,
  686. & bounds, workl, workl(np+1))
  687. c
  688. end if
  689. c
  690. c if (msglvl .gt. 0) then
  691. c call ivout (logfil, 1, nconv, ndigit,
  692. c & '_naup2: no. of "converged" Ritz values at this iter.')
  693. c if (msglvl .gt. 1) then
  694. c kp(1) = nev
  695. c kp(2) = np
  696. c call ivout (logfil, 2, kp, ndigit,
  697. c & '_naup2: NEV and NP are')
  698. c call dvout (logfil, nev, ritzr(np+1), ndigit,
  699. c & '_naup2: "wanted" Ritz values -- real part')
  700. c call dvout (logfil, nev, ritzi(np+1), ndigit,
  701. c & '_naup2: "wanted" Ritz values -- imag part')
  702. c call dvout (logfil, nev, bounds(np+1), ndigit,
  703. c & '_naup2: Ritz estimates of the "wanted" values ')
  704. c end if
  705. c end if
  706. c
  707. if (ishift .eq. 0) then
  708. c
  709. c %-------------------------------------------------------%
  710. c | User specified shifts: reverse comminucation to |
  711. c | compute the shifts. They are returned in the first |
  712. c | 2*NP locations of WORKL. |
  713. c %-------------------------------------------------------%
  714. c
  715. ushift = .true.
  716. ido = 3
  717. go to 9000
  718. end if
  719. c
  720. 50 continue
  721. c
  722. c %------------------------------------%
  723. c | Back from reverse communication; |
  724. c | User specified shifts are returned |
  725. c | in WORKL(1:2*NP) |
  726. c %------------------------------------%
  727. c
  728. ushift = .false.
  729. c
  730. if ( ishift .eq. 0 ) then
  731. c
  732. c %----------------------------------%
  733. c | Move the NP shifts from WORKL to |
  734. c | RITZR, RITZI to free up WORKL |
  735. c | for non-exact shift case. |
  736. c %----------------------------------%
  737. c
  738. call dcopy (np, workl, 1, ritzr, 1)
  739. call dcopy (np, workl(np+1), 1, ritzi, 1)
  740. end if
  741. c
  742. c if (msglvl .gt. 2) then
  743. c call ivout (logfil, 1, np, ndigit,
  744. c & '_naup2: The number of shifts to apply ')
  745. c call dvout (logfil, np, ritzr, ndigit,
  746. c & '_naup2: Real part of the shifts')
  747. c call dvout (logfil, np, ritzi, ndigit,
  748. c & '_naup2: Imaginary part of the shifts')
  749. c if ( ishift .eq. 1 )
  750. c & call dvout (logfil, np, bounds, ndigit,
  751. c & '_naup2: Ritz estimates of the shifts')
  752. c end if
  753. c
  754. c %---------------------------------------------------------%
  755. c | Apply the NP implicit shifts by QR bulge chasing. |
  756. c | Each shift is applied to the whole upper Hessenberg |
  757. c | matrix H. |
  758. c | The first 2*N locations of WORKD are used as workspace. |
  759. c %---------------------------------------------------------%
  760. c
  761. call dnapps (n, nev, np, ritzr, ritzi, v, ldv,
  762. & h, ldh, resid, q, ldq, workl, workd)
  763. c
  764. c %---------------------------------------------%
  765. c | Compute the B-norm of the updated residual. |
  766. c | Keep B*RESID in WORKD(1:N) to be used in |
  767. c | the first step of the next call to dnaitr . |
  768. c %---------------------------------------------%
  769. c
  770. cnorm = .true.
  771. * call arscnd (t2)
  772. if (bmat .eq. 'G') then
  773. nbx = nbx + 1
  774. call dcopy (n, resid, 1, workd(n+1), 1)
  775. ipntr(1) = n + 1
  776. ipntr(2) = 1
  777. ido = 2
  778. c
  779. c %----------------------------------%
  780. c | Exit in order to compute B*RESID |
  781. c %----------------------------------%
  782. c
  783. go to 9000
  784. else if (bmat .eq. 'I') then
  785. call dcopy (n, resid, 1, workd, 1)
  786. end if
  787. c
  788. 100 continue
  789. c
  790. c %----------------------------------%
  791. c | Back from reverse communication; |
  792. c | WORKD(1:N) := B*RESID |
  793. c %----------------------------------%
  794. c
  795. if (bmat .eq. 'G') then
  796. * call arscnd (t3)
  797. tmvbx = tmvbx + (t3 - t2)
  798. end if
  799. c
  800. if (bmat .eq. 'G') then
  801. rnorm = ddot (n, resid, 1, workd, 1)
  802. rnorm = sqrt(abs(rnorm))
  803. else if (bmat .eq. 'I') then
  804. rnorm = dnrm2 (n, resid, 1)
  805. end if
  806. cnorm = .false.
  807. c
  808. if (msglvl .gt. 2) then
  809. c call dvout (logfil, 1, rnorm, ndigit,
  810. c & '_naup2: B-norm of residual for compressed factorization')
  811. call dmout (logfil, nev, nev, h, ldh, ndigit,
  812. & '_naup2: Compressed upper Hessenberg matrix H')
  813. end if
  814. c
  815. go to 1000
  816. c
  817. c %---------------------------------------------------------------%
  818. c | |
  819. c | E N D O F M A I N I T E R A T I O N L O O P |
  820. c | |
  821. c %---------------------------------------------------------------%
  822. c
  823. 1100 continue
  824. c
  825. mxiter = iter
  826. nev = numcnv
  827. c
  828. 1200 continue
  829. ido = 99
  830. c
  831. c %------------%
  832. c | Error Exit |
  833. c %------------%
  834. c
  835. * call arscnd (t1)
  836. tnaup2 = t1 - t0
  837. c
  838. 9000 continue
  839. c
  840. c %---------------%
  841. c | End of dnaup2 |
  842. c %---------------%
  843. c
  844. return
  845. end
  846.  
  847.  

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