Actual source code: slepcsc.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/slepcimpl.h>
12: #include <slepcrg.h>
13: #include <slepcst.h>
15: /*@
16: SlepcSCCompare - Compares two (possibly complex) values according
17: to a certain criterion.
19: Not Collective
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
28: Output Parameter:
29: . res - result of comparison
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.
36: Level: developer
38: .seealso: SlepcSortEigenvalues(), SlepcSC
39: @*/
40: PetscErrorCode SlepcSCCompare(SlepcSC sc,PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *res)
41: {
42: PetscScalar re[2],im[2];
43: PetscInt cin[2];
44: PetscBool inside[2];
46: PetscFunctionBegin;
47: PetscAssertPointer(res,6);
48: #if defined(PETSC_USE_DEBUG)
49: PetscCheck(sc->comparison,PETSC_COMM_SELF,PETSC_ERR_USER,"Undefined comparison function");
50: #endif
51: re[0] = ar; re[1] = br;
52: im[0] = ai; im[1] = bi;
53: if (sc->map) PetscCall((*sc->map)(sc->mapobj,2,re,im));
54: if (sc->rg) {
55: PetscCall(RGCheckInside(sc->rg,2,re,im,cin));
56: inside[0] = PetscNot(cin[0]<0);
57: inside[1] = PetscNot(cin[1]<0);
58: if (inside[0] && !inside[1]) *res = -1;
59: else if (!inside[0] && inside[1]) *res = 1;
60: else PetscCall((*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx));
61: } else PetscCall((*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode SlepcSortEigenvalues_Private(SlepcSC sc,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,PetscInt *perm,PetscBool flg)
66: {
67: PetscScalar re,im;
68: PetscInt i,j,result,tmp;
70: PetscFunctionBegin;
71: /* insertion sort */
72: for (i=n-1;i>=0;i--) {
73: re = eigr[perm[i]];
74: im = eigi[perm[i]];
75: j = i+1;
76: #if !defined(PETSC_USE_COMPLEX)
77: if (im!=0 && (re!=0 || !flg)) {
78: /* complex eigenvalue */
79: i--;
80: im = eigi[perm[i]];
81: }
82: #endif
83: while (j<n) {
84: PetscCall(SlepcSCCompare(sc,re,im,eigr[perm[j]],eigi[perm[j]],&result));
85: if (result<=0) break;
86: #if !defined(PETSC_USE_COMPLEX)
87: /* keep together every complex conjugated eigenpair */
88: if (!im || (!re && flg)) {
89: if (eigi[perm[j]] == 0.0 || (flg && eigr[perm[j]] == 0.0)) {
90: #endif
91: tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = tmp;
92: j++;
93: #if !defined(PETSC_USE_COMPLEX)
94: } else {
95: tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp;
96: j+=2;
97: }
98: } else {
99: if (eigi[perm[j]] == 0.0 || (flg && eigr[perm[j]] == 0.0)) {
100: tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = perm[j-1]; perm[j-1] = tmp;
101: j++;
102: } else {
103: tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = tmp;
104: tmp = perm[j-1]; perm[j-1] = perm[j+1]; perm[j+1] = tmp;
105: j+=2;
106: }
107: }
108: #endif
109: }
110: }
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
113: /*@
114: SlepcSortEigenvalues - Sorts a list of eigenvalues according to the
115: sorting criterion specified in a SlepcSC context.
117: Not Collective
119: Input Parameters:
120: + sc - the sorting criterion context
121: . n - number of eigenvalues in the list
122: . eigr - pointer to the array containing the eigenvalues
123: - eigi - imaginary part of the eigenvalues (only when using real numbers)
125: Output Parameter:
126: . perm - permutation array. Must be initialized to 0:n-1 on input.
128: Notes:
129: The result is a list of indices in the original eigenvalue array
130: corresponding to the first n eigenvalues sorted in the specified
131: criterion.
133: In real scalars, this functions assumes that complex values come in
134: conjugate pairs that are consecutive (including purely imaginary ones).
136: Level: developer
138: .seealso: SlepcSCCompare(), SlepcSC
139: @*/
140: PetscErrorCode SlepcSortEigenvalues(SlepcSC sc,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,PetscInt *perm)
141: {
142: PetscFunctionBegin;
143: PetscAssertPointer(sc,1);
144: PetscAssertPointer(eigr,3);
145: PetscAssertPointer(eigi,4);
146: PetscAssertPointer(perm,5);
147: PetscCall(SlepcSortEigenvalues_Private(sc,n,eigr,eigi,perm,PETSC_FALSE));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@
152: SlepcSortEigenvaluesSpecial - Sorts a list of eigenvalues according to the
153: sorting criterion specified in a SlepcSC context, with a special assumption
154: on the input values.
156: Not Collective
158: Input Parameters:
159: + sc - the sorting criterion context
160: . n - number of eigenvalues in the list
161: . eigr - pointer to the array containing the eigenvalues
162: - eigi - imaginary part of the eigenvalues (only when using real numbers)
164: Output Parameter:
165: . perm - permutation array. Must be initialized to 0:n-1 on input.
167: Notes:
168: The result is a list of indices in the original eigenvalue array
169: corresponding to the first n eigenvalues sorted in the specified
170: criterion.
172: In real scalars, this functions assumes that complex values come in
173: conjugate pairs that are consecutive, but not purely imaginary ones in which
174: case only the one with positive imaginary part appears.
176: Level: developer
178: .seealso: SlepcSCCompare(), SlepcSC
179: @*/
180: PetscErrorCode SlepcSortEigenvaluesSpecial(SlepcSC sc,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,PetscInt *perm)
181: {
182: PetscFunctionBegin;
183: PetscAssertPointer(sc,1);
184: PetscAssertPointer(eigr,3);
185: PetscAssertPointer(eigi,4);
186: PetscAssertPointer(perm,5);
187: PetscCall(SlepcSortEigenvalues_Private(sc,n,eigr,eigi,perm,PETSC_TRUE));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: /*
192: SlepcMap_ST - Gateway function to call STBackTransform from outside ST.
193: */
194: PetscErrorCode SlepcMap_ST(PetscObject obj,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
195: {
196: PetscFunctionBegin;
197: PetscCall(STBackTransform((ST)obj,n,eigr,eigi));
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: PetscErrorCode SlepcCompareLargestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
202: {
203: PetscReal a,b;
205: PetscFunctionBegin;
206: a = SlepcAbsEigenvalue(ar,ai);
207: b = SlepcAbsEigenvalue(br,bi);
208: if (a<b) *result = 1;
209: else if (a>b) *result = -1;
210: else *result = 0;
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: PetscErrorCode SlepcCompareSmallestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
215: {
216: PetscReal a,b;
218: PetscFunctionBegin;
219: a = SlepcAbsEigenvalue(ar,ai);
220: b = SlepcAbsEigenvalue(br,bi);
221: if (a>b) *result = 1;
222: else if (a<b) *result = -1;
223: else *result = 0;
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: PetscErrorCode SlepcCompareLargestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
228: {
229: PetscReal a,b;
231: PetscFunctionBegin;
232: a = PetscRealPart(ar);
233: b = PetscRealPart(br);
234: if (a<b) *result = 1;
235: else if (a>b) *result = -1;
236: else *result = 0;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: PetscErrorCode SlepcCompareSmallestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
241: {
242: PetscReal a,b;
244: PetscFunctionBegin;
245: a = PetscRealPart(ar);
246: b = PetscRealPart(br);
247: if (a>b) *result = 1;
248: else if (a<b) *result = -1;
249: else *result = 0;
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: PetscErrorCode SlepcCompareLargestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
254: {
255: PetscReal a,b;
257: PetscFunctionBegin;
258: #if defined(PETSC_USE_COMPLEX)
259: a = PetscImaginaryPart(ar);
260: b = PetscImaginaryPart(br);
261: #else
262: a = PetscAbsReal(ai);
263: b = PetscAbsReal(bi);
264: #endif
265: if (a<b) *result = 1;
266: else if (a>b) *result = -1;
267: else { /* break the tie by checking the magnitude */
268: a = SlepcAbsEigenvalue(ar,ai);
269: b = SlepcAbsEigenvalue(br,bi);
270: if (a<b) *result = 1;
271: else if (a>b) *result = -1;
272: else *result = 0;
273: }
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: PetscErrorCode SlepcCompareSmallestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
278: {
279: PetscReal a,b;
281: PetscFunctionBegin;
282: #if defined(PETSC_USE_COMPLEX)
283: a = PetscImaginaryPart(ar);
284: b = PetscImaginaryPart(br);
285: #else
286: a = PetscAbsReal(ai);
287: b = PetscAbsReal(bi);
288: #endif
289: if (a>b) *result = 1;
290: else if (a<b) *result = -1;
291: else { /* break the tie by checking the magnitude */
292: a = SlepcAbsEigenvalue(ar,ai);
293: b = SlepcAbsEigenvalue(br,bi);
294: if (a<b) *result = 1;
295: else if (a>b) *result = -1;
296: else *result = 0;
297: }
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: PetscErrorCode SlepcCompareTargetMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
302: {
303: PetscReal a,b;
304: PetscScalar *target = (PetscScalar*)ctx;
306: PetscFunctionBegin;
307: /* complex target only allowed if scalartype=complex */
308: a = SlepcAbsEigenvalue(ar-(*target),ai);
309: b = SlepcAbsEigenvalue(br-(*target),bi);
310: if (a>b) *result = 1;
311: else if (a<b) *result = -1;
312: else *result = 0;
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: PetscErrorCode SlepcCompareTargetReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
317: {
318: PetscReal a,b;
319: PetscScalar *target = (PetscScalar*)ctx;
321: PetscFunctionBegin;
322: a = PetscAbsReal(PetscRealPart(ar-(*target)));
323: b = PetscAbsReal(PetscRealPart(br-(*target)));
324: if (a>b) *result = 1;
325: else if (a<b) *result = -1;
326: else *result = 0;
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: PetscErrorCode SlepcCompareTargetImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
331: {
332: #if defined(PETSC_USE_COMPLEX)
333: PetscReal a,b;
334: PetscScalar *target = (PetscScalar*)ctx;
335: #endif
337: PetscFunctionBegin;
338: #if defined(PETSC_USE_COMPLEX)
339: a = PetscAbsReal(PetscImaginaryPart(ar-(*target)));
340: b = PetscAbsReal(PetscImaginaryPart(br-(*target)));
341: if (a>b) *result = 1;
342: else if (a<b) *result = -1;
343: else *result = 0;
344: #else
345: *result = 0;
346: #endif
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: /*
351: Used in the SVD for computing smallest singular values
352: from the cyclic matrix.
353: */
354: PetscErrorCode SlepcCompareSmallestPosReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
355: {
356: PetscReal a,b;
357: PetscBool aisright,bisright;
359: PetscFunctionBegin;
360: if (PetscRealPart(ar)>0.0) aisright = PETSC_TRUE;
361: else aisright = PETSC_FALSE;
362: if (PetscRealPart(br)>0.0) bisright = PETSC_TRUE;
363: else bisright = PETSC_FALSE;
364: if (aisright == bisright) { /* same sign */
365: a = SlepcAbsEigenvalue(ar,ai);
366: b = SlepcAbsEigenvalue(br,bi);
367: if (a>b) *result = 1;
368: else if (a<b) *result = -1;
369: else *result = 0;
370: } else if (aisright && !bisright) *result = -1; /* 'a' is on the right */
371: else *result = 1; /* 'b' is on the right */
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }