Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/parpack/pdngets.f
5194 views
1
c-----------------------------------------------------------------------
2
c\BeginDoc
3
c
4
c\Name: pdngets
5
c
6
c Message Passing Layer: MPI
7
c
8
c\Description:
9
c Given the eigenvalues of the upper Hessenberg matrix H,
10
c computes the NP shifts AMU that are zeros of the polynomial of
11
c degree NP which filters out components of the unwanted eigenvectors
12
c corresponding to the AMU's based on some given criteria.
13
c
14
c NOTE: call this even in the case of user specified shifts in order
15
c to sort the eigenvalues, and error bounds of H for later use.
16
c
17
c\Usage:
18
c call pdngets
19
c ( COMM, ISHIFT, WHICH, KEV, NP, RITZR, RITZI, BOUNDS, SHIFTR, SHIFTI )
20
c
21
c\Arguments
22
c COMM MPI Communicator for the processor grid. (INPUT)
23
c
24
c ISHIFT Integer. (INPUT)
25
c Method for selecting the implicit shifts at each iteration.
26
c ISHIFT = 0: user specified shifts
27
c ISHIFT = 1: exact shift with respect to the matrix H.
28
c
29
c WHICH Character*2. (INPUT)
30
c Shift selection criteria.
31
c 'LM' -> want the KEV eigenvalues of largest magnitude.
32
c 'SM' -> want the KEV eigenvalues of smallest magnitude.
33
c 'LR' -> want the KEV eigenvalues of largest real part.
34
c 'SR' -> want the KEV eigenvalues of smallest real part.
35
c 'LI' -> want the KEV eigenvalues of largest imaginary part.
36
c 'SI' -> want the KEV eigenvalues of smallest imaginary part.
37
c
38
c KEV Integer. (INPUT/OUTPUT)
39
c INPUT: KEV+NP is the size of the matrix H.
40
c OUTPUT: Possibly increases KEV by one to keep complex conjugate
41
c pairs together.
42
c
43
c NP Integer. (INPUT/OUTPUT)
44
c Number of implicit shifts to be computed.
45
c OUTPUT: Possibly decreases NP by one to keep complex conjugate
46
c pairs together.
47
c
48
c RITZR, Double precision array of length KEV+NP. (INPUT/OUTPUT)
49
c RITZI On INPUT, RITZR and RITZI contain the real and imaginary
50
c parts of the eigenvalues of H.
51
c On OUTPUT, RITZR and RITZI are sorted so that the unwanted
52
c eigenvalues are in the first NP locations and the wanted
53
c portion is in the last KEV locations. When exact shifts are
54
c selected, the unwanted part corresponds to the shifts to
55
c be applied. Also, if ISHIFT .eq. 1, the unwanted eigenvalues
56
c are further sorted so that the ones with largest Ritz values
57
c are first.
58
c
59
c BOUNDS Double precision array of length KEV+NP. (INPUT/OUTPUT)
60
c Error bounds corresponding to the ordering in RITZ.
61
c
62
c SHIFTR, SHIFTI *** USE deprecated as of version 2.1. ***
63
c
64
c
65
c\EndDoc
66
c
67
c-----------------------------------------------------------------------
68
c
69
c\BeginLib
70
c
71
c\Local variables:
72
c xxxxxx real
73
c
74
c\Routines called:
75
c dsortc ARPACK sorting routine.
76
c dcopy Level 1 BLAS that copies one vector to another .
77
c
78
c\Author
79
c Danny Sorensen Phuong Vu
80
c Richard Lehoucq CRPC / Rice University
81
c Dept. of Computational & Houston, Texas
82
c Applied Mathematics
83
c Rice University
84
c Houston, Texas
85
c
86
c\Parallel Modifications
87
c Kristi Maschhoff
88
c
89
c\Revision history:
90
c Starting Point: Serial Code FILE: ngets.F SID: 2.2
91
c
92
c\SCCS Information:
93
c FILE: ngets.F SID: 1.2 DATE OF SID: 2/22/96
94
c
95
c\Remarks
96
c 1. xxxx
97
c
98
c\EndLib
99
c
100
c-----------------------------------------------------------------------
101
c
102
subroutine pdngets
103
& ( comm, ishift, which, kev, np, ritzr, ritzi,
104
& bounds, shiftr, shifti )
105
c
106
c %--------------------%
107
c | MPI Communicator |
108
c %--------------------%
109
c
110
integer comm
111
c
112
c %----------------------------------------------------%
113
c | Include files for debugging and timing information |
114
c %----------------------------------------------------%
115
c
116
include 'debug.h'
117
include 'stat.h'
118
c
119
c %------------------%
120
c | Scalar Arguments |
121
c %------------------%
122
c
123
character*2 which
124
integer ishift, kev, np
125
c
126
c %-----------------%
127
c | Array Arguments |
128
c %-----------------%
129
c
130
Double precision
131
& bounds(kev+np), ritzr(kev+np), ritzi(kev+np),
132
& shiftr(1), shifti(1)
133
c
134
c %------------%
135
c | Parameters |
136
c %------------%
137
c
138
Double precision
139
& one, zero
140
parameter (one = 1.0, zero = 0.0)
141
c
142
c %---------------%
143
c | Local Scalars |
144
c %---------------%
145
c
146
integer msglvl
147
c
148
c %----------------------%
149
c | External Subroutines |
150
c %----------------------%
151
c
152
external dcopy, dsortc, second
153
c
154
c %----------------------%
155
c | Intrinsics Functions |
156
c %----------------------%
157
c
158
intrinsic abs
159
c
160
c %-----------------------%
161
c | Executable Statements |
162
c %-----------------------%
163
c
164
c %-------------------------------%
165
c | Initialize timing statistics |
166
c | & message level for debugging |
167
c %-------------------------------%
168
c
169
call second (t0)
170
msglvl = mngets
171
c
172
c %----------------------------------------------------%
173
c | LM, SM, LR, SR, LI, SI case. |
174
c | Sort the eigenvalues of H into the desired order |
175
c | and apply the resulting order to BOUNDS. |
176
c | The eigenvalues are sorted so that the wanted part |
177
c | are always in the last KEV locations. |
178
c | We first do a pre-processing sort in order to keep |
179
c | complex conjugate pairs together |
180
c %----------------------------------------------------%
181
c
182
if (which .eq. 'LM') then
183
call dsortc ('LR', .true., kev+np, ritzr, ritzi, bounds)
184
else if (which .eq. 'SM') then
185
call dsortc ('SR', .true., kev+np, ritzr, ritzi, bounds)
186
else if (which .eq. 'LR') then
187
call dsortc ('LM', .true., kev+np, ritzr, ritzi, bounds)
188
else if (which .eq. 'SR') then
189
call dsortc ('SM', .true., kev+np, ritzr, ritzi, bounds)
190
else if (which .eq. 'LI') then
191
call dsortc ('LM', .true., kev+np, ritzr, ritzi, bounds)
192
else if (which .eq. 'SI') then
193
call dsortc ('SM', .true., kev+np, ritzr, ritzi, bounds)
194
end if
195
c
196
call dsortc (which, .true., kev+np, ritzr, ritzi, bounds)
197
c
198
c %-------------------------------------------------------%
199
c | Increase KEV by one if the ( ritzr(np),ritzi(np) ) |
200
c | = ( ritzr(np+1),-ritzi(np+1) ) and ritz(np) .ne. zero |
201
c | Accordingly decrease NP by one. In other words keep |
202
c | complex conjugate pairs together. |
203
c %-------------------------------------------------------%
204
c
205
if ( ( ritzr(np+1) - ritzr(np) ) .eq. zero
206
& .and. ( ritzi(np+1) + ritzi(np) ) .eq. zero ) then
207
np = np - 1
208
kev = kev + 1
209
end if
210
c
211
if ( ishift .eq. 1 ) then
212
c
213
c %-------------------------------------------------------%
214
c | Sort the unwanted Ritz values used as shifts so that |
215
c | the ones with largest Ritz estimates are first |
216
c | This will tend to minimize the effects of the |
217
c | forward instability of the iteration when they shifts |
218
c | are applied in subroutine pdnapps. |
219
c | Be careful and use 'SR' since we want to sort BOUNDS! |
220
c %-------------------------------------------------------%
221
c
222
call dsortc ( 'SR', .true., np, bounds, ritzr, ritzi )
223
end if
224
c
225
call second (t1)
226
tngets = tngets + (t1 - t0)
227
c
228
if (msglvl .gt. 0) then
229
call pivout (comm, logfil, 1, [kev], ndigit, '_ngets: KEV is')
230
call pivout (comm, logfil, 1, [np], ndigit, '_ngets: NP is')
231
call pdvout (comm, logfil, kev+np, ritzr, ndigit,
232
& '_ngets: Eigenvalues of current H matrix -- real part')
233
call pdvout (comm, logfil, kev+np, ritzi, ndigit,
234
& '_ngets: Eigenvalues of current H matrix -- imag part')
235
call pdvout (comm, logfil, kev+np, bounds, ndigit,
236
& '_ngets: Ritz estimates of the current KEV+NP Ritz values')
237
end if
238
c
239
return
240
c
241
c %----------------%
242
c | End of pdngets |
243
c %----------------%
244
c
245
end
246
247