LCOV - code coverage report
Current view: top level - sys - slepcsc.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 141 145 97.2 %
Date: 2024-03-28 00:52:16 Functions: 12 12 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : #include <slepc/private/slepcimpl.h>     /*I "slepcsys.h" I*/
      12             : #include <slepcrg.h>
      13             : #include <slepcst.h>
      14             : 
      15             : /*@
      16             :    SlepcSCCompare - Compares two (possibly complex) values according
      17             :    to a certain criterion.
      18             : 
      19             :    Not Collective
      20             : 
      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
      27             : 
      28             :    Output Parameter:
      29             : .  res - result of comparison
      30             : 
      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.
      35             : 
      36             :    Level: developer
      37             : 
      38             : .seealso: SlepcSortEigenvalues(), SlepcSC
      39             : @*/
      40     2457930 : PetscErrorCode SlepcSCCompare(SlepcSC sc,PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *res)
      41             : {
      42     2457930 :   PetscScalar    re[2],im[2];
      43     2457930 :   PetscInt       cin[2];
      44     2457930 :   PetscBool      inside[2];
      45             : 
      46     2457930 :   PetscFunctionBegin;
      47     2457930 :   PetscAssertPointer(res,6);
      48             : #if defined(PETSC_USE_DEBUG)
      49     2457930 :   PetscCheck(sc->comparison,PETSC_COMM_SELF,PETSC_ERR_USER,"Undefined comparison function");
      50             : #endif
      51     2457930 :   re[0] = ar; re[1] = br;
      52     2457930 :   im[0] = ai; im[1] = bi;
      53     2457930 :   if (sc->map) PetscCall((*sc->map)(sc->mapobj,2,re,im));
      54     2457930 :   if (sc->rg) {
      55       65421 :     PetscCall(RGCheckInside(sc->rg,2,re,im,cin));
      56       65421 :     inside[0] = PetscNot(cin[0]<0);
      57       65421 :     inside[1] = PetscNot(cin[1]<0);
      58       65421 :     if (inside[0] && !inside[1]) *res = -1;
      59       59382 :     else if (!inside[0] && inside[1]) *res = 1;
      60       58482 :     else PetscCall((*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx));
      61     2392509 :   } else PetscCall((*sc->comparison)(re[0],im[0],re[1],im[1],res,sc->comparisonctx));
      62     2457930 :   PetscFunctionReturn(PETSC_SUCCESS);
      63             : }
      64             : 
      65             : /*@
      66             :    SlepcSortEigenvalues - Sorts a list of eigenvalues according to the
      67             :    sorting criterion specified in a SlepcSC context.
      68             : 
      69             :    Not Collective
      70             : 
      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)
      76             : 
      77             :    Output Parameter:
      78             : .  perm - permutation array. Must be initialized to 0:n-1 on input.
      79             : 
      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.
      84             : 
      85             :    Level: developer
      86             : 
      87             : .seealso: SlepcSCCompare(), SlepcSC
      88             : @*/
      89        1148 : PetscErrorCode SlepcSortEigenvalues(SlepcSC sc,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,PetscInt *perm)
      90             : {
      91        1148 :   PetscScalar    re,im;
      92        1148 :   PetscInt       i,j,result,tmp;
      93             : 
      94        1148 :   PetscFunctionBegin;
      95        1148 :   PetscAssertPointer(sc,1);
      96        1148 :   PetscAssertPointer(eigr,3);
      97        1148 :   PetscAssertPointer(eigi,4);
      98        1148 :   PetscAssertPointer(perm,5);
      99             :   /* insertion sort */
     100       11492 :   for (i=n-1;i>=0;i--) {
     101       10344 :     re = eigr[perm[i]];
     102       10344 :     im = eigi[perm[i]];
     103       10344 :     j = i+1;
     104             : #if !defined(PETSC_USE_COMPLEX)
     105       10344 :     if (im!=0) {
     106             :       /* complex eigenvalue */
     107         237 :       i--;
     108         237 :       im = eigi[perm[i]];
     109             :     }
     110             : #endif
     111       10872 :     while (j<n) {
     112        9574 :       PetscCall(SlepcSCCompare(sc,re,im,eigr[perm[j]],eigi[perm[j]],&result));
     113        9574 :       if (result<=0) break;
     114             : #if !defined(PETSC_USE_COMPLEX)
     115             :       /* keep together every complex conjugated eigenpair */
     116         528 :       if (!im) {
     117         516 :         if (eigi[perm[j]] == 0.0) {
     118             : #endif
     119         516 :           tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = tmp;
     120         516 :           j++;
     121             : #if !defined(PETSC_USE_COMPLEX)
     122             :         } else {
     123           0 :           tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp;
     124           0 :           j+=2;
     125             :         }
     126             :       } else {
     127          12 :         if (eigi[perm[j]] == 0.0) {
     128           6 :           tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = perm[j-1]; perm[j-1] = tmp;
     129           6 :           j++;
     130             :         } else {
     131           6 :           tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = tmp;
     132           6 :           tmp = perm[j-1]; perm[j-1] = perm[j+1]; perm[j+1] = tmp;
     133           6 :           j+=2;
     134             :         }
     135             :       }
     136             : #endif
     137             :     }
     138             :   }
     139        1148 :   PetscFunctionReturn(PETSC_SUCCESS);
     140             : }
     141             : 
     142             : /*
     143             :    SlepcMap_ST - Gateway function to call STBackTransform from outside ST.
     144             : */
     145     2141178 : PetscErrorCode SlepcMap_ST(PetscObject obj,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
     146             : {
     147     2141178 :   PetscFunctionBegin;
     148     2141178 :   PetscCall(STBackTransform((ST)obj,n,eigr,eigi));
     149     2141178 :   PetscFunctionReturn(PETSC_SUCCESS);
     150             : }
     151             : 
     152      836660 : PetscErrorCode SlepcCompareLargestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
     153             : {
     154      836660 :   PetscReal a,b;
     155             : 
     156      836660 :   PetscFunctionBegin;
     157      836660 :   a = SlepcAbsEigenvalue(ar,ai);
     158      836660 :   b = SlepcAbsEigenvalue(br,bi);
     159      836660 :   if (a<b) *result = 1;
     160      710960 :   else if (a>b) *result = -1;
     161         474 :   else *result = 0;
     162      836660 :   PetscFunctionReturn(PETSC_SUCCESS);
     163             : }
     164             : 
     165       31622 : PetscErrorCode SlepcCompareSmallestMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
     166             : {
     167       31622 :   PetscReal a,b;
     168             : 
     169       31622 :   PetscFunctionBegin;
     170       31622 :   a = SlepcAbsEigenvalue(ar,ai);
     171       31622 :   b = SlepcAbsEigenvalue(br,bi);
     172       31622 :   if (a>b) *result = 1;
     173       15512 :   else if (a<b) *result = -1;
     174          30 :   else *result = 0;
     175       31622 :   PetscFunctionReturn(PETSC_SUCCESS);
     176             : }
     177             : 
     178      306713 : PetscErrorCode SlepcCompareLargestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
     179             : {
     180      306713 :   PetscReal a,b;
     181             : 
     182      306713 :   PetscFunctionBegin;
     183      306713 :   a = PetscRealPart(ar);
     184      306713 :   b = PetscRealPart(br);
     185      306713 :   if (a<b) *result = 1;
     186      252527 :   else if (a>b) *result = -1;
     187         779 :   else *result = 0;
     188      306713 :   PetscFunctionReturn(PETSC_SUCCESS);
     189             : }
     190             : 
     191      164217 : PetscErrorCode SlepcCompareSmallestReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
     192             : {
     193      164217 :   PetscReal a,b;
     194             : 
     195      164217 :   PetscFunctionBegin;
     196      164217 :   a = PetscRealPart(ar);
     197      164217 :   b = PetscRealPart(br);
     198      164217 :   if (a>b) *result = 1;
     199       83181 :   else if (a<b) *result = -1;
     200          69 :   else *result = 0;
     201      164217 :   PetscFunctionReturn(PETSC_SUCCESS);
     202             : }
     203             : 
     204        2345 : PetscErrorCode SlepcCompareLargestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
     205             : {
     206        2345 :   PetscReal a,b;
     207             : 
     208        2345 :   PetscFunctionBegin;
     209             : #if defined(PETSC_USE_COMPLEX)
     210             :   a = PetscImaginaryPart(ar);
     211             :   b = PetscImaginaryPart(br);
     212             : #else
     213        2345 :   a = PetscAbsReal(ai);
     214        2345 :   b = PetscAbsReal(bi);
     215             : #endif
     216        2345 :   if (a<b) *result = 1;
     217        2117 :   else if (a>b) *result = -1;
     218        1697 :   else *result = 0;
     219        2345 :   PetscFunctionReturn(PETSC_SUCCESS);
     220             : }
     221             : 
     222        5766 : PetscErrorCode SlepcCompareSmallestImaginary(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
     223             : {
     224        5766 :   PetscReal a,b;
     225             : 
     226        5766 :   PetscFunctionBegin;
     227             : #if defined(PETSC_USE_COMPLEX)
     228             :   a = PetscImaginaryPart(ar);
     229             :   b = PetscImaginaryPart(br);
     230             : #else
     231        5766 :   a = PetscAbsReal(ai);
     232        5766 :   b = PetscAbsReal(bi);
     233             : #endif
     234        5766 :   if (a>b) *result = 1;
     235        5184 :   else if (a<b) *result = -1;
     236        4544 :   else *result = 0;
     237        5766 :   PetscFunctionReturn(PETSC_SUCCESS);
     238             : }
     239             : 
     240     1003905 : PetscErrorCode SlepcCompareTargetMagnitude(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
     241             : {
     242     1003905 :   PetscReal   a,b;
     243     1003905 :   PetscScalar *target = (PetscScalar*)ctx;
     244             : 
     245     1003905 :   PetscFunctionBegin;
     246             :   /* complex target only allowed if scalartype=complex */
     247     1003905 :   a = SlepcAbsEigenvalue(ar-(*target),ai);
     248     1003905 :   b = SlepcAbsEigenvalue(br-(*target),bi);
     249     1003905 :   if (a>b) *result = 1;
     250      850017 :   else if (a<b) *result = -1;
     251           4 :   else *result = 0;
     252     1003905 :   PetscFunctionReturn(PETSC_SUCCESS);
     253             : }
     254             : 
     255        1532 : PetscErrorCode SlepcCompareTargetReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
     256             : {
     257        1532 :   PetscReal   a,b;
     258        1532 :   PetscScalar *target = (PetscScalar*)ctx;
     259             : 
     260        1532 :   PetscFunctionBegin;
     261        1532 :   a = PetscAbsReal(PetscRealPart(ar-(*target)));
     262        1532 :   b = PetscAbsReal(PetscRealPart(br-(*target)));
     263        1532 :   if (a>b) *result = 1;
     264         496 :   else if (a<b) *result = -1;
     265           0 :   else *result = 0;
     266        1532 :   PetscFunctionReturn(PETSC_SUCCESS);
     267             : }
     268             : 
     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;
     274             : 
     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
     284             : 
     285             : /*
     286             :    Used in the SVD for computing smallest singular values
     287             :    from the cyclic matrix.
     288             : */
     289        8531 : PetscErrorCode SlepcCompareSmallestPosReal(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result,void *ctx)
     290             : {
     291        8531 :   PetscReal a,b;
     292        8531 :   PetscBool aisright,bisright;
     293             : 
     294        8531 :   PetscFunctionBegin;
     295        8531 :   if (PetscRealPart(ar)>0.0) aisright = PETSC_TRUE;
     296        1415 :   else aisright = PETSC_FALSE;
     297        8531 :   if (PetscRealPart(br)>0.0) bisright = PETSC_TRUE;
     298        6163 :   else bisright = PETSC_FALSE;
     299        8531 :   if (aisright == bisright) { /* same sign */
     300        3513 :     a = SlepcAbsEigenvalue(ar,ai);
     301        3513 :     b = SlepcAbsEigenvalue(br,bi);
     302        3513 :     if (a>b) *result = 1;
     303         945 :     else if (a<b) *result = -1;
     304           0 :     else *result = 0;
     305        5018 :   } else if (aisright && !bisright) *result = -1; /* 'a' is on the right */
     306         135 :   else *result = 1;  /* 'b' is on the right */
     307        8531 :   PetscFunctionReturn(PETSC_SUCCESS);
     308             : }

Generated by: LCOV version 1.14