Actual source code: slepcsc.c
slepc-main 2024-11-09
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 *result = 0;
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: PetscErrorCode SlepcCompareSmallestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
223: {
224: PetscReal a,b;
226: PetscFunctionBegin;
227: #if defined(PETSC_USE_COMPLEX)
228: a = PetscImaginaryPart(ar);
229: b = PetscImaginaryPart(br);
230: #else
231: a = PetscAbsReal(ai);
232: b = PetscAbsReal(bi);
233: #endif
234: if (a>b) *result = 1;
235: else if (a<b) *result = -1;
236: else *result = 0;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: PetscErrorCode SlepcCompareTargetMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
241: {
242: PetscReal a,b;
243: PetscScalar *target = (PetscScalar*)ctx;
245: PetscFunctionBegin;
246: /* complex target only allowed if scalartype=complex */
247: a = SlepcAbsEigenvalue(ar-(*target),ai);
248: b = SlepcAbsEigenvalue(br-(*target),bi);
249: if (a>b) *result = 1;
250: else if (a<b) *result = -1;
251: else *result = 0;
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: PetscErrorCode SlepcCompareTargetReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
256: {
257: PetscReal a,b;
258: PetscScalar *target = (PetscScalar*)ctx;
260: PetscFunctionBegin;
261: a = PetscAbsReal(PetscRealPart(ar-(*target)));
262: b = PetscAbsReal(PetscRealPart(br-(*target)));
263: if (a>b) *result = 1;
264: else if (a<b) *result = -1;
265: else *result = 0;
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
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;
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
285: /*
286: Used in the SVD for computing smallest singular values
287: from the cyclic matrix.
288: */
289: PetscErrorCode SlepcCompareSmallestPosReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
290: {
291: PetscReal a,b;
292: PetscBool aisright,bisright;
294: PetscFunctionBegin;
295: if (PetscRealPart(ar)>0.0) aisright = PETSC_TRUE;
296: else aisright = PETSC_FALSE;
297: if (PetscRealPart(br)>0.0) bisright = PETSC_TRUE;
298: else bisright = PETSC_FALSE;
299: if (aisright == bisright) { /* same sign */
300: a = SlepcAbsEigenvalue(ar,ai);
301: b = SlepcAbsEigenvalue(br,bi);
302: if (a>b) *result = 1;
303: else if (a<b) *result = -1;
304: else *result = 0;
305: } else if (aisright && !bisright) *result = -1; /* 'a' is on the right */
306: else *result = 1; /* 'b' is on the right */
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }