LCOV - code coverage report
Current view: top level - eps/interface - epsdefault.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 195 202 96.5 %
Date: 2024-11-21 00:34:55 Functions: 13 13 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             :    This file contains some simple default routines for common operations
      12             : */
      13             : 
      14             : #include <slepc/private/epsimpl.h>   /*I "slepceps.h" I*/
      15             : #include <slepcvec.h>
      16             : 
      17         618 : PetscErrorCode EPSBackTransform_Default(EPS eps)
      18             : {
      19         618 :   PetscFunctionBegin;
      20         618 :   PetscCall(STBackTransform(eps->st,eps->nconv,eps->eigr,eps->eigi));
      21         618 :   PetscFunctionReturn(PETSC_SUCCESS);
      22             : }
      23             : 
      24             : /*
      25             :   EPSComputeVectors_Hermitian - Copies the Lanczos vectors as eigenvectors
      26             :   using purification for generalized eigenproblems.
      27             :  */
      28         277 : PetscErrorCode EPSComputeVectors_Hermitian(EPS eps)
      29             : {
      30         277 :   PetscBool      iscayley,indef;
      31         277 :   Mat            B,C;
      32             : 
      33         277 :   PetscFunctionBegin;
      34         277 :   if (eps->purify) {
      35          71 :     PetscCall(EPS_Purify(eps,eps->nconv));
      36          71 :     PetscCall(BVNormalize(eps->V,NULL));
      37             :   } else {
      38             :     /* In the case of Cayley transform, eigenvectors need to be B-normalized */
      39         206 :     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley));
      40         206 :     if (iscayley && eps->isgeneralized) {
      41           1 :       PetscCall(STGetMatrix(eps->st,1,&B));
      42           1 :       PetscCall(BVGetMatrix(eps->V,&C,&indef));
      43           1 :       PetscCheck(!indef,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"The inner product should not be indefinite");
      44           1 :       PetscCall(BVSetMatrix(eps->V,B,PETSC_FALSE));
      45           1 :       PetscCall(BVNormalize(eps->V,NULL));
      46           1 :       PetscCall(BVSetMatrix(eps->V,C,PETSC_FALSE));  /* restore original matrix */
      47             :     }
      48             :   }
      49         277 :   PetscFunctionReturn(PETSC_SUCCESS);
      50             : }
      51             : 
      52             : /*
      53             :   EPSComputeVectors_Indefinite - similar to the Schur version but
      54             :   for indefinite problems
      55             :  */
      56          16 : PetscErrorCode EPSComputeVectors_Indefinite(EPS eps)
      57             : {
      58          16 :   PetscInt       n;
      59          16 :   Mat            X;
      60             : 
      61          16 :   PetscFunctionBegin;
      62          16 :   PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
      63          16 :   PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
      64          16 :   PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
      65          16 :   PetscCall(BVMultInPlace(eps->V,X,0,n));
      66          16 :   PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
      67             : 
      68             :   /* purification */
      69          16 :   if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));
      70             : 
      71             :   /* normalization */
      72          16 :   PetscCall(BVNormalize(eps->V,eps->eigi));
      73          16 :   PetscFunctionReturn(PETSC_SUCCESS);
      74             : }
      75             : 
      76             : /*
      77             :   EPSComputeVectors_Twosided - Adjust left eigenvectors in generalized problems: y = B^-* y.
      78             :  */
      79          18 : PetscErrorCode EPSComputeVectors_Twosided(EPS eps)
      80             : {
      81          18 :   PetscInt       i;
      82          18 :   Vec            w,y;
      83             : 
      84          18 :   PetscFunctionBegin;
      85          18 :   if (!eps->twosided || !eps->isgeneralized) PetscFunctionReturn(PETSC_SUCCESS);
      86           3 :   PetscCall(EPSSetWorkVecs(eps,1));
      87           3 :   w = eps->work[0];
      88          10 :   for (i=0;i<eps->nconv;i++) {
      89           7 :     PetscCall(BVCopyVec(eps->W,i,w));
      90           7 :     PetscCall(BVGetColumn(eps->W,i,&y));
      91           7 :     PetscCall(STMatSolveHermitianTranspose(eps->st,w,y));
      92           7 :     PetscCall(BVRestoreColumn(eps->W,i,&y));
      93             :   }
      94           3 :   PetscFunctionReturn(PETSC_SUCCESS);
      95             : }
      96             : 
      97             : /*
      98             :   EPSComputeVectors_Schur - Compute eigenvectors from the vectors
      99             :   provided by the eigensolver. This version is intended for solvers
     100             :   that provide Schur vectors. Given the partial Schur decomposition
     101             :   OP*V=V*T, the following steps are performed:
     102             :       1) compute eigenvectors of T: T*Z=Z*D
     103             :       2) compute eigenvectors of OP: X=V*Z
     104             :  */
     105         301 : PetscErrorCode EPSComputeVectors_Schur(EPS eps)
     106             : {
     107         301 :   PetscInt       i;
     108         301 :   Mat            Z;
     109         301 :   Vec            z;
     110             : 
     111         301 :   PetscFunctionBegin;
     112         301 :   if (eps->ishermitian) {
     113          49 :     if (eps->isgeneralized && !eps->ispositive) PetscCall(EPSComputeVectors_Indefinite(eps));
     114          49 :     else PetscCall(EPSComputeVectors_Hermitian(eps));
     115          49 :     PetscFunctionReturn(PETSC_SUCCESS);
     116             :   }
     117             : 
     118             :   /* right eigenvectors */
     119         252 :   PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
     120             : 
     121             :   /* V = V * Z */
     122         252 :   PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
     123         252 :   PetscCall(BVMultInPlace(eps->V,Z,0,eps->nconv));
     124         252 :   PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
     125             : 
     126             :   /* Purify eigenvectors */
     127         252 :   if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));
     128             : 
     129             :   /* Fix eigenvectors if balancing was used */
     130         252 :   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
     131          74 :     for (i=0;i<eps->nconv;i++) {
     132          62 :       PetscCall(BVGetColumn(eps->V,i,&z));
     133          62 :       PetscCall(VecPointwiseDivide(z,z,eps->D));
     134          62 :       PetscCall(BVRestoreColumn(eps->V,i,&z));
     135             :     }
     136             :   }
     137             : 
     138             :   /* normalize eigenvectors (when using purification or balancing) */
     139         252 :   if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) PetscCall(BVNormalize(eps->V,eps->eigi));
     140             : 
     141             :   /* left eigenvectors */
     142         252 :   if (eps->twosided) {
     143          13 :     PetscCall(DSVectors(eps->ds,DS_MAT_Y,NULL,NULL));
     144             :     /* W = W * Z */
     145          13 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Y,&Z));
     146          13 :     PetscCall(BVMultInPlace(eps->W,Z,0,eps->nconv));
     147          13 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Y,&Z));
     148             :     /* Fix left eigenvectors if balancing was used */
     149          13 :     if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
     150          10 :       for (i=0;i<eps->nconv;i++) {
     151           8 :         PetscCall(BVGetColumn(eps->W,i,&z));
     152           8 :         PetscCall(VecPointwiseMult(z,z,eps->D));
     153           8 :         PetscCall(BVRestoreColumn(eps->W,i,&z));
     154             :       }
     155             :     }
     156          13 :     PetscCall(EPSComputeVectors_Twosided(eps));
     157             :     /* normalize */
     158          13 :     PetscCall(BVNormalize(eps->W,eps->eigi));
     159             : #if !defined(PETSC_USE_COMPLEX)
     160          60 :     for (i=0;i<eps->nconv-1;i++) {
     161          47 :       if (eps->eigi[i] != 0.0) {
     162           2 :         if (eps->eigi[i] > 0.0) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
     163           2 :         i++;
     164             :       }
     165             :     }
     166             : #endif
     167             :   }
     168         252 :   PetscFunctionReturn(PETSC_SUCCESS);
     169             : }
     170             : 
     171             : /*@
     172             :    EPSSetWorkVecs - Sets a number of work vectors into an EPS object.
     173             : 
     174             :    Collective
     175             : 
     176             :    Input Parameters:
     177             : +  eps - eigensolver context
     178             : -  nw  - number of work vectors to allocate
     179             : 
     180             :    Developer Notes:
     181             :    This is SLEPC_EXTERN because it may be required by user plugin EPS
     182             :    implementations.
     183             : 
     184             :    Level: developer
     185             : 
     186             : .seealso: EPSSetUp()
     187             : @*/
     188        4050 : PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
     189             : {
     190        4050 :   Vec            t;
     191             : 
     192        4050 :   PetscFunctionBegin;
     193        4050 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     194       12150 :   PetscValidLogicalCollectiveInt(eps,nw,2);
     195        4050 :   PetscCheck(nw>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
     196        4050 :   if (eps->nwork < nw) {
     197         687 :     PetscCall(VecDestroyVecs(eps->nwork,&eps->work));
     198         687 :     eps->nwork = nw;
     199         687 :     PetscCall(BVGetColumn(eps->V,0,&t));
     200         687 :     PetscCall(VecDuplicateVecs(t,nw,&eps->work));
     201         687 :     PetscCall(BVRestoreColumn(eps->V,0,&t));
     202             :   }
     203        4050 :   PetscFunctionReturn(PETSC_SUCCESS);
     204             : }
     205             : 
     206             : /*
     207             :   EPSSetWhichEigenpairs_Default - Sets the default value for which,
     208             :   depending on the ST.
     209             :  */
     210         114 : PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
     211             : {
     212         114 :   PetscBool      target;
     213             : 
     214         114 :   PetscFunctionBegin;
     215         114 :   PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,""));
     216         114 :   if (target) eps->which = EPS_TARGET_MAGNITUDE;
     217          95 :   else eps->which = EPS_LARGEST_MAGNITUDE;
     218         114 :   PetscFunctionReturn(PETSC_SUCCESS);
     219             : }
     220             : 
     221             : /*
     222             :   EPSConvergedRelative - Checks convergence relative to the eigenvalue.
     223             : */
     224       46496 : PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
     225             : {
     226       46496 :   PetscReal w;
     227             : 
     228       46496 :   PetscFunctionBegin;
     229       46496 :   w = SlepcAbsEigenvalue(eigr,eigi);
     230       46496 :   *errest = (w!=0.0)? res/w: PETSC_MAX_REAL;
     231       46496 :   PetscFunctionReturn(PETSC_SUCCESS);
     232             : }
     233             : 
     234             : /*
     235             :   EPSConvergedAbsolute - Checks convergence absolutely.
     236             : */
     237        4151 : PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
     238             : {
     239        4151 :   PetscFunctionBegin;
     240        4151 :   *errest = res;
     241        4151 :   PetscFunctionReturn(PETSC_SUCCESS);
     242             : }
     243             : 
     244             : /*
     245             :   EPSConvergedNorm - Checks convergence relative to the eigenvalue and
     246             :   the matrix norms.
     247             : */
     248        5026 : PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
     249             : {
     250        5026 :   PetscReal w;
     251             : 
     252        5026 :   PetscFunctionBegin;
     253        5026 :   w = SlepcAbsEigenvalue(eigr,eigi);
     254        5026 :   *errest = res / (eps->nrma + w*eps->nrmb);
     255        5026 :   PetscFunctionReturn(PETSC_SUCCESS);
     256             : }
     257             : 
     258             : /*@C
     259             :    EPSStoppingBasic - Default routine to determine whether the outer eigensolver
     260             :    iteration must be stopped.
     261             : 
     262             :    Collective
     263             : 
     264             :    Input Parameters:
     265             : +  eps    - eigensolver context obtained from EPSCreate()
     266             : .  its    - current number of iterations
     267             : .  max_it - maximum number of iterations
     268             : .  nconv  - number of currently converged eigenpairs
     269             : .  nev    - number of requested eigenpairs
     270             : -  ctx    - context (not used here)
     271             : 
     272             :    Output Parameter:
     273             : .  reason - result of the stopping test
     274             : 
     275             :    Notes:
     276             :    A positive value of reason indicates that the iteration has finished successfully
     277             :    (converged), and a negative value indicates an error condition (diverged). If
     278             :    the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
     279             :    (zero).
     280             : 
     281             :    EPSStoppingBasic() will stop if all requested eigenvalues are converged, or if
     282             :    the maximum number of iterations has been reached.
     283             : 
     284             :    Use EPSSetStoppingTest() to provide your own test instead of using this one.
     285             : 
     286             :    Level: advanced
     287             : 
     288             : .seealso: EPSSetStoppingTest(), EPSConvergedReason, EPSGetConvergedReason()
     289             : @*/
     290       46773 : PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
     291             : {
     292       46773 :   PetscFunctionBegin;
     293       46773 :   *reason = EPS_CONVERGED_ITERATING;
     294       46773 :   if (nconv >= nev) {
     295         736 :     PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
     296         736 :     *reason = EPS_CONVERGED_TOL;
     297       46037 :   } else if (its >= max_it) {
     298          14 :     *reason = EPS_DIVERGED_ITS;
     299          14 :     PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
     300             :   }
     301       46773 :   PetscFunctionReturn(PETSC_SUCCESS);
     302             : }
     303             : 
     304             : /*
     305             :   EPSComputeRitzVector - Computes the current Ritz vector.
     306             : 
     307             :   Simple case (complex scalars or real scalars with Zi=NULL):
     308             :     x = V*Zr  (V is a basis of nv vectors, Zr has length nv)
     309             : 
     310             :   Split case:
     311             :     x = V*Zr  y = V*Zi  (Zr and Zi have length nv)
     312             : */
     313         385 : PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
     314             : {
     315         385 :   PetscInt       l,k;
     316         385 :   PetscReal      norm;
     317             : #if !defined(PETSC_USE_COMPLEX)
     318         385 :   Vec            z;
     319             : #endif
     320             : 
     321         385 :   PetscFunctionBegin;
     322             :   /* compute eigenvector */
     323         385 :   PetscCall(BVGetActiveColumns(V,&l,&k));
     324         385 :   PetscCall(BVSetActiveColumns(V,0,k));
     325         385 :   PetscCall(BVMultVec(V,1.0,0.0,x,Zr));
     326             : 
     327             :   /* purify eigenvector if necessary */
     328         385 :   if (eps->purify) {
     329          30 :     PetscCall(STApply(eps->st,x,y));
     330          30 :     if (eps->ishermitian) PetscCall(BVNormVec(eps->V,y,NORM_2,&norm));
     331          18 :     else PetscCall(VecNorm(y,NORM_2,&norm));
     332          30 :     PetscCall(VecScale(y,1.0/norm));
     333          30 :     PetscCall(VecCopy(y,x));
     334             :   }
     335             :   /* fix eigenvector if balancing is used */
     336         385 :   if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(x,x,eps->D));
     337             : #if !defined(PETSC_USE_COMPLEX)
     338             :   /* compute imaginary part of eigenvector */
     339         385 :   if (Zi) {
     340          25 :     PetscCall(BVMultVec(V,1.0,0.0,y,Zi));
     341          25 :     if (eps->ispositive) {
     342           0 :       PetscCall(BVCreateVec(V,&z));
     343           0 :       PetscCall(STApply(eps->st,y,z));
     344           0 :       PetscCall(VecNorm(z,NORM_2,&norm));
     345           0 :       PetscCall(VecScale(z,1.0/norm));
     346           0 :       PetscCall(VecCopy(z,y));
     347           0 :       PetscCall(VecDestroy(&z));
     348             :     }
     349          25 :     if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(y,y,eps->D));
     350             :   } else
     351             : #endif
     352         360 :     PetscCall(VecSet(y,0.0));
     353             : 
     354             :   /* normalize eigenvectors (when using balancing) */
     355         385 :   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
     356             : #if !defined(PETSC_USE_COMPLEX)
     357          12 :     if (Zi) PetscCall(VecNormalizeComplex(x,y,PETSC_TRUE,NULL));
     358             :     else
     359             : #endif
     360           0 :     PetscCall(VecNormalize(x,NULL));
     361             :   }
     362         385 :   PetscCall(BVSetActiveColumns(V,l,k));
     363         385 :   PetscFunctionReturn(PETSC_SUCCESS);
     364             : }
     365             : 
     366             : /*
     367             :   EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
     368             :   diagonal matrix to be applied for balancing in non-Hermitian problems.
     369             : */
     370          13 : PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
     371             : {
     372          13 :   Vec               z,p,r;
     373          13 :   PetscInt          i,j;
     374          13 :   PetscReal         norma;
     375          13 :   PetscScalar       *pz,*pD;
     376          13 :   const PetscScalar *pr,*pp;
     377          13 :   PetscRandom       rand;
     378             : 
     379          13 :   PetscFunctionBegin;
     380          13 :   PetscCall(EPSSetWorkVecs(eps,3));
     381          13 :   PetscCall(BVGetRandomContext(eps->V,&rand));
     382          13 :   r = eps->work[0];
     383          13 :   p = eps->work[1];
     384          13 :   z = eps->work[2];
     385          13 :   PetscCall(VecSet(eps->D,1.0));
     386             : 
     387          78 :   for (j=0;j<eps->balance_its;j++) {
     388             : 
     389             :     /* Build a random vector of +-1's */
     390          65 :     PetscCall(VecSetRandom(z,rand));
     391          65 :     PetscCall(VecGetArray(z,&pz));
     392        5095 :     for (i=0;i<eps->nloc;i++) {
     393        5030 :       if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
     394        2394 :       else pz[i]=1.0;
     395             :     }
     396          65 :     PetscCall(VecRestoreArray(z,&pz));
     397             : 
     398             :     /* Compute p=DA(D\z) */
     399          65 :     PetscCall(VecPointwiseDivide(r,z,eps->D));
     400          65 :     PetscCall(STApply(eps->st,r,p));
     401          65 :     PetscCall(VecPointwiseMult(p,p,eps->D));
     402          65 :     if (eps->balance == EPS_BALANCE_TWOSIDE) {
     403          40 :       if (j==0) {
     404             :         /* Estimate the matrix inf-norm */
     405           8 :         PetscCall(VecAbs(p));
     406           8 :         PetscCall(VecMax(p,NULL,&norma));
     407             :       }
     408             :       /* Compute r=D\(A'Dz) */
     409          40 :       PetscCall(VecPointwiseMult(z,z,eps->D));
     410          40 :       PetscCall(STApplyHermitianTranspose(eps->st,z,r));
     411          40 :       PetscCall(VecPointwiseDivide(r,r,eps->D));
     412             :     }
     413             : 
     414             :     /* Adjust values of D */
     415          65 :     PetscCall(VecGetArrayRead(r,&pr));
     416          65 :     PetscCall(VecGetArrayRead(p,&pp));
     417          65 :     PetscCall(VecGetArray(eps->D,&pD));
     418        5095 :     for (i=0;i<eps->nloc;i++) {
     419        5030 :       if (eps->balance == EPS_BALANCE_TWOSIDE) {
     420        2900 :         if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
     421        2814 :           pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
     422             :       } else {
     423        2130 :         if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
     424             :       }
     425             :     }
     426          65 :     PetscCall(VecRestoreArrayRead(r,&pr));
     427          65 :     PetscCall(VecRestoreArrayRead(p,&pp));
     428          65 :     PetscCall(VecRestoreArray(eps->D,&pD));
     429             :   }
     430          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     431             : }

Generated by: LCOV version 1.14