Actual source code: slepcsc.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcsys.h" I*/
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: {
43: PetscScalar re[2],im[2];
44: PetscInt cin[2];
45: PetscBool inside[2];
49: #if defined(PETSC_USE_DEBUG)
50: if (!sc->comparison) SETERRQ(PETSC_COMM_SELF,1,"Undefined comparison function");
51: #endif
52: re[0] = ar; re[1] = br;
53: im[0] = ai; im[1] = bi;
54: if (sc->map) {
55: (*sc->map)(sc->mapobj,2,re,im);
56: }
57: if (sc->rg) {
58: RGCheckInside(sc->rg,2,re,im,cin);
59: inside[0] = PetscNot(cin[0]<0);
60: inside[1] = PetscNot(cin[1]<0);
61: if (inside[0] && !inside[1]) *res = -1;
62: else if (!inside[0] && inside[1]) *res = 1;
63: else {
64: (*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx);
65: }
66: } else {
67: (*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx);
68: }
69: return(0);
70: }
72: /*@
73: SlepcSortEigenvalues - Sorts a list of eigenvalues according to the
74: sorting criterion specified in a SlepcSC context.
76: Not Collective
78: Input Parameters:
79: + sc - the sorting criterion context
80: . n - number of eigenvalues in the list
81: . eigr - pointer to the array containing the eigenvalues
82: - eigi - imaginary part of the eigenvalues (only when using real numbers)
84: Input/Output Parameter:
85: . perm - permutation array. Must be initialized to 0:n-1 on input.
87: Note:
88: The result is a list of indices in the original eigenvalue array
89: corresponding to the first n eigenvalues sorted in the specified
90: criterion.
92: Level: developer
94: .seealso: SlepcSCCompare(), SlepcSC
95: @*/
96: PetscErrorCode SlepcSortEigenvalues(SlepcSC sc,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,PetscInt *perm)
97: {
99: PetscScalar re,im;
100: PetscInt i,j,result,tmp;
107: /* insertion sort */
108: for (i=n-1;i>=0;i--) {
109: re = eigr[perm[i]];
110: im = eigi[perm[i]];
111: j = i+1;
112: #if !defined(PETSC_USE_COMPLEX)
113: if (im!=0) {
114: /* complex eigenvalue */
115: i--;
116: im = eigi[perm[i]];
117: }
118: #endif
119: while (j<n) {
120: SlepcSCCompare(sc,re,im,eigr[perm[j]],eigi[perm[j]],&result);
121: if (result<0) break;
122: #if !defined(PETSC_USE_COMPLEX)
123: /* keep together every complex conjugated eigenpair */
124: if (!im) {
125: if (eigi[perm[j]] == 0.0) {
126: #endif
127: tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = tmp;
128: j++;
129: #if !defined(PETSC_USE_COMPLEX)
130: } else {
131: tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp;
132: j+=2;
133: }
134: } else {
135: if (eigi[perm[j]] == 0.0) {
136: tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = perm[j-1]; perm[j-1] = tmp;
137: j++;
138: } else {
139: tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = tmp;
140: tmp = perm[j-1]; perm[j-1] = perm[j+1]; perm[j+1] = tmp;
141: j+=2;
142: }
143: }
144: #endif
145: }
146: }
147: return(0);
148: }
150: /*
151: SlepcMap_ST - Gateway function to call STBackTransform from outside ST.
152: */
153: PetscErrorCode SlepcMap_ST(PetscObject obj,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
154: {
158: STBackTransform((ST)obj,n,eigr,eigi);
159: return(0);
160: }
162: PetscErrorCode SlepcCompareLargestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
163: {
164: PetscReal a,b;
167: a = SlepcAbsEigenvalue(ar,ai);
168: b = SlepcAbsEigenvalue(br,bi);
169: if (a<b) *result = 1;
170: else if (a>b) *result = -1;
171: else *result = 0;
172: return(0);
173: }
175: PetscErrorCode SlepcCompareSmallestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
176: {
177: PetscReal a,b;
180: a = SlepcAbsEigenvalue(ar,ai);
181: b = SlepcAbsEigenvalue(br,bi);
182: if (a>b) *result = 1;
183: else if (a<b) *result = -1;
184: else *result = 0;
185: return(0);
186: }
188: PetscErrorCode SlepcCompareLargestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
189: {
190: PetscReal a,b;
193: a = PetscRealPart(ar);
194: b = PetscRealPart(br);
195: if (a<b) *result = 1;
196: else if (a>b) *result = -1;
197: else *result = 0;
198: return(0);
199: }
201: PetscErrorCode SlepcCompareSmallestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
202: {
203: PetscReal a,b;
206: a = PetscRealPart(ar);
207: b = PetscRealPart(br);
208: if (a>b) *result = 1;
209: else if (a<b) *result = -1;
210: else *result = 0;
211: return(0);
212: }
214: PetscErrorCode SlepcCompareLargestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
215: {
216: PetscReal a,b;
219: #if defined(PETSC_USE_COMPLEX)
220: a = PetscImaginaryPart(ar);
221: b = PetscImaginaryPart(br);
222: #else
223: a = PetscAbsReal(ai);
224: b = PetscAbsReal(bi);
225: #endif
226: if (a<b) *result = 1;
227: else if (a>b) *result = -1;
228: else *result = 0;
229: return(0);
230: }
232: PetscErrorCode SlepcCompareSmallestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
233: {
234: PetscReal a,b;
237: #if defined(PETSC_USE_COMPLEX)
238: a = PetscImaginaryPart(ar);
239: b = PetscImaginaryPart(br);
240: #else
241: a = PetscAbsReal(ai);
242: b = PetscAbsReal(bi);
243: #endif
244: if (a>b) *result = 1;
245: else if (a<b) *result = -1;
246: else *result = 0;
247: return(0);
248: }
250: PetscErrorCode SlepcCompareTargetMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
251: {
252: PetscReal a,b;
253: PetscScalar *target = (PetscScalar*)ctx;
256: /* complex target only allowed if scalartype=complex */
257: a = SlepcAbsEigenvalue(ar-(*target),ai);
258: b = SlepcAbsEigenvalue(br-(*target),bi);
259: if (a>b) *result = 1;
260: else if (a<b) *result = -1;
261: else *result = 0;
262: return(0);
263: }
265: PetscErrorCode SlepcCompareTargetReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
266: {
267: PetscReal a,b;
268: PetscScalar *target = (PetscScalar*)ctx;
271: a = PetscAbsReal(PetscRealPart(ar-(*target)));
272: b = PetscAbsReal(PetscRealPart(br-(*target)));
273: if (a>b) *result = 1;
274: else if (a<b) *result = -1;
275: else *result = 0;
276: return(0);
277: }
279: #if defined(PETSC_USE_COMPLEX)
280: PetscErrorCode SlepcCompareTargetImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
281: {
282: PetscReal a,b;
283: PetscScalar *target = (PetscScalar*)ctx;
286: a = PetscAbsReal(PetscImaginaryPart(ar-(*target)));
287: b = PetscAbsReal(PetscImaginaryPart(br-(*target)));
288: if (a>b) *result = 1;
289: else if (a<b) *result = -1;
290: else *result = 0;
291: return(0);
292: }
293: #endif
295: /*
296: Used in the SVD for computing smallest singular values
297: from the cyclic matrix.
298: */
299: PetscErrorCode SlepcCompareSmallestPosReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
300: {
301: PetscReal a,b;
302: PetscBool aisright,bisright;
305: if (PetscRealPart(ar)>0.0) aisright = PETSC_TRUE;
306: else aisright = PETSC_FALSE;
307: if (PetscRealPart(br)>0.0) bisright = PETSC_TRUE;
308: else bisright = PETSC_FALSE;
309: if (aisright == bisright) { /* same sign */
310: a = SlepcAbsEigenvalue(ar,ai);
311: b = SlepcAbsEigenvalue(br,bi);
312: if (a>b) *result = 1;
313: else if (a<b) *result = -1;
314: else *result = 0;
315: } else if (aisright && !bisright) *result = -1; /* 'a' is on the right */
316: else *result = 1; /* 'b' is on the right */
317: return(0);
318: }