Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/arpack/dsgets.f
5177 views
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 second 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, second
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 second (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 second (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