Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/arpack/dngets.f
5191 views
1
c-----------------------------------------------------------------------
2
c\BeginDoc
3
c
4
c\Name: dngets
5
c
6
c\Description:
7
c Given the eigenvalues of the upper Hessenberg 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: call this even in the case of user specified shifts in order
13
c to sort the eigenvalues, and error bounds of H for later use.
14
c
15
c\Usage:
16
c call dngets
17
c ( ISHIFT, WHICH, KEV, NP, RITZR, RITZI, BOUNDS, SHIFTR, SHIFTI )
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' -> want the KEV eigenvalues of largest magnitude.
28
c 'SM' -> want the KEV eigenvalues of smallest magnitude.
29
c 'LR' -> want the KEV eigenvalues of largest real part.
30
c 'SR' -> want the KEV eigenvalues of smallest real part.
31
c 'LI' -> want the KEV eigenvalues of largest imaginary part.
32
c 'SI' -> want the KEV eigenvalues of smallest imaginary part.
33
c
34
c KEV Integer. (INPUT/OUTPUT)
35
c INPUT: KEV+NP is the size of the matrix H.
36
c OUTPUT: Possibly increases KEV by one to keep complex conjugate
37
c pairs together.
38
c
39
c NP Integer. (INPUT/OUTPUT)
40
c Number of implicit shifts to be computed.
41
c OUTPUT: Possibly decreases NP by one to keep complex conjugate
42
c pairs together.
43
c
44
c RITZR, Double precision array of length KEV+NP. (INPUT/OUTPUT)
45
c RITZI On INPUT, RITZR and RITZI contain the real and imaginary
46
c parts of the eigenvalues of H.
47
c On OUTPUT, RITZR and RITZI are sorted so that the unwanted
48
c eigenvalues are in the first NP locations and the wanted
49
c portion is in the last KEV locations. When exact shifts are
50
c selected, the unwanted part corresponds to the shifts to
51
c be applied. Also, if ISHIFT .eq. 1, the unwanted eigenvalues
52
c are further sorted so that the ones with largest Ritz values
53
c are first.
54
c
55
c BOUNDS Double precision array of length KEV+NP. (INPUT/OUTPUT)
56
c Error bounds corresponding to the ordering in RITZ.
57
c
58
c SHIFTR, SHIFTI *** USE deprecated as of version 2.1. ***
59
c
60
c
61
c\EndDoc
62
c
63
c-----------------------------------------------------------------------
64
c
65
c\BeginLib
66
c
67
c\Local variables:
68
c xxxxxx real
69
c
70
c\Routines called:
71
c dsortc ARPACK sorting routine.
72
c dcopy Level 1 BLAS that copies one vector to another .
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/92: Version ' 2.1'
84
c
85
c\SCCS Information: @(#)
86
c FILE: ngets.F SID: 2.3 DATE OF SID: 4/20/96 RELEASE: 2
87
c
88
c\Remarks
89
c 1. xxxx
90
c
91
c\EndLib
92
c
93
c-----------------------------------------------------------------------
94
c
95
subroutine dngets ( ishift, which, kev, np, ritzr, ritzi, bounds,
96
& shiftr, shifti )
97
c
98
c %----------------------------------------------------%
99
c | Include files for debugging and timing information |
100
c %----------------------------------------------------%
101
c
102
include 'debug.h'
103
include 'stat.h'
104
c
105
c %------------------%
106
c | Scalar Arguments |
107
c %------------------%
108
c
109
character*2 which
110
integer ishift, kev, np
111
c
112
c %-----------------%
113
c | Array Arguments |
114
c %-----------------%
115
c
116
Double precision
117
& bounds(kev+np), ritzr(kev+np), ritzi(kev+np),
118
& shiftr(1), shifti(1)
119
c
120
c %------------%
121
c | Parameters |
122
c %------------%
123
c
124
Double precision
125
& one, zero
126
parameter (one = 1.0, zero = 0.0)
127
c
128
c %---------------%
129
c | Local Scalars |
130
c %---------------%
131
c
132
integer msglvl
133
c
134
c %----------------------%
135
c | External Subroutines |
136
c %----------------------%
137
c
138
external dcopy, dsortc, second
139
c
140
c %----------------------%
141
c | Intrinsics Functions |
142
c %----------------------%
143
c
144
intrinsic abs
145
c
146
c %-----------------------%
147
c | Executable Statements |
148
c %-----------------------%
149
c
150
c %-------------------------------%
151
c | Initialize timing statistics |
152
c | & message level for debugging |
153
c %-------------------------------%
154
c
155
call second (t0)
156
msglvl = mngets
157
c
158
c %----------------------------------------------------%
159
c | LM, SM, LR, SR, LI, SI case. |
160
c | Sort the eigenvalues of H into the desired order |
161
c | and apply the resulting order to BOUNDS. |
162
c | The eigenvalues are sorted so that the wanted part |
163
c | are always in the last KEV locations. |
164
c | We first do a pre-processing sort in order to keep |
165
c | complex conjugate pairs together |
166
c %----------------------------------------------------%
167
c
168
if (which .eq. 'LM') then
169
call dsortc ('LR', .true., kev+np, ritzr, ritzi, bounds)
170
else if (which .eq. 'SM') then
171
call dsortc ('SR', .true., kev+np, ritzr, ritzi, bounds)
172
else if (which .eq. 'LR') then
173
call dsortc ('LM', .true., kev+np, ritzr, ritzi, bounds)
174
else if (which .eq. 'SR') then
175
call dsortc ('SM', .true., kev+np, ritzr, ritzi, bounds)
176
else if (which .eq. 'LI') then
177
call dsortc ('LM', .true., kev+np, ritzr, ritzi, bounds)
178
else if (which .eq. 'SI') then
179
call dsortc ('SM', .true., kev+np, ritzr, ritzi, bounds)
180
end if
181
c
182
call dsortc (which, .true., kev+np, ritzr, ritzi, bounds)
183
c
184
c %-------------------------------------------------------%
185
c | Increase KEV by one if the ( ritzr(np),ritzi(np) ) |
186
c | = ( ritzr(np+1),-ritzi(np+1) ) and ritz(np) .ne. zero |
187
c | Accordingly decrease NP by one. In other words keep |
188
c | complex conjugate pairs together. |
189
c %-------------------------------------------------------%
190
c
191
if ( ( ritzr(np+1) - ritzr(np) ) .eq. zero
192
& .and. ( ritzi(np+1) + ritzi(np) ) .eq. zero ) then
193
np = np - 1
194
kev = kev + 1
195
end if
196
c
197
if ( ishift .eq. 1 ) then
198
c
199
c %-------------------------------------------------------%
200
c | Sort the unwanted Ritz values used as shifts so that |
201
c | the ones with largest Ritz estimates are first |
202
c | This will tend to minimize the effects of the |
203
c | forward instability of the iteration when they shifts |
204
c | are applied in subroutine dnapps. |
205
c | Be careful and use 'SR' since we want to sort BOUNDS! |
206
c %-------------------------------------------------------%
207
c
208
call dsortc ( 'SR', .true., np, bounds, ritzr, ritzi )
209
end if
210
c
211
call second (t1)
212
tngets = tngets + (t1 - t0)
213
c
214
if (msglvl .gt. 0) then
215
call ivout (logfil, 1, [kev], ndigit, '_ngets: KEV is')
216
call ivout (logfil, 1, [np], ndigit, '_ngets: NP is')
217
call dvout (logfil, kev+np, ritzr, ndigit,
218
& '_ngets: Eigenvalues of current H matrix -- real part')
219
call dvout (logfil, kev+np, ritzi, ndigit,
220
& '_ngets: Eigenvalues of current H matrix -- imag part')
221
call dvout (logfil, kev+np, bounds, ndigit,
222
& '_ngets: Ritz estimates of the current KEV+NP Ritz values')
223
end if
224
c
225
return
226
c
227
c %---------------%
228
c | End of dngets |
229
c %---------------%
230
c
231
end
232
233