Numérotation des lignes :

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

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