Télécharger dsgets.eso

Retour à la liste

Numérotation des lignes :

  1. C DSGETS SOURCE BP208322 15/10/13 21:15:54 8670
  2. c-----------------------------------------------------------------------
  3. c\BeginDoc
  4. c
  5. c\Name: dsgets
  6. c
  7. c\Description:
  8. c Given the eigenvalues of the symmetric tridiagonal 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: This is called even in the case of user specified shifts in
  14. c order to sort the eigenvalues, and error bounds of H for later use.
  15. c
  16. c\Usage:
  17. c call dsgets
  18. c ( ISHIFT, WHICH, KEV, NP, RITZ, BOUNDS, SHIFTS )
  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' -> KEV eigenvalues of largest magnitude are retained.
  29. c 'SM' -> KEV eigenvalues of smallest magnitude are retained.
  30. c 'LA' -> KEV eigenvalues of largest value are retained.
  31. c 'SA' -> KEV eigenvalues of smallest value are retained.
  32. c 'BE' -> KEV eigenvalues, half from each end of the spectrum.
  33. c If KEV is odd, compute one more from the high end.
  34. c
  35. c KEV Integer. (INPUT)
  36. c KEV+NP is the size of the matrix H.
  37. c
  38. c NP Integer. (INPUT)
  39. c Number of implicit shifts to be computed.
  40. c
  41. c RITZ Double precision array of length KEV+NP. (INPUT/OUTPUT)
  42. c On INPUT, RITZ contains the eigenvalues of H.
  43. c On OUTPUT, RITZ are sorted so that the unwanted eigenvalues
  44. c are in the first NP locations and the wanted part is in
  45. c the last KEV locations. When exact shifts are selected, the
  46. c unwanted part corresponds to the shifts to be applied.
  47. c
  48. c BOUNDS Double precision array of length KEV+NP. (INPUT/OUTPUT)
  49. c Error bounds corresponding to the ordering in RITZ.
  50. c
  51. c SHIFTS Double precision array of length NP. (INPUT/OUTPUT)
  52. c On INPUT: contains the user specified shifts if ISHIFT = 0.
  53. c On OUTPUT: contains the shifts sorted into decreasing order
  54. c of magnitude with respect to the Ritz estimates contained in
  55. c BOUNDS. If ISHIFT = 0, SHIFTS is not modified on exit.
  56. c
  57. c\EndDoc
  58. c
  59. c-----------------------------------------------------------------------
  60. c
  61. c\BeginLib
  62. c
  63. c\Local variables:
  64. c xxxxxx real
  65. c
  66. c\Routines called:
  67. c dsortr ARPACK utility sorting routine.
  68. c ivout ARPACK utility routine that prints integers.
  69. c arscnd ARPACK utility routine for timing.
  70. c dvout ARPACK utility routine that prints vectors.
  71. c dcopy Level 1 BLAS that copies one vector to another.
  72. c dswap Level 1 BLAS that swaps the contents of two vectors.
  73. c
  74. c\Author
  75. c Danny Sorensen Phuong Vu
  76. c Richard Lehoucq CRPC / Rice University
  77. c Dept. of Computational & Houston, Texas
  78. c Applied Mathematics
  79. c Rice University
  80. c Houston, Texas
  81. c
  82. c\Revision history:
  83. c xx/xx/93: Version ' 2.1'
  84. c
  85. c\SCCS Information: @(#)
  86. c FILE: sgets.F SID: 2.4 DATE OF SID: 4/19/96 RELEASE: 2
  87. c
  88. c\Remarks
  89. c
  90. c\EndLib
  91. c
  92. c-----------------------------------------------------------------------
  93. c
  94. subroutine dsgets ( ishift, which, kev, np, ritz, bounds, shifts )
  95. c
  96. c %----------------------------------------------------%
  97. c | Include files for debugging and timing information |
  98. -INC TARTRAK
  99. c %----------------------------------------------------%
  100. c
  101. c
  102. c %------------------%
  103. c | Scalar Arguments |
  104. c %------------------%
  105. c
  106. character*2 which
  107. integer ishift, kev, np
  108. c
  109. c %-----------------%
  110. c | Array Arguments |
  111. c %-----------------%
  112. c
  113. REAL*8
  114. & bounds(kev+np), ritz(kev+np), shifts(np)
  115. c
  116. c %------------%
  117. c | Parameters |
  118. c %------------%
  119. c
  120. REAL*8
  121. & one, zero
  122. parameter (one = 1.0D+0, zero = 0.0D+0)
  123. c
  124. c %---------------%
  125. c | Local Scalars |
  126. c %---------------%
  127. c
  128. integer kevd2, msglvl
  129. c
  130. c %----------------------%
  131. c | External Subroutines |
  132. c %----------------------%
  133. c
  134. external dswap, dcopy, dsortr, arscnd
  135. c
  136. c %---------------------%
  137. **c | Intrinsic Functions |
  138. **c %---------------------%
  139. **c
  140. ** intrinsic max, min
  141. **c
  142. **c %-----------------------%
  143. **c | Executable Statements |
  144. c %-----------------------%
  145. c
  146. c %-------------------------------%
  147. c | Initialize timing statistics |
  148. c | & message level for debugging |
  149. c %-------------------------------%
  150. c
  151. * call arscnd (t0)
  152. msglvl = msgets
  153. c
  154. if (which .eq. 'BE') then
  155. c
  156. c %-----------------------------------------------------%
  157. c | Both ends of the spectrum are requested. |
  158. c | Sort the eigenvalues into algebraically increasing |
  159. c | order first then swap high end of the spectrum next |
  160. c | to low end in appropriate locations. |
  161. c | NOTE: when np < floor(kev/2) be careful not to swap |
  162. c | overlapping locations. |
  163. c %-----------------------------------------------------%
  164. c
  165. call dsortr ('LA', .true., kev+np, ritz, bounds)
  166. kevd2 = kev / 2
  167. if ( kev .gt. 1 ) then
  168. call dswap ( min(kevd2,np), ritz, 1,
  169. & ritz( max(kevd2,np)+1 ), 1)
  170. call dswap ( min(kevd2,np), bounds, 1,
  171. & bounds( max(kevd2,np)+1 ), 1)
  172. end if
  173. c
  174. else
  175. c
  176. c %----------------------------------------------------%
  177. c | LM, SM, LA, SA case. |
  178. c | Sort the eigenvalues of H into the desired order |
  179. c | and apply the resulting order to BOUNDS. |
  180. c | The eigenvalues are sorted so that the wanted part |
  181. c | are always in the last KEV locations. |
  182. c %----------------------------------------------------%
  183. c
  184. call dsortr (which, .true., kev+np, ritz, bounds)
  185. end if
  186. c
  187. if (ishift .eq. 1 .and. np .gt. 0) then
  188. c
  189. c %-------------------------------------------------------%
  190. c | Sort the unwanted Ritz values used as shifts so that |
  191. c | the ones with largest Ritz estimates are first. |
  192. c | This will tend to minimize the effects of the |
  193. c | forward instability of the iteration when the shifts |
  194. c | are applied in subroutine dsapps. |
  195. c %-------------------------------------------------------%
  196. c
  197. call dsortr ('SM', .true., np, bounds, ritz)
  198. call dcopy (np, ritz, 1, shifts, 1)
  199. end if
  200. c
  201. * call arscnd (t1)
  202. tsgets = tsgets + (t1 - t0)
  203. c
  204. if (msglvl .gt. 0) then
  205. call ivout (logfil, 1, kev, ndigit, '_sgets: KEV is')
  206. call ivout (logfil, 1, np, ndigit, '_sgets: NP is')
  207. call dvout (logfil, kev+np, ritz, ndigit,
  208. & '_sgets: Eigenvalues of current H matrix')
  209. call dvout (logfil, kev+np, bounds, ndigit,
  210. & '_sgets: Associated Ritz estimates')
  211. end if
  212. c
  213. return
  214. c
  215. c %---------------%
  216. c | End of dsgets |
  217. c %---------------%
  218. c
  219. end
  220.  
  221.  

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