Actual source code: slepcsc.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }