Télécharger dsgets.eso

Retour à la liste

Numérotation des lignes :

  1. C DSGETS SOURCE GF238795 18/02/01 21:15:21 9724
  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. real*8 T0,T1
  109. c
  110. c %-----------------%
  111. c | Array Arguments |
  112. c %-----------------%
  113. c
  114. REAL*8
  115. & bounds(kev+np), ritz(kev+np), shifts(np)
  116. c
  117. c %------------%
  118. c | Parameters |
  119. c %------------%
  120. c
  121. REAL*8
  122. & one, zero
  123. parameter (one = 1.0D+0, zero = 0.0D+0)
  124. c
  125. c %---------------%
  126. c | Local Scalars |
  127. c %---------------%
  128. c
  129. integer kevd2, msglvl
  130. c
  131. c %----------------------%
  132. c | External Subroutines |
  133. c %----------------------%
  134. c
  135. external dswap, dcopy, dsortr, arscnd
  136. c
  137. c %---------------------%
  138. **c | Intrinsic Functions |
  139. **c %---------------------%
  140. **c
  141. ** intrinsic max, min
  142. **c
  143. **c %-----------------------%
  144. **c | Executable Statements |
  145. c %-----------------------%
  146. T0=0.D0
  147. T1=0.D0
  148. c
  149. c %-------------------------------%
  150. c | Initialize timing statistics |
  151. c | & message level for debugging |
  152. c %-------------------------------%
  153. c
  154. * call arscnd (t0)
  155. msglvl = msgets
  156. c
  157. if (which .eq. 'BE') then
  158. c
  159. c %-----------------------------------------------------%
  160. c | Both ends of the spectrum are requested. |
  161. c | Sort the eigenvalues into algebraically increasing |
  162. c | order first then swap high end of the spectrum next |
  163. c | to low end in appropriate locations. |
  164. c | NOTE: when np < floor(kev/2) be careful not to swap |
  165. c | overlapping locations. |
  166. c %-----------------------------------------------------%
  167. c
  168. call dsortr ('LA', .true., kev+np, ritz, bounds)
  169. kevd2 = kev / 2
  170. if ( kev .gt. 1 ) then
  171. call dswap ( min(kevd2,np), ritz, 1,
  172. & ritz( max(kevd2,np)+1 ), 1)
  173. call dswap ( min(kevd2,np), bounds, 1,
  174. & bounds( max(kevd2,np)+1 ), 1)
  175. end if
  176. c
  177. else
  178. c
  179. c %----------------------------------------------------%
  180. c | LM, SM, LA, SA case. |
  181. c | Sort the eigenvalues of H into the desired order |
  182. c | and apply the resulting order to BOUNDS. |
  183. c | The eigenvalues are sorted so that the wanted part |
  184. c | are always in the last KEV locations. |
  185. c %----------------------------------------------------%
  186. c
  187. call dsortr (which, .true., kev+np, ritz, bounds)
  188. end if
  189. c
  190. if (ishift .eq. 1 .and. np .gt. 0) then
  191. c
  192. c %-------------------------------------------------------%
  193. c | Sort the unwanted Ritz values used as shifts so that |
  194. c | the ones with largest Ritz estimates are first. |
  195. c | This will tend to minimize the effects of the |
  196. c | forward instability of the iteration when the shifts |
  197. c | are applied in subroutine dsapps. |
  198. c %-------------------------------------------------------%
  199. c
  200. call dsortr ('SM', .true., np, bounds, ritz)
  201. call dcopy (np, ritz, 1, shifts, 1)
  202. end if
  203. c
  204. * call arscnd (t1)
  205. tsgets = tsgets + (t1 - t0)
  206. c
  207. if (msglvl .gt. 0) then
  208. call ivout (logfil, 1, kev, ndigit, '_sgets: KEV is')
  209. call ivout (logfil, 1, np, ndigit, '_sgets: NP is')
  210. call dvout (logfil, kev+np, ritz, ndigit,
  211. & '_sgets: Eigenvalues of current H matrix')
  212. call dvout (logfil, kev+np, bounds, ndigit,
  213. & '_sgets: Associated Ritz estimates')
  214. end if
  215. c
  216. return
  217. c
  218. c %---------------%
  219. c | End of dsgets |
  220. c %---------------%
  221. c
  222. end
  223.  
  224.  
  225.  

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