Mercurial > octave
comparison libcruft/arpack/src/dsgets.f @ 12274:9f5d2ef078e8 release-3-4-x
import ARPACK sources to libcruft from Debian package libarpack2 2.1+parpack96.dfsg-3+b1
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 28 Jan 2011 14:04:33 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
12273:83133b5bf392 | 12274:9f5d2ef078e8 |
---|---|
1 c----------------------------------------------------------------------- | |
2 c\BeginDoc | |
3 c | |
4 c\Name: dsgets | |
5 c | |
6 c\Description: | |
7 c Given the eigenvalues of the symmetric tridiagonal matrix H, | |
8 c computes the NP shifts AMU that are zeros of the polynomial of | |
9 c degree NP which filters out components of the unwanted eigenvectors | |
10 c corresponding to the AMU's based on some given criteria. | |
11 c | |
12 c NOTE: This is called even in the case of user specified shifts in | |
13 c order to sort the eigenvalues, and error bounds of H for later use. | |
14 c | |
15 c\Usage: | |
16 c call dsgets | |
17 c ( ISHIFT, WHICH, KEV, NP, RITZ, BOUNDS, SHIFTS ) | |
18 c | |
19 c\Arguments | |
20 c ISHIFT Integer. (INPUT) | |
21 c Method for selecting the implicit shifts at each iteration. | |
22 c ISHIFT = 0: user specified shifts | |
23 c ISHIFT = 1: exact shift with respect to the matrix H. | |
24 c | |
25 c WHICH Character*2. (INPUT) | |
26 c Shift selection criteria. | |
27 c 'LM' -> KEV eigenvalues of largest magnitude are retained. | |
28 c 'SM' -> KEV eigenvalues of smallest magnitude are retained. | |
29 c 'LA' -> KEV eigenvalues of largest value are retained. | |
30 c 'SA' -> KEV eigenvalues of smallest value are retained. | |
31 c 'BE' -> KEV eigenvalues, half from each end of the spectrum. | |
32 c If KEV is odd, compute one more from the high end. | |
33 c | |
34 c KEV Integer. (INPUT) | |
35 c KEV+NP is the size of the matrix H. | |
36 c | |
37 c NP Integer. (INPUT) | |
38 c Number of implicit shifts to be computed. | |
39 c | |
40 c RITZ Double precision array of length KEV+NP. (INPUT/OUTPUT) | |
41 c On INPUT, RITZ contains the eigenvalues of H. | |
42 c On OUTPUT, RITZ are sorted so that the unwanted eigenvalues | |
43 c are in the first NP locations and the wanted part is in | |
44 c the last KEV locations. When exact shifts are selected, the | |
45 c unwanted part corresponds to the shifts to be applied. | |
46 c | |
47 c BOUNDS Double precision array of length KEV+NP. (INPUT/OUTPUT) | |
48 c Error bounds corresponding to the ordering in RITZ. | |
49 c | |
50 c SHIFTS Double precision array of length NP. (INPUT/OUTPUT) | |
51 c On INPUT: contains the user specified shifts if ISHIFT = 0. | |
52 c On OUTPUT: contains the shifts sorted into decreasing order | |
53 c of magnitude with respect to the Ritz estimates contained in | |
54 c BOUNDS. If ISHIFT = 0, SHIFTS is not modified on exit. | |
55 c | |
56 c\EndDoc | |
57 c | |
58 c----------------------------------------------------------------------- | |
59 c | |
60 c\BeginLib | |
61 c | |
62 c\Local variables: | |
63 c xxxxxx real | |
64 c | |
65 c\Routines called: | |
66 c dsortr ARPACK utility sorting routine. | |
67 c ivout ARPACK utility routine that prints integers. | |
68 c arscnd ARPACK utility routine for timing. | |
69 c dvout ARPACK utility routine that prints vectors. | |
70 c dcopy Level 1 BLAS that copies one vector to another. | |
71 c dswap Level 1 BLAS that swaps the contents of two vectors. | |
72 c | |
73 c\Author | |
74 c Danny Sorensen Phuong Vu | |
75 c Richard Lehoucq CRPC / Rice University | |
76 c Dept. of Computational & Houston, Texas | |
77 c Applied Mathematics | |
78 c Rice University | |
79 c Houston, Texas | |
80 c | |
81 c\Revision history: | |
82 c xx/xx/93: Version ' 2.1' | |
83 c | |
84 c\SCCS Information: @(#) | |
85 c FILE: sgets.F SID: 2.4 DATE OF SID: 4/19/96 RELEASE: 2 | |
86 c | |
87 c\Remarks | |
88 c | |
89 c\EndLib | |
90 c | |
91 c----------------------------------------------------------------------- | |
92 c | |
93 subroutine dsgets ( ishift, which, kev, np, ritz, bounds, shifts ) | |
94 c | |
95 c %----------------------------------------------------% | |
96 c | Include files for debugging and timing information | | |
97 c %----------------------------------------------------% | |
98 c | |
99 include 'debug.h' | |
100 include 'stat.h' | |
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 Double precision | |
114 & bounds(kev+np), ritz(kev+np), shifts(np) | |
115 c | |
116 c %------------% | |
117 c | Parameters | | |
118 c %------------% | |
119 c | |
120 Double precision | |
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 |