Actual source code: slepcsc.c
slepc-3.23.3 2025-09-08
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: /*@
66: SlepcSortEigenvalues - Sorts a list of eigenvalues according to the
67: sorting criterion specified in a SlepcSC context.
69: Not Collective
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)
77: Output Parameter:
78: . perm - permutation array. Must be initialized to 0:n-1 on input.
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.
85: Level: developer
87: .seealso: SlepcSCCompare(), SlepcSC
88: @*/
89: PetscErrorCode SlepcSortEigenvalues(SlepcSC sc,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,PetscInt *perm)
90: {
91: PetscScalar re,im;
92: PetscInt i,j,result,tmp;
94: PetscFunctionBegin;
95: PetscAssertPointer(sc,1);
96: PetscAssertPointer(eigr,3);
97: PetscAssertPointer(eigi,4);
98: PetscAssertPointer(perm,5);
99: /* insertion sort */
100: for (i=n-1;i>=0;i--) {
101: re = eigr[perm[i]];
102: im = eigi[perm[i]];
103: 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: while (j<n) {
112: PetscCall(SlepcSCCompare(sc,re,im,eigr[perm[j]],eigi[perm[j]],&result));
113: 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: tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = tmp;
120: 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: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: /*
143: SlepcMap_ST - Gateway function to call STBackTransform from outside ST.
144: */
145: PetscErrorCode SlepcMap_ST(PetscObject obj,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
146: {
147: PetscFunctionBegin;
148: PetscCall(STBackTransform((ST)obj,n,eigr,eigi));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: PetscErrorCode SlepcCompareLargestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
153: {
154: PetscReal a,b;
156: PetscFunctionBegin;
157: a = SlepcAbsEigenvalue(ar,ai);
158: b = SlepcAbsEigenvalue(br,bi);
159: if (a<b) *result = 1;
160: else if (a>b) *result = -1;
161: else *result = 0;
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: PetscErrorCode SlepcCompareSmallestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
166: {
167: PetscReal a,b;
169: PetscFunctionBegin;
170: a = SlepcAbsEigenvalue(ar,ai);
171: b = SlepcAbsEigenvalue(br,bi);
172: if (a>b) *result = 1;
173: else if (a<b) *result = -1;
174: else *result = 0;
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: PetscErrorCode SlepcCompareLargestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
179: {
180: PetscReal a,b;
182: PetscFunctionBegin;
183: a = PetscRealPart(ar);
184: b = PetscRealPart(br);
185: if (a<b) *result = 1;
186: else if (a>b) *result = -1;
187: else *result = 0;
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: PetscErrorCode SlepcCompareSmallestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
192: {
193: PetscReal a,b;
195: PetscFunctionBegin;
196: a = PetscRealPart(ar);
197: b = PetscRealPart(br);
198: if (a>b) *result = 1;
199: else if (a<b) *result = -1;
200: else *result = 0;
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: PetscErrorCode SlepcCompareLargestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
205: {
206: PetscReal a,b;
208: PetscFunctionBegin;
209: #if defined(PETSC_USE_COMPLEX)
210: a = PetscImaginaryPart(ar);
211: b = PetscImaginaryPart(br);
212: #else
213: a = PetscAbsReal(ai);
214: b = PetscAbsReal(bi);
215: #endif
216: if (a<b) *result = 1;
217: else if (a>b) *result = -1;
218: else { /* break the tie by checking the magnitude */
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: }
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: PetscErrorCode SlepcCompareSmallestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
229: {
230: PetscReal a,b;
232: PetscFunctionBegin;
233: #if defined(PETSC_USE_COMPLEX)
234: a = PetscImaginaryPart(ar);
235: b = PetscImaginaryPart(br);
236: #else
237: a = PetscAbsReal(ai);
238: b = PetscAbsReal(bi);
239: #endif
240: if (a>b) *result = 1;
241: else if (a<b) *result = -1;
242: else { /* break the tie by checking the magnitude */
243: a = SlepcAbsEigenvalue(ar,ai);
244: b = SlepcAbsEigenvalue(br,bi);
245: if (a<b) *result = 1;
246: else if (a>b) *result = -1;
247: else *result = 0;
248: }
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: PetscErrorCode SlepcCompareTargetMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
253: {
254: PetscReal a,b;
255: PetscScalar *target = (PetscScalar*)ctx;
257: PetscFunctionBegin;
258: /* complex target only allowed if scalartype=complex */
259: a = SlepcAbsEigenvalue(ar-(*target),ai);
260: b = SlepcAbsEigenvalue(br-(*target),bi);
261: if (a>b) *result = 1;
262: else if (a<b) *result = -1;
263: else *result = 0;
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: PetscErrorCode SlepcCompareTargetReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
268: {
269: PetscReal a,b;
270: PetscScalar *target = (PetscScalar*)ctx;
272: PetscFunctionBegin;
273: a = PetscAbsReal(PetscRealPart(ar-(*target)));
274: b = PetscAbsReal(PetscRealPart(br-(*target)));
275: if (a>b) *result = 1;
276: else if (a<b) *result = -1;
277: else *result = 0;
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: PetscErrorCode SlepcCompareTargetImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
282: {
283: #if defined(PETSC_USE_COMPLEX)
284: PetscReal a,b;
285: PetscScalar *target = (PetscScalar*)ctx;
286: #endif
288: PetscFunctionBegin;
289: #if defined(PETSC_USE_COMPLEX)
290: a = PetscAbsReal(PetscImaginaryPart(ar-(*target)));
291: b = PetscAbsReal(PetscImaginaryPart(br-(*target)));
292: if (a>b) *result = 1;
293: else if (a<b) *result = -1;
294: else *result = 0;
295: #else
296: *result = 0;
297: #endif
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: /*
302: Used in the SVD for computing smallest singular values
303: from the cyclic matrix.
304: */
305: PetscErrorCode SlepcCompareSmallestPosReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
306: {
307: PetscReal a,b;
308: PetscBool aisright,bisright;
310: PetscFunctionBegin;
311: if (PetscRealPart(ar)>0.0) aisright = PETSC_TRUE;
312: else aisright = PETSC_FALSE;
313: if (PetscRealPart(br)>0.0) bisright = PETSC_TRUE;
314: else bisright = PETSC_FALSE;
315: if (aisright == bisright) { /* same sign */
316: a = SlepcAbsEigenvalue(ar,ai);
317: b = SlepcAbsEigenvalue(br,bi);
318: if (a>b) *result = 1;
319: else if (a<b) *result = -1;
320: else *result = 0;
321: } else if (aisright && !bisright) *result = -1; /* 'a' is on the right */
322: else *result = 1; /* 'b' is on the right */
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }