Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : #include <slepc/private/slepcimpl.h> /*I "slepcsys.h" I*/
12 : #include <slepcrg.h>
13 : #include <slepcst.h>
14 :
15 : /*@
16 : SlepcSCCompare - Compares two (possibly complex) values according
17 : to a certain criterion.
18 :
19 : Not Collective
20 :
21 : Input Parameters:
22 : + sc - the sorting criterion context
23 : . ar - real part of the 1st value
24 : . ai - imaginary part of the 1st value
25 : . br - real part of the 2nd value
26 : - bi - imaginary part of the 2nd value
27 :
28 : Output Parameter:
29 : . res - result of comparison
30 :
31 : Notes:
32 : Returns an integer less than, equal to, or greater than zero if the first
33 : value is considered to be respectively less than, equal to, or greater
34 : than the second one.
35 :
36 : Level: developer
37 :
38 : .seealso: SlepcSortEigenvalues(), SlepcSC
39 : @*/
40 2368872 : PetscErrorCode SlepcSCCompare(SlepcSC sc,PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *res)
41 : {
42 2368872 : PetscScalar re[2],im[2];
43 2368872 : PetscInt cin[2];
44 2368872 : PetscBool inside[2];
45 :
46 2368872 : PetscFunctionBegin;
47 2368872 : PetscAssertPointer(res,6);
48 : #if defined(PETSC_USE_DEBUG)
49 2368872 : PetscCheck(sc->comparison,PETSC_COMM_SELF,PETSC_ERR_USER,"Undefined comparison function");
50 : #endif
51 2368872 : re[0] = ar; re[1] = br;
52 2368872 : im[0] = ai; im[1] = bi;
53 2368872 : if (sc->map) PetscCall((*sc->map)(sc->mapobj,2,re,im));
54 2368872 : if (sc->rg) {
55 65407 : PetscCall(RGCheckInside(sc->rg,2,re,im,cin));
56 65407 : inside[0] = PetscNot(cin[0]<0);
57 65407 : inside[1] = PetscNot(cin[1]<0);
58 65407 : if (inside[0] && !inside[1]) *res = -1;
59 59386 : else if (!inside[0] && inside[1]) *res = 1;
60 58488 : else PetscCall((*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx));
61 2303465 : } else PetscCall((*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx));
62 2368872 : PetscFunctionReturn(PETSC_SUCCESS);
63 : }
64 :
65 : /*@
66 : SlepcSortEigenvalues - Sorts a list of eigenvalues according to the
67 : sorting criterion specified in a SlepcSC context.
68 :
69 : Not Collective
70 :
71 : Input Parameters:
72 : + sc - the sorting criterion context
73 : . n - number of eigenvalues in the list
74 : . eigr - pointer to the array containing the eigenvalues
75 : - eigi - imaginary part of the eigenvalues (only when using real numbers)
76 :
77 : Output Parameter:
78 : . perm - permutation array. Must be initialized to 0:n-1 on input.
79 :
80 : Note:
81 : The result is a list of indices in the original eigenvalue array
82 : corresponding to the first n eigenvalues sorted in the specified
83 : criterion.
84 :
85 : Level: developer
86 :
87 : .seealso: SlepcSCCompare(), SlepcSC
88 : @*/
89 1177 : PetscErrorCode SlepcSortEigenvalues(SlepcSC sc,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,PetscInt *perm)
90 : {
91 1177 : PetscScalar re,im;
92 1177 : PetscInt i,j,result,tmp;
93 :
94 1177 : PetscFunctionBegin;
95 1177 : PetscAssertPointer(sc,1);
96 1177 : PetscAssertPointer(eigr,3);
97 1177 : PetscAssertPointer(eigi,4);
98 1177 : PetscAssertPointer(perm,5);
99 : /* insertion sort */
100 11773 : for (i=n-1;i>=0;i--) {
101 10596 : re = eigr[perm[i]];
102 10596 : im = eigi[perm[i]];
103 10596 : j = i+1;
104 : #if !defined(PETSC_USE_COMPLEX)
105 10596 : if (im!=0) {
106 : /* complex eigenvalue */
107 237 : i--;
108 237 : im = eigi[perm[i]];
109 : }
110 : #endif
111 11141 : while (j<n) {
112 9813 : PetscCall(SlepcSCCompare(sc,re,im,eigr[perm[j]],eigi[perm[j]],&result));
113 9813 : if (result<=0) break;
114 : #if !defined(PETSC_USE_COMPLEX)
115 : /* keep together every complex conjugated eigenpair */
116 545 : if (!im) {
117 533 : if (eigi[perm[j]] == 0.0) {
118 : #endif
119 533 : tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = tmp;
120 533 : j++;
121 : #if !defined(PETSC_USE_COMPLEX)
122 : } else {
123 0 : tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp;
124 0 : j+=2;
125 : }
126 : } else {
127 12 : if (eigi[perm[j]] == 0.0) {
128 6 : tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = perm[j-1]; perm[j-1] = tmp;
129 6 : j++;
130 : } else {
131 6 : tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = tmp;
132 6 : tmp = perm[j-1]; perm[j-1] = perm[j+1]; perm[j+1] = tmp;
133 6 : j+=2;
134 : }
135 : }
136 : #endif
137 : }
138 : }
139 1177 : PetscFunctionReturn(PETSC_SUCCESS);
140 : }
141 :
142 : /*
143 : SlepcMap_ST - Gateway function to call STBackTransform from outside ST.
144 : */
145 2053815 : PetscErrorCode SlepcMap_ST(PetscObject obj,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
146 : {
147 2053815 : PetscFunctionBegin;
148 2053815 : PetscCall(STBackTransform((ST)obj,n,eigr,eigi));
149 2053815 : PetscFunctionReturn(PETSC_SUCCESS);
150 : }
151 :
152 854859 : PetscErrorCode SlepcCompareLargestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
153 : {
154 854859 : PetscReal a,b;
155 :
156 854859 : PetscFunctionBegin;
157 854859 : a = SlepcAbsEigenvalue(ar,ai);
158 854859 : b = SlepcAbsEigenvalue(br,bi);
159 854859 : if (a<b) *result = 1;
160 728248 : else if (a>b) *result = -1;
161 516 : else *result = 0;
162 854859 : PetscFunctionReturn(PETSC_SUCCESS);
163 : }
164 :
165 46595 : PetscErrorCode SlepcCompareSmallestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
166 : {
167 46595 : PetscReal a,b;
168 :
169 46595 : PetscFunctionBegin;
170 46595 : a = SlepcAbsEigenvalue(ar,ai);
171 46595 : b = SlepcAbsEigenvalue(br,bi);
172 46595 : if (a>b) *result = 1;
173 24976 : else if (a<b) *result = -1;
174 30 : else *result = 0;
175 46595 : PetscFunctionReturn(PETSC_SUCCESS);
176 : }
177 :
178 292090 : PetscErrorCode SlepcCompareLargestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
179 : {
180 292090 : PetscReal a,b;
181 :
182 292090 : PetscFunctionBegin;
183 292090 : a = PetscRealPart(ar);
184 292090 : b = PetscRealPart(br);
185 292090 : if (a<b) *result = 1;
186 236288 : else if (a>b) *result = -1;
187 787 : else *result = 0;
188 292090 : PetscFunctionReturn(PETSC_SUCCESS);
189 : }
190 :
191 163591 : PetscErrorCode SlepcCompareSmallestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
192 : {
193 163591 : PetscReal a,b;
194 :
195 163591 : PetscFunctionBegin;
196 163591 : a = PetscRealPart(ar);
197 163591 : b = PetscRealPart(br);
198 163591 : if (a>b) *result = 1;
199 83483 : else if (a<b) *result = -1;
200 69 : else *result = 0;
201 163591 : PetscFunctionReturn(PETSC_SUCCESS);
202 : }
203 :
204 2345 : PetscErrorCode SlepcCompareLargestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
205 : {
206 2345 : PetscReal a,b;
207 :
208 2345 : PetscFunctionBegin;
209 : #if defined(PETSC_USE_COMPLEX)
210 : a = PetscImaginaryPart(ar);
211 : b = PetscImaginaryPart(br);
212 : #else
213 2345 : a = PetscAbsReal(ai);
214 2345 : b = PetscAbsReal(bi);
215 : #endif
216 2345 : if (a<b) *result = 1;
217 2117 : else if (a>b) *result = -1;
218 1697 : else *result = 0;
219 2345 : PetscFunctionReturn(PETSC_SUCCESS);
220 : }
221 :
222 5766 : PetscErrorCode SlepcCompareSmallestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
223 : {
224 5766 : PetscReal a,b;
225 :
226 5766 : PetscFunctionBegin;
227 : #if defined(PETSC_USE_COMPLEX)
228 : a = PetscImaginaryPart(ar);
229 : b = PetscImaginaryPart(br);
230 : #else
231 5766 : a = PetscAbsReal(ai);
232 5766 : b = PetscAbsReal(bi);
233 : #endif
234 5766 : if (a>b) *result = 1;
235 5184 : else if (a<b) *result = -1;
236 4544 : else *result = 0;
237 5766 : PetscFunctionReturn(PETSC_SUCCESS);
238 : }
239 :
240 896579 : PetscErrorCode SlepcCompareTargetMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
241 : {
242 896579 : PetscReal a,b;
243 896579 : PetscScalar *target = (PetscScalar*)ctx;
244 :
245 896579 : PetscFunctionBegin;
246 : /* complex target only allowed if scalartype=complex */
247 896579 : a = SlepcAbsEigenvalue(ar-(*target),ai);
248 896579 : b = SlepcAbsEigenvalue(br-(*target),bi);
249 896579 : if (a>b) *result = 1;
250 748651 : else if (a<b) *result = -1;
251 4 : else *result = 0;
252 896579 : PetscFunctionReturn(PETSC_SUCCESS);
253 : }
254 :
255 1532 : PetscErrorCode SlepcCompareTargetReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
256 : {
257 1532 : PetscReal a,b;
258 1532 : PetscScalar *target = (PetscScalar*)ctx;
259 :
260 1532 : PetscFunctionBegin;
261 1532 : a = PetscAbsReal(PetscRealPart(ar-(*target)));
262 1532 : b = PetscAbsReal(PetscRealPart(br-(*target)));
263 1532 : if (a>b) *result = 1;
264 496 : else if (a<b) *result = -1;
265 0 : else *result = 0;
266 1532 : PetscFunctionReturn(PETSC_SUCCESS);
267 : }
268 :
269 : #if defined(PETSC_USE_COMPLEX)
270 : PetscErrorCode SlepcCompareTargetImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
271 : {
272 : PetscReal a,b;
273 : PetscScalar *target = (PetscScalar*)ctx;
274 :
275 : PetscFunctionBegin;
276 : a = PetscAbsReal(PetscImaginaryPart(ar-(*target)));
277 : b = PetscAbsReal(PetscImaginaryPart(br-(*target)));
278 : if (a>b) *result = 1;
279 : else if (a<b) *result = -1;
280 : else *result = 0;
281 : PetscFunctionReturn(PETSC_SUCCESS);
282 : }
283 : #endif
284 :
285 : /*
286 : Used in the SVD for computing smallest singular values
287 : from the cyclic matrix.
288 : */
289 8531 : PetscErrorCode SlepcCompareSmallestPosReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
290 : {
291 8531 : PetscReal a,b;
292 8531 : PetscBool aisright,bisright;
293 :
294 8531 : PetscFunctionBegin;
295 8531 : if (PetscRealPart(ar)>0.0) aisright = PETSC_TRUE;
296 1415 : else aisright = PETSC_FALSE;
297 8531 : if (PetscRealPart(br)>0.0) bisright = PETSC_TRUE;
298 6163 : else bisright = PETSC_FALSE;
299 8531 : if (aisright == bisright) { /* same sign */
300 3513 : a = SlepcAbsEigenvalue(ar,ai);
301 3513 : b = SlepcAbsEigenvalue(br,bi);
302 3513 : if (a>b) *result = 1;
303 945 : else if (a<b) *result = -1;
304 0 : else *result = 0;
305 5018 : } else if (aisright && !bisright) *result = -1; /* 'a' is on the right */
306 135 : else *result = 1; /* 'b' is on the right */
307 8531 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
|
