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 2006117 : PetscErrorCode SlepcSCCompare(SlepcSC sc,PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *res)
41 : {
42 2006117 : PetscScalar re[2],im[2];
43 2006117 : PetscInt cin[2];
44 2006117 : PetscBool inside[2];
45 :
46 2006117 : PetscFunctionBegin;
47 2006117 : PetscAssertPointer(res,6);
48 : #if defined(PETSC_USE_DEBUG)
49 2006117 : PetscCheck(sc->comparison,PETSC_COMM_SELF,PETSC_ERR_USER,"Undefined comparison function");
50 : #endif
51 2006117 : re[0] = ar; re[1] = br;
52 2006117 : im[0] = ai; im[1] = bi;
53 2006117 : if (sc->map) PetscCall((*sc->map)(sc->mapobj,2,re,im));
54 2006117 : if (sc->rg) {
55 116308 : PetscCall(RGCheckInside(sc->rg,2,re,im,cin));
56 116308 : inside[0] = PetscNot(cin[0]<0);
57 116308 : inside[1] = PetscNot(cin[1]<0);
58 116308 : if (inside[0] && !inside[1]) *res = -1;
59 103808 : else if (!inside[0] && inside[1]) *res = 1;
60 102385 : else PetscCall((*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx));
61 1889809 : } else PetscCall((*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx));
62 2006117 : 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 1186 : PetscErrorCode SlepcSortEigenvalues(SlepcSC sc,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,PetscInt *perm)
90 : {
91 1186 : PetscScalar re,im;
92 1186 : PetscInt i,j,result,tmp;
93 :
94 1186 : PetscFunctionBegin;
95 1186 : PetscAssertPointer(sc,1);
96 1186 : PetscAssertPointer(eigr,3);
97 1186 : PetscAssertPointer(eigi,4);
98 1186 : PetscAssertPointer(perm,5);
99 : /* insertion sort */
100 12255 : for (i=n-1;i>=0;i--) {
101 11069 : re = eigr[perm[i]];
102 11069 : im = eigi[perm[i]];
103 11069 : j = i+1;
104 : #if !defined(PETSC_USE_COMPLEX)
105 : if (im!=0) {
106 : /* complex eigenvalue */
107 : i--;
108 : im = eigi[perm[i]];
109 : }
110 : #endif
111 12352 : while (j<n) {
112 10894 : PetscCall(SlepcSCCompare(sc,re,im,eigr[perm[j]],eigi[perm[j]],&result));
113 10894 : if (result<=0) break;
114 : #if !defined(PETSC_USE_COMPLEX)
115 : /* keep together every complex conjugated eigenpair */
116 : if (!im) {
117 : if (eigi[perm[j]] == 0.0) {
118 : #endif
119 1283 : tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = tmp;
120 1283 : j++;
121 : #if !defined(PETSC_USE_COMPLEX)
122 : } else {
123 : tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp;
124 : j+=2;
125 : }
126 : } else {
127 : if (eigi[perm[j]] == 0.0) {
128 : tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = perm[j-1]; perm[j-1] = tmp;
129 : j++;
130 : } else {
131 : tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = tmp;
132 : tmp = perm[j-1]; perm[j-1] = perm[j+1]; perm[j+1] = tmp;
133 : j+=2;
134 : }
135 : }
136 : #endif
137 : }
138 : }
139 1186 : PetscFunctionReturn(PETSC_SUCCESS);
140 : }
141 :
142 : /*
143 : SlepcMap_ST - Gateway function to call STBackTransform from outside ST.
144 : */
145 1689279 : PetscErrorCode SlepcMap_ST(PetscObject obj,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
146 : {
147 1689279 : PetscFunctionBegin;
148 1689279 : PetscCall(STBackTransform((ST)obj,n,eigr,eigi));
149 1689279 : PetscFunctionReturn(PETSC_SUCCESS);
150 : }
151 :
152 926517 : PetscErrorCode SlepcCompareLargestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
153 : {
154 926517 : PetscReal a,b;
155 :
156 926517 : PetscFunctionBegin;
157 926517 : a = SlepcAbsEigenvalue(ar,ai);
158 926517 : b = SlepcAbsEigenvalue(br,bi);
159 926517 : if (a<b) *result = 1;
160 787548 : else if (a>b) *result = -1;
161 2616 : else *result = 0;
162 926517 : PetscFunctionReturn(PETSC_SUCCESS);
163 : }
164 :
165 9921 : PetscErrorCode SlepcCompareSmallestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
166 : {
167 9921 : PetscReal a,b;
168 :
169 9921 : PetscFunctionBegin;
170 9921 : a = SlepcAbsEigenvalue(ar,ai);
171 9921 : b = SlepcAbsEigenvalue(br,bi);
172 9921 : if (a>b) *result = 1;
173 8696 : else if (a<b) *result = -1;
174 0 : else *result = 0;
175 9921 : PetscFunctionReturn(PETSC_SUCCESS);
176 : }
177 :
178 321608 : PetscErrorCode SlepcCompareLargestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
179 : {
180 321608 : PetscReal a,b;
181 :
182 321608 : PetscFunctionBegin;
183 321608 : a = PetscRealPart(ar);
184 321608 : b = PetscRealPart(br);
185 321608 : if (a<b) *result = 1;
186 258930 : else if (a>b) *result = -1;
187 781 : else *result = 0;
188 321608 : PetscFunctionReturn(PETSC_SUCCESS);
189 : }
190 :
191 173097 : PetscErrorCode SlepcCompareSmallestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
192 : {
193 173097 : PetscReal a,b;
194 :
195 173097 : PetscFunctionBegin;
196 173097 : a = PetscRealPart(ar);
197 173097 : b = PetscRealPart(br);
198 173097 : if (a>b) *result = 1;
199 92830 : else if (a<b) *result = -1;
200 56 : else *result = 0;
201 173097 : PetscFunctionReturn(PETSC_SUCCESS);
202 : }
203 :
204 4174 : PetscErrorCode SlepcCompareLargestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
205 : {
206 4174 : PetscReal a,b;
207 :
208 4174 : PetscFunctionBegin;
209 : #if defined(PETSC_USE_COMPLEX)
210 4174 : a = PetscImaginaryPart(ar);
211 4174 : b = PetscImaginaryPart(br);
212 : #else
213 : a = PetscAbsReal(ai);
214 : b = PetscAbsReal(bi);
215 : #endif
216 4174 : if (a<b) *result = 1;
217 3207 : else if (a>b) *result = -1;
218 0 : else *result = 0;
219 4174 : PetscFunctionReturn(PETSC_SUCCESS);
220 : }
221 :
222 4417 : PetscErrorCode SlepcCompareSmallestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
223 : {
224 4417 : PetscReal a,b;
225 :
226 4417 : PetscFunctionBegin;
227 : #if defined(PETSC_USE_COMPLEX)
228 4417 : a = PetscImaginaryPart(ar);
229 4417 : b = PetscImaginaryPart(br);
230 : #else
231 : a = PetscAbsReal(ai);
232 : b = PetscAbsReal(bi);
233 : #endif
234 4417 : if (a>b) *result = 1;
235 3489 : else if (a<b) *result = -1;
236 0 : else *result = 0;
237 4417 : PetscFunctionReturn(PETSC_SUCCESS);
238 : }
239 :
240 405066 : PetscErrorCode SlepcCompareTargetMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
241 : {
242 405066 : PetscReal a,b;
243 405066 : PetscScalar *target = (PetscScalar*)ctx;
244 :
245 405066 : PetscFunctionBegin;
246 : /* complex target only allowed if scalartype=complex */
247 405066 : a = SlepcAbsEigenvalue(ar-(*target),ai);
248 405066 : b = SlepcAbsEigenvalue(br-(*target),bi);
249 405066 : if (a>b) *result = 1;
250 269990 : else if (a<b) *result = -1;
251 15 : else *result = 0;
252 405066 : PetscFunctionReturn(PETSC_SUCCESS);
253 : }
254 :
255 3301 : PetscErrorCode SlepcCompareTargetReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
256 : {
257 3301 : PetscReal a,b;
258 3301 : PetscScalar *target = (PetscScalar*)ctx;
259 :
260 3301 : PetscFunctionBegin;
261 3301 : a = PetscAbsReal(PetscRealPart(ar-(*target)));
262 3301 : b = PetscAbsReal(PetscRealPart(br-(*target)));
263 3301 : if (a>b) *result = 1;
264 888 : else if (a<b) *result = -1;
265 2 : else *result = 0;
266 3301 : PetscFunctionReturn(PETSC_SUCCESS);
267 : }
268 :
269 : #if defined(PETSC_USE_COMPLEX)
270 2841 : PetscErrorCode SlepcCompareTargetImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
271 : {
272 2841 : PetscReal a,b;
273 2841 : PetscScalar *target = (PetscScalar*)ctx;
274 :
275 2841 : PetscFunctionBegin;
276 2841 : a = PetscAbsReal(PetscImaginaryPart(ar-(*target)));
277 2841 : b = PetscAbsReal(PetscImaginaryPart(br-(*target)));
278 2841 : if (a>b) *result = 1;
279 2118 : else if (a<b) *result = -1;
280 0 : else *result = 0;
281 2841 : 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 9542 : PetscErrorCode SlepcCompareSmallestPosReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
290 : {
291 9542 : PetscReal a,b;
292 9542 : PetscBool aisright,bisright;
293 :
294 9542 : PetscFunctionBegin;
295 9542 : if (PetscRealPart(ar)>0.0) aisright = PETSC_TRUE;
296 1543 : else aisright = PETSC_FALSE;
297 9542 : if (PetscRealPart(br)>0.0) bisright = PETSC_TRUE;
298 6874 : else bisright = PETSC_FALSE;
299 9542 : if (aisright == bisright) { /* same sign */
300 3913 : a = SlepcAbsEigenvalue(ar,ai);
301 3913 : b = SlepcAbsEigenvalue(br,bi);
302 3913 : if (a>b) *result = 1;
303 1032 : else if (a<b) *result = -1;
304 0 : else *result = 0;
305 5629 : } else if (aisright && !bisright) *result = -1; /* 'a' is on the right */
306 149 : else *result = 1; /* 'b' is on the right */
307 9542 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
|
