Actual source code: slepcsc.c

slepc-3.21.1 2024-04-26
Report Typos and Errors
```  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.
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: }
```