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 PetscDefined(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 !PetscDefined(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 !PetscDefined(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 !PetscDefined(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 scalars)
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 scalars)
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: /*@
202: SlepcCompareLargestMagnitude - An eigenvalue comparison function used to sort with respect
203: to largest magnitude.
205: Logically Collective
207: Input Parameters:
208: + ar - real part of the 1st eigenvalue
209: . ai - imaginary part of the 1st eigenvalue
210: . br - real part of the 2nd eigenvalue
211: . bi - imaginary part of the 2nd eigenvalue
212: - ctx - user-defined context, not used here
214: Output Parameter:
215: . result - result of comparison
217: Note:
218: The result is 1 if $|\lambda_1|<|\lambda_2|$.
220: Level: developer
222: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
223: @*/
224: PetscErrorCode SlepcCompareLargestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
225: {
226: PetscReal a,b;
228: PetscFunctionBegin;
229: a = SlepcAbsEigenvalue(ar,ai);
230: b = SlepcAbsEigenvalue(br,bi);
231: if (a<b) *result = 1;
232: else if (a>b) *result = -1;
233: else *result = 0;
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: /*@
238: SlepcCompareSmallestMagnitude - An eigenvalue comparison function used to sort with respect
239: to smallest magnitude.
241: Logically Collective
243: Input Parameters:
244: + ar - real part of the 1st eigenvalue
245: . ai - imaginary part of the 1st eigenvalue
246: . br - real part of the 2nd eigenvalue
247: . bi - imaginary part of the 2nd eigenvalue
248: - ctx - user-defined context, not used here
250: Output Parameter:
251: . result - result of comparison
253: Note:
254: The result is 1 if $|\lambda_1|>|\lambda_2|$.
256: Level: developer
258: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
259: @*/
260: PetscErrorCode SlepcCompareSmallestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
261: {
262: PetscReal a,b;
264: PetscFunctionBegin;
265: a = SlepcAbsEigenvalue(ar,ai);
266: b = SlepcAbsEigenvalue(br,bi);
267: if (a>b) *result = 1;
268: else if (a<b) *result = -1;
269: else *result = 0;
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*@
274: SlepcCompareLargestReal - An eigenvalue comparison function used to sort with respect
275: to largest real part.
277: Logically Collective
279: Input Parameters:
280: + ar - real part of the 1st eigenvalue
281: . ai - imaginary part of the 1st eigenvalue
282: . br - real part of the 2nd eigenvalue
283: . bi - imaginary part of the 2nd eigenvalue
284: - ctx - user-defined context, not used here
286: Output Parameter:
287: . result - result of comparison
289: Note:
290: The result is 1 if $\mathrm{Re}(\lambda_1)<\mathrm{Re}(\lambda_2)$.
292: Level: developer
294: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
295: @*/
296: PetscErrorCode SlepcCompareLargestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
297: {
298: PetscReal a,b;
300: PetscFunctionBegin;
301: a = PetscRealPart(ar);
302: b = PetscRealPart(br);
303: if (a<b) *result = 1;
304: else if (a>b) *result = -1;
305: else *result = 0;
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*@
310: SlepcCompareSmallestReal - An eigenvalue comparison function used to sort with respect
311: to smallest real part.
313: Logically Collective
315: Input Parameters:
316: + ar - real part of the 1st eigenvalue
317: . ai - imaginary part of the 1st eigenvalue
318: . br - real part of the 2nd eigenvalue
319: . bi - imaginary part of the 2nd eigenvalue
320: - ctx - user-defined context, not used here
322: Output Parameter:
323: . result - result of comparison
325: Note:
326: The result is 1 if $\mathrm{Re}(\lambda_1)>\mathrm{Re}(\lambda_2)$.
328: Level: developer
330: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
331: @*/
332: PetscErrorCode SlepcCompareSmallestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
333: {
334: PetscReal a,b;
336: PetscFunctionBegin;
337: a = PetscRealPart(ar);
338: b = PetscRealPart(br);
339: if (a>b) *result = 1;
340: else if (a<b) *result = -1;
341: else *result = 0;
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /*@
346: SlepcCompareLargestImaginary - An eigenvalue comparison function used to sort with respect
347: to largest imaginary real part.
349: Logically Collective
351: Input Parameters:
352: + ar - real part of the 1st eigenvalue
353: . ai - imaginary part of the 1st eigenvalue
354: . br - real part of the 2nd eigenvalue
355: . bi - imaginary part of the 2nd eigenvalue
356: - ctx - user-defined context, not used here
358: Output Parameter:
359: . result - result of comparison
361: Note:
362: In complex scalars, the result is 1 if $\mathrm{Im}(\lambda_1)<\mathrm{Im}(\lambda_2)$.
363: In real scalars, the result is 1 if $|\mathrm{Im}(\lambda_1)|<|\mathrm{Im}(\lambda_2)|$.
364: If the two values are equal, then $|\lambda_1|<|\lambda_2|$ is used to break the tie.
366: Level: developer
368: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
369: @*/
370: PetscErrorCode SlepcCompareLargestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
371: {
372: PetscReal a,b;
374: PetscFunctionBegin;
375: #if PetscDefined(USE_COMPLEX)
376: a = PetscImaginaryPart(ar);
377: b = PetscImaginaryPart(br);
378: #else
379: a = PetscAbsReal(ai);
380: b = PetscAbsReal(bi);
381: #endif
382: if (a<b) *result = 1;
383: else if (a>b) *result = -1;
384: else { /* break the tie by checking the magnitude */
385: a = SlepcAbsEigenvalue(ar,ai);
386: b = SlepcAbsEigenvalue(br,bi);
387: if (a<b) *result = 1;
388: else if (a>b) *result = -1;
389: else *result = 0;
390: }
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /*@
395: SlepcCompareSmallestImaginary - An eigenvalue comparison function used to sort with respect
396: to smallest imaginary real part.
398: Logically Collective
400: Input Parameters:
401: + ar - real part of the 1st eigenvalue
402: . ai - imaginary part of the 1st eigenvalue
403: . br - real part of the 2nd eigenvalue
404: . bi - imaginary part of the 2nd eigenvalue
405: - ctx - user-defined context, not used here
407: Output Parameter:
408: . result - result of comparison
410: Notes:
411: In complex scalars, the result is 1 if $\mathrm{Im}(\lambda_1)>\mathrm{Im}(\lambda_2)$.
412: In real scalars, the result is 1 if $|\mathrm{Im}(\lambda_1)|>|\mathrm{Im}(\lambda_2)|$.
413: If the two values are equal, then $|\lambda_1|<|\lambda_2|$ is used to break the tie.
415: Level: developer
417: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
418: @*/
419: PetscErrorCode SlepcCompareSmallestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
420: {
421: PetscReal a,b;
423: PetscFunctionBegin;
424: #if PetscDefined(USE_COMPLEX)
425: a = PetscImaginaryPart(ar);
426: b = PetscImaginaryPart(br);
427: #else
428: a = PetscAbsReal(ai);
429: b = PetscAbsReal(bi);
430: #endif
431: if (a>b) *result = 1;
432: else if (a<b) *result = -1;
433: else { /* break the tie by checking the magnitude */
434: a = SlepcAbsEigenvalue(ar,ai);
435: b = SlepcAbsEigenvalue(br,bi);
436: if (a<b) *result = 1;
437: else if (a>b) *result = -1;
438: else *result = 0;
439: }
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: /*@
444: SlepcCompareTargetMagnitude - An eigenvalue comparison function used to sort with respect
445: to a target (in magnitude).
447: Logically Collective
449: Input Parameters:
450: + ar - real part of the 1st eigenvalue
451: . ai - imaginary part of the 1st eigenvalue
452: . br - real part of the 2nd eigenvalue
453: . bi - imaginary part of the 2nd eigenvalue
454: - ctx - user-defined context, contains the target value (a `PetscScalar`)
456: Output Parameter:
457: . result - result of comparison
459: Note:
460: The result is 1 if $|\lambda_1-\tau|<|\lambda_2-\tau|$ where $\tau$ is the target.
462: Level: developer
464: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
465: @*/
466: PetscErrorCode SlepcCompareTargetMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
467: {
468: PetscReal a,b;
469: PetscScalar *target = (PetscScalar*)ctx;
471: PetscFunctionBegin;
472: /* complex target only allowed if scalartype=complex */
473: a = SlepcAbsEigenvalue(ar-(*target),ai);
474: b = SlepcAbsEigenvalue(br-(*target),bi);
475: if (a>b) *result = 1;
476: else if (a<b) *result = -1;
477: else *result = 0;
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@
482: SlepcCompareTargetReal - An eigenvalue comparison function used to sort with respect
483: to a target (along the real axis).
485: Logically Collective
487: Input Parameters:
488: + ar - real part of the 1st eigenvalue
489: . ai - imaginary part of the 1st eigenvalue
490: . br - real part of the 2nd eigenvalue
491: . bi - imaginary part of the 2nd eigenvalue
492: - ctx - user-defined context, contains the target value (a `PetscScalar`)
494: Output Parameter:
495: . result - result of comparison
497: Note:
498: The result is 1 if $\mathrm{Re}(\lambda_1-\tau)<\mathrm{Re}(\lambda_2-\tau)$ where
499: $\tau$ is the target.
501: Level: developer
503: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
504: @*/
505: PetscErrorCode SlepcCompareTargetReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
506: {
507: PetscReal a,b;
508: PetscScalar *target = (PetscScalar*)ctx;
510: PetscFunctionBegin;
511: a = PetscAbsReal(PetscRealPart(ar-(*target)));
512: b = PetscAbsReal(PetscRealPart(br-(*target)));
513: if (a>b) *result = 1;
514: else if (a<b) *result = -1;
515: else *result = 0;
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: /*@
520: SlepcCompareTargetImaginary - An eigenvalue comparison function used to sort with respect
521: to a target (along the imaginary axis).
523: Logically Collective
525: Input Parameters:
526: + ar - real part of the 1st eigenvalue
527: . ai - imaginary part of the 1st eigenvalue
528: . br - real part of the 2nd eigenvalue
529: . bi - imaginary part of the 2nd eigenvalue
530: - ctx - user-defined context, contains the target value (a `PetscScalar`)
532: Output Parameter:
533: . result - result of comparison
535: Notes:
536: In complex scalars, the result is 1 if $\mathrm{Im}(\lambda_1-\tau)<\mathrm{Im}(\lambda_2-\tau)$
537: where $\tau$ is the target.
538: In real scalars, the result is always zero because sorting with respect to target is
539: not supported (the target is always real in this case).
541: Level: developer
543: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
544: @*/
545: PetscErrorCode SlepcCompareTargetImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
546: {
547: #if PetscDefined(USE_COMPLEX)
548: PetscReal a,b;
549: PetscScalar *target = (PetscScalar*)ctx;
550: #endif
552: PetscFunctionBegin;
553: #if PetscDefined(USE_COMPLEX)
554: a = PetscAbsReal(PetscImaginaryPart(ar-(*target)));
555: b = PetscAbsReal(PetscImaginaryPart(br-(*target)));
556: if (a>b) *result = 1;
557: else if (a<b) *result = -1;
558: else *result = 0;
559: #else
560: *result = 0;
561: #endif
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: /*@
566: SlepcCompareSmallestPosReal - An eigenvalue comparison function used to sort according
567: to the smallest positive real part.
569: Logically Collective
571: Input Parameters:
572: + ar - real part of the 1st eigenvalue
573: . ai - imaginary part of the 1st eigenvalue
574: . br - real part of the 2nd eigenvalue
575: . bi - imaginary part of the 2nd eigenvalue
576: - ctx - user-defined context, not used here
578: Output Parameter:
579: . result - result of comparison
581: Note:
582: This sorting criterion is used in the SVD for computing smallest singular values
583: from the cyclic matrix. In that case, the computed values come in pairs $\pm\lambda$.
584: If the two values to compare have the same sign, the preferred one is the one with
585: smallest magnitude. Otherwise, the prefered one is the rightmost (positive sign).
587: Level: developer
589: .seealso: `SlepcEigenvalueComparisonFn`, `SlepcSCCompare()`, `SlepcSC`
590: @*/
591: PetscErrorCode SlepcCompareSmallestPosReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,PetscCtx ctx)
592: {
593: PetscReal a,b;
594: PetscBool aisright,bisright;
596: PetscFunctionBegin;
597: if (PetscRealPart(ar)>0.0) aisright = PETSC_TRUE;
598: else aisright = PETSC_FALSE;
599: if (PetscRealPart(br)>0.0) bisright = PETSC_TRUE;
600: else bisright = PETSC_FALSE;
601: if (aisright == bisright) { /* same sign */
602: a = SlepcAbsEigenvalue(ar,ai);
603: b = SlepcAbsEigenvalue(br,bi);
604: if (a>b) *result = 1;
605: else if (a<b) *result = -1;
606: else *result = 0;
607: } else if (aisright && !bisright) *result = -1; /* 'a' is on the right */
608: else *result = 1; /* 'b' is on the right */
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }