Télécharger dnaup2.eso

Retour à la liste

Numérotation des lignes :

  1. C DNAUP2 SOURCE BP208322 20/02/06 21:15:20 10512
  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. -> deleted by BP in 2020
  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, ITRAK, info )
  179. c
  180. c %----------------------------------------------------%
  181. c | Include files for debugging and timing information |
  182. c -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. INTEGER ITRAK(5)
  202. REAL*8
  203. & bounds(nev+np), h(ldh,nev+np), q(ldq,nev+np), resid(n),
  204. & ritzi(nev+np), ritzr(nev+np), v(ldv,nev+np),
  205. & workd(3*n), workl( (nev+np)*(nev+np+3) )
  206. c
  207. c %------------%
  208. c | Parameters |
  209. c %------------%
  210. c
  211. REAL*8
  212. & one, zero
  213. parameter (one = 1.0D+0 , zero = 0.0D+0 )
  214. c
  215. c %---------------%
  216. c | Local Scalars |
  217. c %---------------%
  218. c
  219. character wprime*2
  220. logical cnorm , getv0, initv, update, ushift
  221. integer ierr , iter , j , kplusp, msglvl, nconv,
  222. & nevbef, nev0 , np0 , nptemp, numcnv
  223. REAL*8
  224. & rnorm , temp , eps23
  225. save cnorm , getv0, initv, update, ushift,
  226. c & rnorm , iter , eps23, kplusp, msglvl, nconv ,
  227. & rnorm , iter , eps23, kplusp, nconv ,
  228. & nevbef, nev0 , np0 , numcnv
  229. parameter (msglvl=0)
  230. c
  231. c c %-----------------------%
  232. c c | Local array arguments |
  233. c c %-----------------------%
  234. c c
  235. c integer kp(4)
  236. c
  237. c %----------------------%
  238. c | External Subroutines |
  239. c %----------------------%
  240. c
  241. external dcopy , dgetv0 , dnaitr , dnconv , dneigh ,
  242. c
  243. c %--------------------%
  244. c | External Functions |
  245. c %--------------------%
  246. c
  247. REAL*8
  248. external ddot , dnrm2 , dlapy2 , dlamch
  249. c
  250. c %---------------------%
  251. **c | Intrinsic Functions |
  252. **c %---------------------%
  253. **c
  254. ** intrinsic min, max, abs, sqrt
  255. **c
  256. **c %-----------------------%
  257. **c | Executable Statements |
  258. c %-----------------------%
  259. NOPX =ITRAK(1)
  260. NBX =ITRAK(2)
  261. NRORTH=ITRAK(3)
  262. NITREF=ITRAK(4)
  263. NRSTRT=ITRAK(5)
  264. c
  265. if (ido .eq. 0) then
  266. c
  267. * call arscnd (t0)
  268. c msglvl = mnaup2
  269. c
  270. c %-------------------------------------%
  271. c | Get the machine dependent constant. |
  272. c %-------------------------------------%
  273. c
  274. eps23 = dlamch ('Epsilon-Machine')
  275. eps23 = eps23**(2.0D+0 / 3.0D+0 )
  276. c
  277. nev0 = nev
  278. np0 = np
  279. c
  280. c %-------------------------------------%
  281. c | kplusp is the bound on the largest |
  282. c | Lanczos factorization built. |
  283. c | nconv is the current number of |
  284. c | "converged" eigenvlues. |
  285. c | iter is the counter on the current |
  286. c | iteration step. |
  287. c %-------------------------------------%
  288. c
  289. kplusp = nev + np
  290. nconv = 0
  291. iter = 0
  292. c
  293. c %---------------------------------------%
  294. c | Set flags for computing the first NEV |
  295. c | steps of the Arnoldi factorization. |
  296. c %---------------------------------------%
  297. c
  298. getv0 = .true.
  299. update = .false.
  300. ushift = .false.
  301. cnorm = .false.
  302. c
  303. if (info .ne. 0) then
  304. c
  305. c %--------------------------------------------%
  306. c | User provides the initial residual vector. |
  307. c %--------------------------------------------%
  308. c
  309. initv = .true.
  310. info = 0
  311. else
  312. initv = .false.
  313. end if
  314. end if
  315. c
  316. c %---------------------------------------------%
  317. c | Get a possibly random starting vector and |
  318. c | force it into the range of the operator OP. |
  319. c %---------------------------------------------%
  320. c
  321. 10 continue
  322. c
  323. if (getv0) then
  324. call dgetv0 (ido, bmat, 1, initv, n, 1, v, ldv, resid, rnorm,
  325. & ipntr, workd, nopx,nbx, info)
  326. c
  327. if (ido .ne. 99) go to 9000
  328. c
  329. if (rnorm .eq. zero) then
  330. c
  331. c %-----------------------------------------%
  332. c | The initial vector is zero. Error exit. |
  333. c %-----------------------------------------%
  334. c
  335. info = -9
  336. go to 1100
  337. end if
  338. getv0 = .false.
  339. ido = 0
  340. end if
  341. c
  342. c %-----------------------------------%
  343. c | Back from reverse communication : |
  344. c | continue with update step |
  345. c %-----------------------------------%
  346. c
  347. if (update) go to 20
  348. c
  349. c %-------------------------------------------%
  350. c | Back from computing user specified shifts |
  351. c %-------------------------------------------%
  352. c
  353. if (ushift) go to 50
  354. c
  355. c %-------------------------------------%
  356. c | Back from computing residual norm |
  357. c | at the end of the current iteration |
  358. c %-------------------------------------%
  359. c
  360. if (cnorm) go to 100
  361. c
  362. c %----------------------------------------------------------%
  363. c | Compute the first NEV steps of the Arnoldi factorization |
  364. c %----------------------------------------------------------%
  365. c
  366. call dnaitr (ido, bmat, n, 0, nev, mode, resid, rnorm, v, ldv,
  367. & h, ldh, ipntr, workd, ITRAK, info)
  368. c
  369. c %---------------------------------------------------%
  370. c | ido .ne. 99 implies use of reverse communication |
  371. c | to compute operations involving OP and possibly B |
  372. c %---------------------------------------------------%
  373. c
  374. if (ido .ne. 99) go to 9000
  375. c
  376. if (info .gt. 0) then
  377. np = info
  378. mxiter = iter
  379. info = -9999
  380. go to 1200
  381. end if
  382. c
  383. c %--------------------------------------------------------------%
  384. c | |
  385. c | M A I N ARNOLDI I T E R A T I O N L O O P |
  386. c | Each iteration implicitly restarts the Arnoldi |
  387. c | factorization in place. |
  388. c | |
  389. c %--------------------------------------------------------------%
  390. c
  391. 1000 continue
  392. c
  393. iter = iter + 1
  394.  
  395. if (msglvl .gt. 0) then
  396. call ivout ( 1, iter, ndigit,
  397. & '_naup2: **** Start of major iteration number ****')
  398. end if
  399.  
  400. c %-----------------------------------------------------------%
  401. c | Compute NP additional steps of the Arnoldi factorization. |
  402. c | Adjust NP since NEV might have been updated by last call |
  403. c | to the shift application routine dnapps . |
  404. c %-----------------------------------------------------------%
  405. c
  406. np = kplusp - nev
  407. c
  408. c if (msglvl .gt. 1) then
  409. c call ivout ( 1, nev, ndigit,
  410. c & '_naup2: The length of the current Arnoldi factorization')
  411. c call ivout ( 1, np, ndigit,
  412. c & '_naup2: Extend the Arnoldi factorization by')
  413. c end if
  414. c
  415. c %-----------------------------------------------------------%
  416. c | Compute NP additional steps of the Arnoldi factorization. |
  417. c %-----------------------------------------------------------%
  418. c
  419. ido = 0
  420. 20 continue
  421. update = .true.
  422. c
  423. call dnaitr (ido , bmat, n , nev, np , mode , resid,
  424. & rnorm, v , ldv, h , ldh, ipntr, workd,
  425. & ITRAK, info)
  426. c
  427. c %---------------------------------------------------%
  428. c | ido .ne. 99 implies use of reverse communication |
  429. c | to compute operations involving OP and possibly B |
  430. c %---------------------------------------------------%
  431. c
  432. if (ido .ne. 99) go to 9000
  433. c
  434. if (info .gt. 0) then
  435. np = info
  436. mxiter = iter
  437. info = -9999
  438. go to 1200
  439. end if
  440. update = .false.
  441. c
  442. c if (msglvl .gt. 1) then
  443. c call dvout ( 1, rnorm, ndigit,
  444. c & '_naup2: Corresponding B-norm of the residual')
  445. c end if
  446. c
  447. c %--------------------------------------------------------%
  448. c | Compute the eigenvalues and corresponding error bounds |
  449. c | of the current upper Hessenberg matrix. |
  450. c %--------------------------------------------------------%
  451. c
  452. call dneigh (rnorm, kplusp, h, ldh, ritzr, ritzi, bounds,
  453. & q, ldq, workl, ierr)
  454. c
  455. if (ierr .ne. 0) then
  456. info = -8
  457. go to 1200
  458. end if
  459. c
  460. c %----------------------------------------------------%
  461. c | Make a copy of eigenvalues and corresponding error |
  462. c | bounds obtained from dneigh . |
  463. c %----------------------------------------------------%
  464. c
  465. call dcopy (kplusp, ritzr, 1, workl(kplusp**2+1), 1)
  466. call dcopy (kplusp, ritzi, 1, workl(kplusp**2+kplusp+1), 1)
  467. call dcopy (kplusp, bounds, 1, workl(kplusp**2+2*kplusp+1), 1)
  468. c
  469. c %---------------------------------------------------%
  470. c | Select the wanted Ritz values and their bounds |
  471. c | to be used in the convergence test. |
  472. c | The wanted part of the spectrum and corresponding |
  473. c | error bounds are in the last NEV loc. of RITZR, |
  474. c | RITZI and BOUNDS respectively. The variables NEV |
  475. c | and NP may be updated if the NEV-th wanted Ritz |
  476. c | value has a non zero imaginary part. In this case |
  477. c | NEV is increased by one and NP decreased by one. |
  478. c | NOTE: The last two arguments of dngets are no |
  479. c | longer used as of version 2.1. |
  480. c %---------------------------------------------------%
  481. c
  482. nev = nev0
  483. np = np0
  484. numcnv = nev
  485. call dngets (ishift, which, nev, np, ritzr, ritzi,
  486. & bounds, workl, workl(np+1))
  487. if (nev .eq. nev0+1) numcnv = nev0+1
  488. c
  489. c %-------------------%
  490. c | Convergence test. |
  491. c %-------------------%
  492. c
  493. call dcopy (nev, bounds(np+1), 1, workl(2*np+1), 1)
  494. call dnconv (nev, ritzr(np+1), ritzi(np+1), workl(2*np+1),
  495. & tol, nconv)
  496. c
  497. c if (msglvl .gt. 2) then
  498. c kp(1) = nev
  499. c kp(2) = np
  500. c kp(3) = numcnv
  501. c kp(4) = nconv
  502. c call ivout ( 4, kp, ndigit,
  503. c & '_naup2: NEV, NP, NUMCNV, NCONV are')
  504. c call dvout ( kplusp, ritzr, ndigit,
  505. c & '_naup2: Real part of the eigenvalues of H')
  506. c call dvout ( kplusp, ritzi, ndigit,
  507. c & '_naup2: Imaginary part of the eigenvalues of H')
  508. c call dvout ( kplusp, bounds, ndigit,
  509. c & '_naup2: Ritz estimates of the current NCV Ritz values')
  510. c end if
  511. c
  512. c %---------------------------------------------------------%
  513. c | Count the number of unwanted Ritz values that have zero |
  514. c | Ritz estimates. If any Ritz estimates are equal to zero |
  515. c | then a leading block of H of order equal to at least |
  516. c | the number of Ritz values with zero Ritz estimates has |
  517. c | split off. None of these Ritz values may be removed by |
  518. c | shifting. Decrease NP the number of shifts to apply. If |
  519. c | no shifts may be applied, then prepare to exit |
  520. c %---------------------------------------------------------%
  521. c
  522. nptemp = np
  523. do 30 j=1, nptemp
  524. if (bounds(j) .eq. zero) then
  525. np = np - 1
  526. nev = nev + 1
  527. end if
  528. 30 continue
  529. c
  530. if ( (nconv .ge. numcnv) .or.
  531. & (iter .gt. mxiter) .or.
  532. & (np .eq. 0) ) then
  533.  
  534. c if (msglvl .gt. 4) then
  535. c call dvout ( kplusp, workl(kplusp**2+1), ndigit,
  536. c & '_naup2: Real part of the eig computed by _neigh:')
  537. c call dvout ( kplusp, workl(kplusp**2+kplusp+1),
  538. c & ndigit,
  539. c & '_naup2: Imag part of the eig computed by _neigh:')
  540. c call dvout ( kplusp, workl(kplusp**2+kplusp*2+1),
  541. c & ndigit,
  542. c & '_naup2: Ritz eistmates computed by _neigh:')
  543. c end if
  544. c
  545. c %------------------------------------------------%
  546. c | Prepare to exit. Put the converged Ritz values |
  547. c | and corresponding bounds in RITZ(1:NCONV) and |
  548. c | BOUNDS(1:NCONV) respectively. Then sort. Be |
  549. c | careful when NCONV > NP |
  550. c %------------------------------------------------%
  551. c
  552. c %------------------------------------------%
  553. c | Use h( 3,1 ) as storage to communicate |
  554. c | rnorm to _neupd if needed |
  555. c %------------------------------------------%
  556.  
  557. h(3,1) = rnorm
  558. c
  559. c %----------------------------------------------%
  560. c | To be consistent with dngets , we first do a |
  561. c | pre-processing sort in order to keep complex |
  562. c | conjugate pairs together. This is similar |
  563. c | to the pre-processing sort used in dngets |
  564. c | except that the sort is done in the opposite |
  565. c | order. |
  566. c %----------------------------------------------%
  567. c
  568. if (which .eq. 'LM') wprime = 'SR'
  569. if (which .eq. 'SM') wprime = 'LR'
  570. if (which .eq. 'LR') wprime = 'SM'
  571. if (which .eq. 'SR') wprime = 'LM'
  572. if (which .eq. 'LI') wprime = 'SM'
  573. if (which .eq. 'SI') wprime = 'LM'
  574. c
  575. call dsortc (wprime, .true., kplusp, ritzr, ritzi, bounds)
  576. c
  577. c %----------------------------------------------%
  578. c | Now sort Ritz values so that converged Ritz |
  579. c | values appear within the first NEV locations |
  580. c | of ritzr, ritzi and bounds, and the most |
  581. c | desired one appears at the front. |
  582. c %----------------------------------------------%
  583. c
  584. if (which .eq. 'LM') wprime = 'SM'
  585. if (which .eq. 'SM') wprime = 'LM'
  586. if (which .eq. 'LR') wprime = 'SR'
  587. if (which .eq. 'SR') wprime = 'LR'
  588. if (which .eq. 'LI') wprime = 'SI'
  589. if (which .eq. 'SI') wprime = 'LI'
  590. c
  591. call dsortc (wprime, .true., kplusp, ritzr, ritzi, bounds)
  592. c
  593. c %--------------------------------------------------%
  594. c | Scale the Ritz estimate of each Ritz value |
  595. c | by 1 / max(eps23,magnitude of the Ritz value). |
  596. c %--------------------------------------------------%
  597. c
  598. do 35 j = 1, numcnv
  599. temp = max(eps23,dlapy2 (ritzr(j),
  600. & ritzi(j)))
  601. bounds(j) = bounds(j)/temp
  602. 35 continue
  603. c
  604. c %----------------------------------------------------%
  605. c | Sort the Ritz values according to the scaled Ritz |
  606. c | esitmates. This will push all the converged ones |
  607. c | towards the front of ritzr, ritzi, bounds |
  608. c | (in the case when NCONV < NEV.) |
  609. c %----------------------------------------------------%
  610. c
  611. wprime = 'LR'
  612. call dsortc (wprime, .true., numcnv, bounds, ritzr, ritzi)
  613. c
  614. c %----------------------------------------------%
  615. c | Scale the Ritz estimate back to its original |
  616. c | value. |
  617. c %----------------------------------------------%
  618. c
  619. do 40 j = 1, numcnv
  620. temp = max(eps23, dlapy2 (ritzr(j),
  621. & ritzi(j)))
  622. bounds(j) = bounds(j)*temp
  623. 40 continue
  624. c
  625. c %------------------------------------------------%
  626. c | Sort the converged Ritz values again so that |
  627. c | the "threshold" value appears at the front of |
  628. c | ritzr, ritzi and bound. |
  629. c %------------------------------------------------%
  630. c
  631. call dsortc (which, .true., nconv, ritzr, ritzi, bounds)
  632. c
  633. c if (msglvl .gt. 1) then
  634. c call dvout ( kplusp, ritzr, ndigit,
  635. c & '_naup2: Sorted real part of the eigenvalues')
  636. c call dvout ( kplusp, ritzi, ndigit,
  637. c & '_naup2: Sorted imaginary part of the eigenvalues')
  638. c call dvout ( kplusp, bounds, ndigit,
  639. c & '_naup2: Sorted ritz estimates.')
  640. c end if
  641.  
  642. c %------------------------------------%
  643. c | Max iterations have been exceeded. |
  644. c %------------------------------------%
  645. c
  646. if (iter .gt. mxiter .and. nconv .lt. numcnv) info = 1
  647. c
  648. c %---------------------%
  649. c | No shifts to apply. |
  650. c %---------------------%
  651. c
  652. if (np .eq. 0 .and. nconv .lt. numcnv) info = 2
  653. c
  654. np = nconv
  655. go to 1100
  656. c
  657. else if ( (nconv .lt. numcnv) .and. (ishift .eq. 1) ) then
  658. c
  659. c %-------------------------------------------------%
  660. c | Do not have all the requested eigenvalues yet. |
  661. c | To prevent possible stagnation, adjust the size |
  662. c | of NEV. |
  663. c %-------------------------------------------------%
  664. c
  665. nevbef = nev
  666. nev = nev + min(nconv, np/2)
  667. if (nev .eq. 1 .and. kplusp .ge. 6) then
  668. nev = kplusp / 2
  669. else if (nev .eq. 1 .and. kplusp .gt. 3) then
  670. nev = 2
  671. end if
  672. c %---- Scipy fix ------------------------------------------------
  673. c | We must keep nev below this value, as otherwise we can get
  674. c | np == 0 (note that dngets below can bump nev by 1). If np == 0,
  675. c | the next call to `dnaitr` will write out-of-bounds.
  676. c |
  677. if (nev .gt. kplusp - 2) then
  678. nev = kplusp - 2
  679. end if
  680. c |
  681. c %---- Scipy fix end --------------------------------------------
  682. c
  683. np = kplusp - nev
  684. c
  685. c %---------------------------------------%
  686. c | If the size of NEV was just increased |
  687. c | resort the eigenvalues. |
  688. c %---------------------------------------%
  689. c
  690. if (nevbef .lt. nev)
  691. & call dngets (ishift, which, nev, np, ritzr, ritzi,
  692. & bounds, workl, workl(np+1))
  693. c
  694. end if
  695. c
  696. if (msglvl .gt. 0) then
  697. call ivout ( 1, nconv, ndigit,
  698. & '_naup2: no. of "converged" Ritz values at this iter.')
  699. c if (msglvl .gt. 1) then
  700. c kp(1) = nev
  701. c kp(2) = np
  702. c call ivout ( 2, kp, ndigit,
  703. c & '_naup2: NEV and NP are')
  704. c call dvout ( nev, ritzr(np+1), ndigit,
  705. c & '_naup2: "wanted" Ritz values -- real part')
  706. c call dvout ( nev, ritzi(np+1), ndigit,
  707. c & '_naup2: "wanted" Ritz values -- imag part')
  708. c call dvout ( nev, bounds(np+1), ndigit,
  709. c & '_naup2: Ritz estimates of the "wanted" values ')
  710. c end if
  711. end if
  712. c
  713. if (ishift .eq. 0) then
  714. c
  715. c %-------------------------------------------------------%
  716. c | User specified shifts: reverse comminucation to |
  717. c | compute the shifts. They are returned in the first |
  718. c | 2*NP locations of WORKL. |
  719. c %-------------------------------------------------------%
  720. c
  721. ushift = .true.
  722. ido = 3
  723. go to 9000
  724. end if
  725. c
  726. 50 continue
  727. c
  728. c %------------------------------------%
  729. c | Back from reverse communication; |
  730. c | User specified shifts are returned |
  731. c | in WORKL(1:2*NP) |
  732. c %------------------------------------%
  733. c
  734. ushift = .false.
  735. c
  736. if ( ishift .eq. 0 ) then
  737. c
  738. c %----------------------------------%
  739. c | Move the NP shifts from WORKL to |
  740. c | RITZR, RITZI to free up WORKL |
  741. c | for non-exact shift case. |
  742. c %----------------------------------%
  743. c
  744. call dcopy (np, workl, 1, ritzr, 1)
  745. call dcopy (np, workl(np+1), 1, ritzi, 1)
  746. end if
  747. c
  748. c if (msglvl .gt. 2) then
  749. c call ivout ( 1, np, ndigit,
  750. c & '_naup2: The number of shifts to apply ')
  751. c call dvout ( np, ritzr, ndigit,
  752. c & '_naup2: Real part of the shifts')
  753. c call dvout ( np, ritzi, ndigit,
  754. c & '_naup2: Imaginary part of the shifts')
  755. c if ( ishift .eq. 1 )
  756. c & call dvout ( np, bounds, ndigit,
  757. c & '_naup2: Ritz estimates of the shifts')
  758. c end if
  759. c
  760. c %---------------------------------------------------------%
  761. c | Apply the NP implicit shifts by QR bulge chasing. |
  762. c | Each shift is applied to the whole upper Hessenberg |
  763. c | matrix H. |
  764. c | The first 2*N locations of WORKD are used as workspace. |
  765. c %---------------------------------------------------------%
  766. c
  767. call dnapps (n, nev, np, ritzr, ritzi, v, ldv,
  768. & h, ldh, resid, q, ldq, workl, workd)
  769. c
  770. c %---------------------------------------------%
  771. c | Compute the B-norm of the updated residual. |
  772. c | Keep B*RESID in WORKD(1:N) to be used in |
  773. c | the first step of the next call to dnaitr . |
  774. c %---------------------------------------------%
  775. c
  776. cnorm = .true.
  777. * call arscnd (t2)
  778. if (bmat .eq. 'G') then
  779. nbx = nbx + 1
  780. call dcopy (n, resid, 1, workd(n+1), 1)
  781. ipntr(1) = n + 1
  782. ipntr(2) = 1
  783. ido = 2
  784. c
  785. c %----------------------------------%
  786. c | Exit in order to compute B*RESID |
  787. c %----------------------------------%
  788. c
  789. go to 9000
  790. else if (bmat .eq. 'I') then
  791. call dcopy (n, resid, 1, workd, 1)
  792. end if
  793. c
  794. 100 continue
  795. c
  796. c %----------------------------------%
  797. c | Back from reverse communication; |
  798. c | WORKD(1:N) := B*RESID |
  799. c %----------------------------------%
  800. c
  801. c if (bmat .eq. 'G') then
  802. * call arscnd (t3)
  803. c tmvbx = tmvbx + (t3 - t2)
  804. c end if
  805. c
  806. if (bmat .eq. 'G') then
  807. rnorm = ddot (n, resid, 1, workd, 1)
  808. rnorm = sqrt(abs(rnorm))
  809. else if (bmat .eq. 'I') then
  810. rnorm = dnrm2 (n, resid, 1)
  811. end if
  812. cnorm = .false.
  813. c
  814. c if (msglvl .gt. 2) then
  815. c c call dvout ( 1, rnorm, ndigit,
  816. c c & '_naup2: B-norm of residual for compressed factorization')
  817. c call dmout (logfil, nev, nev, h, ldh, ndigit,
  818. c & '_naup2: Compressed upper Hessenberg matrix H')
  819. c end if
  820. c
  821. go to 1000
  822. c
  823. c %---------------------------------------------------------------%
  824. c | |
  825. c | E N D O F M A I N I T E R A T I O N L O O P |
  826. c | |
  827. c %---------------------------------------------------------------%
  828. c
  829. 1100 continue
  830. c
  831. mxiter = iter
  832. nev = numcnv
  833. c
  834. 1200 continue
  835. ido = 99
  836. c
  837. c %------------%
  838. c | Error Exit |
  839. c %------------%
  840. c
  841. * call arscnd (t1)
  842. c tnaup2 = t1 - t0
  843. c
  844. 9000 continue
  845. ITRAK(1)=NOPX
  846. ITRAK(2)=NBX
  847. ITRAK(3)=NRORTH
  848. ITRAK(4)=NITREF
  849. ITRAK(5)=NRSTRT
  850. c
  851. c %---------------%
  852. c | End of dnaup2 |
  853. c %---------------%
  854. c
  855. end
  856.  
  857.  
  858.  

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