Télécharger dngets.eso

Retour à la liste

Numérotation des lignes :

  1. C DNGETS SOURCE BP208322 15/10/13 21:15:47 8670
  2. c-----------------------------------------------------------------------
  3. c\BeginDoc
  4. c
  5. c\Name: dngets
  6. c
  7. c\Description:
  8. c Given the eigenvalues of the upper Hessenberg matrix H,
  9. c computes the NP shifts AMU that are zeros of the polynomial of
  10. c degree NP which filters out components of the unwanted eigenvectors
  11. c corresponding to the AMU's based on some given criteria.
  12. c
  13. c NOTE: call this even in the case of user specified shifts in order
  14. c to sort the eigenvalues, and error bounds of H for later use.
  15. c
  16. c\Usage:
  17. c call dngets
  18. c ( ISHIFT, WHICH, KEV, NP, RITZR, RITZI, BOUNDS, SHIFTR, SHIFTI )
  19. c
  20. c\Arguments
  21. c ISHIFT Integer. (INPUT)
  22. c Method for selecting the implicit shifts at each iteration.
  23. c ISHIFT = 0: user specified shifts
  24. c ISHIFT = 1: exact shift with respect to the matrix H.
  25. c
  26. c WHICH Character*2. (INPUT)
  27. c Shift selection criteria.
  28. c 'LM' -> want the KEV eigenvalues of largest magnitude.
  29. c 'SM' -> want the KEV eigenvalues of smallest magnitude.
  30. c 'LR' -> want the KEV eigenvalues of largest real part.
  31. c 'SR' -> want the KEV eigenvalues of smallest real part.
  32. c 'LI' -> want the KEV eigenvalues of largest imaginary part.
  33. c 'SI' -> want the KEV eigenvalues of smallest imaginary part.
  34. c
  35. c KEV Integer. (INPUT/OUTPUT)
  36. c INPUT: KEV+NP is the size of the matrix H.
  37. c OUTPUT: Possibly increases KEV by one to keep complex conjugate
  38. c pairs together.
  39. c
  40. c NP Integer. (INPUT/OUTPUT)
  41. c Number of implicit shifts to be computed.
  42. c OUTPUT: Possibly decreases NP by one to keep complex conjugate
  43. c pairs together.
  44. c
  45. c RITZR, Double precision array of length KEV+NP. (INPUT/OUTPUT)
  46. c RITZI On INPUT, RITZR and RITZI contain the real and imaginary
  47. c parts of the eigenvalues of H.
  48. c On OUTPUT, RITZR and RITZI are sorted so that the unwanted
  49. c eigenvalues are in the first NP locations and the wanted
  50. c portion is in the last KEV locations. When exact shifts are
  51. c selected, the unwanted part corresponds to the shifts to
  52. c be applied. Also, if ISHIFT .eq. 1, the unwanted eigenvalues
  53. c are further sorted so that the ones with largest Ritz values
  54. c are first.
  55. c
  56. c BOUNDS Double precision array of length KEV+NP. (INPUT/OUTPUT)
  57. c Error bounds corresponding to the ordering in RITZ.
  58. c
  59. c SHIFTR, SHIFTI *** USE deprecated as of version 2.1. ***
  60. c
  61. c
  62. c\EndDoc
  63. c
  64. c-----------------------------------------------------------------------
  65. c
  66. c\BeginLib
  67. c
  68. c\Local variables:
  69. c xxxxxx real
  70. c
  71. c\Routines called:
  72. c dsortc ARPACK sorting routine.
  73. c dcopy Level 1 BLAS that copies one vector to another .
  74. c
  75. c\Author
  76. c Danny Sorensen Phuong Vu
  77. c Richard Lehoucq CRPC / Rice University
  78. c Dept. of Computational & Houston, Texas
  79. c Applied Mathematics
  80. c Rice University
  81. c Houston, Texas
  82. c
  83. c\Revision history:
  84. c xx/xx/92: Version ' 2.1'
  85. c
  86. c\SCCS Information: @(#)
  87. c FILE: ngets.F SID: 2.3 DATE OF SID: 4/20/96 RELEASE: 2
  88. c
  89. c\Remarks
  90. c 1. xxxx
  91. c
  92. c\EndLib
  93. c
  94. c-----------------------------------------------------------------------
  95. c
  96. subroutine dngets ( ishift, which, kev, np, ritzr, ritzi, bounds,
  97. & shiftr, shifti )
  98. c
  99. c %----------------------------------------------------%
  100. c | Include files for debugging and timing information |
  101. -INC TARTRAK
  102. c %----------------------------------------------------%
  103. c
  104. c
  105. c %------------------%
  106. c | Scalar Arguments |
  107. c %------------------%
  108. c
  109. character*2 which
  110. integer ishift, kev, np
  111. c
  112. c %-----------------%
  113. c | Array Arguments |
  114. c %-----------------%
  115. c
  116. REAL*8
  117. & bounds(kev+np), ritzr(kev+np), ritzi(kev+np),
  118. & shiftr(1), shifti(1)
  119. c
  120. c %------------%
  121. c | Parameters |
  122. c %------------%
  123. c
  124. REAL*8
  125. & one, zero
  126. parameter (one = 1.0, zero = 0.0)
  127. c
  128. c %---------------%
  129. c | Local Scalars |
  130. c %---------------%
  131. c
  132. integer msglvl
  133. c
  134. c %----------------------%
  135. c | External Subroutines |
  136. c %----------------------%
  137. c
  138. external dcopy, dsortc, arscnd
  139. c
  140. c %----------------------%
  141. **c | Intrinsics Functions |
  142. **c %----------------------%
  143. **c
  144. ** intrinsic abs
  145. **c
  146. **c %-----------------------%
  147. **c | Executable Statements |
  148. c %-----------------------%
  149. c
  150. c %-------------------------------%
  151. c | Initialize timing statistics |
  152. c | & message level for debugging |
  153. c %-------------------------------%
  154. c
  155. * call arscnd (t0)
  156. msglvl = mngets
  157. c
  158. c %----------------------------------------------------%
  159. c | LM, SM, LR, SR, LI, SI case. |
  160. c | Sort the eigenvalues of H into the desired order |
  161. c | and apply the resulting order to BOUNDS. |
  162. c | The eigenvalues are sorted so that the wanted part |
  163. c | are always in the last KEV locations. |
  164. c | We first do a pre-processing sort in order to keep |
  165. c | complex conjugate pairs together |
  166. c %----------------------------------------------------%
  167. c
  168. if (which .eq. 'LM') then
  169. call dsortc ('LR', .true., kev+np, ritzr, ritzi, bounds)
  170. else if (which .eq. 'SM') then
  171. call dsortc ('SR', .true., kev+np, ritzr, ritzi, bounds)
  172. else if (which .eq. 'LR') then
  173. call dsortc ('LM', .true., kev+np, ritzr, ritzi, bounds)
  174. else if (which .eq. 'SR') then
  175. call dsortc ('SM', .true., kev+np, ritzr, ritzi, bounds)
  176. else if (which .eq. 'LI') then
  177. call dsortc ('LM', .true., kev+np, ritzr, ritzi, bounds)
  178. else if (which .eq. 'SI') then
  179. call dsortc ('SM', .true., kev+np, ritzr, ritzi, bounds)
  180. end if
  181. c
  182. call dsortc (which, .true., kev+np, ritzr, ritzi, bounds)
  183. c
  184. c %-------------------------------------------------------%
  185. c | Increase KEV by one if the ( ritzr(np),ritzi(np) ) |
  186. c | = ( ritzr(np+1),-ritzi(np+1) ) and ritz(np) .ne. zero |
  187. c | Accordingly decrease NP by one. In other words keep |
  188. c | complex conjugate pairs together. |
  189. c %-------------------------------------------------------%
  190. c
  191. if ( ( ritzr(np+1) - ritzr(np) ) .eq. zero
  192. & .and. ( ritzi(np+1) + ritzi(np) ) .eq. zero ) then
  193. np = np - 1
  194. kev = kev + 1
  195. end if
  196. c
  197. if ( ishift .eq. 1 ) then
  198. c
  199. c %-------------------------------------------------------%
  200. c | Sort the unwanted Ritz values used as shifts so that |
  201. c | the ones with largest Ritz estimates are first |
  202. c | This will tend to minimize the effects of the |
  203. c | forward instability of the iteration when they shifts |
  204. c | are applied in subroutine dnapps. |
  205. c | Be careful and use 'SR' since we want to sort BOUNDS! |
  206. c %-------------------------------------------------------%
  207. c
  208. call dsortc ( 'SR', .true., np, bounds, ritzr, ritzi )
  209. end if
  210. c
  211. * call arscnd (t1)
  212. tngets = tngets + (t1 - t0)
  213. c
  214. if (msglvl .gt. 0) then
  215. call ivout (logfil, 1, kev, ndigit, '_ngets: KEV is')
  216. call ivout (logfil, 1, np, ndigit, '_ngets: NP is')
  217. call dvout (logfil, kev+np, ritzr, ndigit,
  218. & '_ngets: Eigenvalues of current H matrix -- real part')
  219. call dvout (logfil, kev+np, ritzi, ndigit,
  220. & '_ngets: Eigenvalues of current H matrix -- imag part')
  221. call dvout (logfil, kev+np, bounds, ndigit,
  222. & '_ngets: Ritz estimates of the current KEV+NP Ritz values')
  223. end if
  224. c
  225. return
  226. c
  227. c %---------------%
  228. c | End of dngets |
  229. c %---------------%
  230. c
  231. end
  232.  
  233.  

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