Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/parpack/pdsgets.f
5191 views
1
c-----------------------------------------------------------------------
2
c\BeginDoc
3
c
4
c\Name: pdsgets
5
c
6
c Message Passing Layer: MPI
7
c
8
c\Description:
9
c Given the eigenvalues of the symmetric tridiagonal 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: This is called even in the case of user specified shifts in
15
c order to sort the eigenvalues, and error bounds of H for later use.
16
c
17
c\Usage:
18
c call pdsgets
19
c ( COMM, ISHIFT, WHICH, KEV, NP, RITZ, BOUNDS, SHIFTS )
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' -> KEV eigenvalues of largest magnitude are retained.
32
c 'SM' -> KEV eigenvalues of smallest magnitude are retained.
33
c 'LA' -> KEV eigenvalues of largest value are retained.
34
c 'SA' -> KEV eigenvalues of smallest value are retained.
35
c 'BE' -> KEV eigenvalues, half from each end of the spectrum.
36
c If KEV is odd, compute one more from the high end.
37
c
38
c KEV Integer. (INPUT)
39
c KEV+NP is the size of the matrix H.
40
c
41
c NP Integer. (INPUT)
42
c Number of implicit shifts to be computed.
43
c
44
c RITZ Double precision array of length KEV+NP. (INPUT/OUTPUT)
45
c On INPUT, RITZ contains the eigenvalues of H.
46
c On OUTPUT, RITZ are sorted so that the unwanted eigenvalues
47
c are in the first NP locations and the wanted part is in
48
c the last KEV locations. When exact shifts are selected, the
49
c unwanted part corresponds to the shifts to be applied.
50
c
51
c BOUNDS Double precision array of length KEV+NP. (INPUT/OUTPUT)
52
c Error bounds corresponding to the ordering in RITZ.
53
c
54
c SHIFTS Double precision array of length NP. (INPUT/OUTPUT)
55
c On INPUT: contains the user specified shifts if ISHIFT = 0.
56
c On OUTPUT: contains the shifts sorted into decreasing order
57
c of magnitude with respect to the Ritz estimates contained in
58
c BOUNDS. If ISHIFT = 0, SHIFTS is not modified on exit.
59
c
60
c\EndDoc
61
c
62
c-----------------------------------------------------------------------
63
c
64
c\BeginLib
65
c
66
c\Local variables:
67
c xxxxxx real
68
c
69
c\Routines called:
70
c dsortr ARPACK utility sorting routine.
71
c pivout Parallel ARPACK utility routine that prints integers.
72
c second ARPACK utility routine for timing.
73
c pdvout Parallel ARPACK utility routine that prints vectors.
74
c dcopy Level 1 BLAS that copies one vector to another.
75
c dswap Level 1 BLAS that swaps the contents of two vectors.
76
c
77
c\Author
78
c Danny Sorensen Phuong Vu
79
c Richard Lehoucq CRPC / Rice University
80
c Dept. of Computational & Houston, Texas
81
c Applied Mathematics
82
c Rice University
83
c Houston, Texas
84
c
85
c\Parallel Modifications
86
c Kristi Maschhoff
87
c
88
c\Revision history:
89
c Starting Point: Serial Code FILE: sgets.F SID: 2.3
90
c
91
c\SCCS Information:
92
c FILE: sgets.F SID: 1.2 DATE OF SID: 2/22/96
93
c
94
c\Remarks
95
c
96
c\EndLib
97
c
98
c-----------------------------------------------------------------------
99
c
100
subroutine pdsgets
101
& ( comm, ishift, which, kev, np, ritz, bounds, shifts )
102
c
103
c %--------------------%
104
c | MPI Communicator |
105
c %--------------------%
106
c
107
integer comm
108
c
109
c %----------------------------------------------------%
110
c | Include files for debugging and timing information |
111
c %----------------------------------------------------%
112
c
113
include 'debug.h'
114
include 'stat.h'
115
c
116
c %------------------%
117
c | Scalar Arguments |
118
c %------------------%
119
c
120
character*2 which
121
integer ishift, kev, np
122
c
123
c %-----------------%
124
c | Array Arguments |
125
c %-----------------%
126
c
127
Double precision
128
& bounds(kev+np), ritz(kev+np), shifts(np)
129
c
130
c %------------%
131
c | Parameters |
132
c %------------%
133
c
134
Double precision
135
& one, zero
136
parameter (one = 1.0, zero = 0.0)
137
c
138
c %---------------%
139
c | Local Scalars |
140
c %---------------%
141
c
142
integer kevd2, msglvl
143
c
144
c %----------------------%
145
c | External Subroutines |
146
c %----------------------%
147
c
148
external dswap, dcopy, dsortr, second
149
c
150
c %---------------------%
151
c | Intrinsic Functions |
152
c %---------------------%
153
c
154
intrinsic max, min
155
c
156
c %-----------------------%
157
c | Executable Statements |
158
c %-----------------------%
159
c
160
c %-------------------------------%
161
c | Initialize timing statistics |
162
c | & message level for debugging |
163
c %-------------------------------%
164
c
165
call second (t0)
166
msglvl = msgets
167
c
168
if (which .eq. 'BE') then
169
c
170
c %-----------------------------------------------------%
171
c | Both ends of the spectrum are requested. |
172
c | Sort the eigenvalues into algebraically increasing |
173
c | order first then swap high end of the spectrum next |
174
c | to low end in appropriate locations. |
175
c | NOTE: when np < floor(kev/2) be careful not to swap |
176
c | overlapping locations. |
177
c %-----------------------------------------------------%
178
c
179
call dsortr ('LA', .true., kev+np, ritz, bounds)
180
kevd2 = kev / 2
181
if ( kev .gt. 1 ) then
182
call dswap ( min(kevd2,np), ritz, 1,
183
& ritz( max(kevd2,np)+1 ), 1)
184
call dswap ( min(kevd2,np), bounds, 1,
185
& bounds( max(kevd2,np)+1 ), 1)
186
end if
187
c
188
else
189
c
190
c %----------------------------------------------------%
191
c | LM, SM, LA, SA case. |
192
c | Sort the eigenvalues of H into the desired order |
193
c | and apply the resulting order to BOUNDS. |
194
c | The eigenvalues are sorted so that the wanted part |
195
c | are always in the last KEV locations. |
196
c %----------------------------------------------------%
197
c
198
call dsortr (which, .true., kev+np, ritz, bounds)
199
end if
200
c
201
if (ishift .eq. 1 .and. np .gt. 0) then
202
c
203
c %-------------------------------------------------------%
204
c | Sort the unwanted Ritz values used as shifts so that |
205
c | the ones with largest Ritz estimates are first. |
206
c | This will tend to minimize the effects of the |
207
c | forward instability of the iteration when the shifts |
208
c | are applied in subroutine pdsapps. |
209
c %-------------------------------------------------------%
210
c
211
call dsortr ('SM', .true., np, bounds, ritz)
212
call dcopy (np, ritz, 1, shifts, 1)
213
end if
214
c
215
call second (t1)
216
tsgets = tsgets + (t1 - t0)
217
c
218
if (msglvl .gt. 0) then
219
call pivout (comm, logfil, 1, [kev], ndigit, '_sgets: KEV is')
220
call pivout (comm, logfil, 1, [np], ndigit, '_sgets: NP is')
221
call pdvout (comm, logfil, kev+np, ritz, ndigit,
222
& '_sgets: Eigenvalues of current H matrix')
223
call pdvout (comm, logfil, kev+np, bounds, ndigit,
224
& '_sgets: Associated Ritz estimates')
225
end if
226
c
227
return
228
c
229
c %----------------%
230
c | End of pdsgets |
231
c %----------------%
232
c
233
end
234
235