LCOV - code coverage report
Current view: top level - eps/impls/ciss - ciss.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 740 802 92.3 %
Date: 2024-03-28 00:28:38 Functions: 38 39 97.4 %
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             :    SLEPc eigensolver: "ciss"
      12             : 
      13             :    Method: Contour Integral Spectral Slicing
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Contour integral based on Sakurai-Sugiura method to construct a
      18             :        subspace, with various eigenpair extractions (Rayleigh-Ritz,
      19             :        explicit moment).
      20             : 
      21             :    Based on code contributed by Y. Maeda, T. Sakurai.
      22             : 
      23             :    References:
      24             : 
      25             :        [1] T. Sakurai and H. Sugiura, "A projection method for generalized
      26             :            eigenvalue problems", J. Comput. Appl. Math. 159:119-128, 2003.
      27             : 
      28             :        [2] T. Sakurai and H. Tadano, "CIRR: a Rayleigh-Ritz type method with
      29             :            contour integral for generalized eigenvalue problems", Hokkaido
      30             :            Math. J. 36:745-757, 2007.
      31             : */
      32             : 
      33             : #include <slepc/private/epsimpl.h>                /*I "slepceps.h" I*/
      34             : #include <slepc/private/slepccontour.h>
      35             : #include <slepcblaslapack.h>
      36             : 
      37             : typedef struct {
      38             :   /* user parameters */
      39             :   PetscInt          N;          /* number of integration points (32) */
      40             :   PetscInt          L;          /* block size (16) */
      41             :   PetscInt          M;          /* moment degree (N/4 = 4) */
      42             :   PetscReal         delta;      /* threshold of singular value (1e-12) */
      43             :   PetscInt          L_max;      /* maximum number of columns of the source matrix V */
      44             :   PetscReal         spurious_threshold; /* discard spurious eigenpairs */
      45             :   PetscBool         isreal;     /* A and B are real */
      46             :   PetscInt          npart;      /* number of partitions */
      47             :   PetscInt          refine_inner;
      48             :   PetscInt          refine_blocksize;
      49             :   EPSCISSQuadRule   quad;
      50             :   EPSCISSExtraction extraction;
      51             :   PetscBool         usest;
      52             :   /* private data */
      53             :   SlepcContourData  contour;
      54             :   PetscReal         *sigma;     /* threshold for numerical rank */
      55             :   PetscScalar       *weight;
      56             :   PetscScalar       *omega;
      57             :   PetscScalar       *pp;
      58             :   BV                V;
      59             :   BV                S;
      60             :   BV                pV;
      61             :   BV                Y;
      62             :   PetscBool         useconj;
      63             :   PetscBool         usest_set;  /* whether the user set the usest flag or not */
      64             :   PetscObjectId     rgid;
      65             :   PetscObjectState  rgstate;
      66             : } EPS_CISS;
      67             : 
      68             : /*
      69             :   Set up KSP solvers for every integration point, only called if !ctx->usest
      70             : */
      71          18 : static PetscErrorCode EPSCISSSetUp(EPS eps,Mat A,Mat B,Mat Pa,Mat Pb)
      72             : {
      73          18 :   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
      74          18 :   SlepcContourData contour;
      75          18 :   PetscInt         i,p_id,nsplit;
      76          18 :   Mat              Amat,Pmat;
      77          18 :   MatStructure     str,strp;
      78             : 
      79          18 :   PetscFunctionBegin;
      80          18 :   if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
      81          18 :   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
      82          18 :   contour = ctx->contour;
      83          18 :   PetscCall(STGetMatStructure(eps->st,&str));
      84          18 :   PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,&strp));
      85         386 :   for (i=0;i<contour->npoints;i++) {
      86         368 :     p_id = i*contour->subcomm->n + contour->subcomm->color;
      87         368 :     PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&Amat));
      88         368 :     if (B) PetscCall(MatAXPY(Amat,-ctx->omega[p_id],B,str));
      89         228 :     else PetscCall(MatShift(Amat,-ctx->omega[p_id]));
      90         368 :     if (nsplit) {
      91          64 :       PetscCall(MatDuplicate(Pa,MAT_COPY_VALUES,&Pmat));
      92          64 :       if (Pb) PetscCall(MatAXPY(Pmat,-ctx->omega[p_id],Pb,strp));
      93           0 :       else PetscCall(MatShift(Pmat,-ctx->omega[p_id]));
      94         304 :     } else Pmat = Amat;
      95         368 :     PetscCall(EPS_KSPSetOperators(contour->ksp[i],Amat,Amat));
      96         368 :     PetscCall(MatDestroy(&Amat));
      97         368 :     if (nsplit) PetscCall(MatDestroy(&Pmat));
      98             :   }
      99          18 :   PetscFunctionReturn(PETSC_SUCCESS);
     100             : }
     101             : 
     102             : /*
     103             :   Y_i = (A-z_i B)^{-1}BV for every integration point, Y=[Y_i] is in the context
     104             : */
     105          37 : static PetscErrorCode EPSCISSSolve(EPS eps,Mat B,BV V,PetscInt L_start,PetscInt L_end)
     106             : {
     107          37 :   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
     108          37 :   SlepcContourData contour;
     109          37 :   PetscInt         i,p_id;
     110          37 :   Mat              MV,BMV=NULL,MC;
     111          37 :   KSP              ksp;
     112             : 
     113          37 :   PetscFunctionBegin;
     114          37 :   if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
     115          37 :   contour = ctx->contour;
     116          37 :   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
     117          37 :   PetscCall(BVSetActiveColumns(V,L_start,L_end));
     118          37 :   PetscCall(BVGetMat(V,&MV));
     119         941 :   for (i=0;i<contour->npoints;i++) {
     120         904 :     p_id = i*contour->subcomm->n + contour->subcomm->color;
     121         904 :     if (ctx->usest)  {
     122         504 :       PetscCall(STSetShift(eps->st,ctx->omega[p_id]));
     123         504 :       PetscCall(STGetKSP(eps->st,&ksp));
     124         400 :     } else ksp = contour->ksp[i];
     125         904 :     PetscCall(BVSetActiveColumns(ctx->Y,i*ctx->L+L_start,i*ctx->L+L_end));
     126         904 :     PetscCall(BVGetMat(ctx->Y,&MC));
     127         904 :     if (B) {
     128         324 :       if (!i) {
     129          13 :         PetscCall(MatProductCreate(B,MV,NULL,&BMV));
     130          13 :         PetscCall(MatProductSetType(BMV,MATPRODUCT_AB));
     131          13 :         PetscCall(MatProductSetFromOptions(BMV));
     132          13 :         PetscCall(MatProductSymbolic(BMV));
     133             :       }
     134         324 :       PetscCall(MatProductNumeric(BMV));
     135         324 :       PetscCall(KSPMatSolve(ksp,BMV,MC));
     136         580 :     } else PetscCall(KSPMatSolve(ksp,MV,MC));
     137         904 :     PetscCall(BVRestoreMat(ctx->Y,&MC));
     138         904 :     if (ctx->usest && i<contour->npoints-1) PetscCall(KSPReset(ksp));
     139             :   }
     140          37 :   PetscCall(MatDestroy(&BMV));
     141          37 :   PetscCall(BVRestoreMat(V,&MV));
     142          37 :   PetscFunctionReturn(PETSC_SUCCESS);
     143             : }
     144             : 
     145          70 : static PetscErrorCode rescale_eig(EPS eps,PetscInt nv)
     146             : {
     147          70 :   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
     148          70 :   PetscInt       i;
     149          70 :   PetscScalar    center;
     150          70 :   PetscReal      radius,a,b,c,d,rgscale;
     151             : #if defined(PETSC_USE_COMPLEX)
     152          70 :   PetscReal      start_ang,end_ang,vscale,theta;
     153             : #endif
     154          70 :   PetscBool      isring,isellipse,isinterval;
     155             : 
     156          70 :   PetscFunctionBegin;
     157          70 :   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse));
     158          70 :   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring));
     159          70 :   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval));
     160          70 :   PetscCall(RGGetScale(eps->rg,&rgscale));
     161          70 :   if (isinterval) {
     162          20 :     PetscCall(RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d));
     163          20 :     if (c==d) {
     164         226 :       for (i=0;i<nv;i++) {
     165             : #if defined(PETSC_USE_COMPLEX)
     166         210 :         eps->eigr[i] = PetscRealPart(eps->eigr[i]);
     167             : #else
     168             :         eps->eigi[i] = 0;
     169             : #endif
     170             :       }
     171             :     }
     172             :   }
     173          70 :   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
     174          14 :     if (isellipse) {
     175           6 :       PetscCall(RGEllipseGetParameters(eps->rg,&center,&radius,NULL));
     176          54 :       for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
     177           8 :     } else if (isinterval) {
     178           4 :       PetscCall(RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d));
     179           4 :       if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
     180          14 :         for (i=0;i<nv;i++) {
     181          12 :           if (c==d) eps->eigr[i] = ((b-a)*(eps->eigr[i]+1.0)/2.0+a)*rgscale;
     182          12 :           if (a==b) {
     183             : #if defined(PETSC_USE_COMPLEX)
     184           0 :             eps->eigr[i] = ((d-c)*(eps->eigr[i]+1.0)/2.0+c)*rgscale*PETSC_i;
     185             : #else
     186             :             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
     187             : #endif
     188             :           }
     189             :         }
     190             :       } else {
     191           2 :         center = (b+a)/2.0+(d+c)/2.0*PETSC_PI;
     192           4 :         radius = PetscSqrtReal(PetscPowRealInt((b-a)/2.0,2)+PetscPowRealInt((d-c)/2.0,2));
     193          14 :         for (i=0;i<nv;i++) eps->eigr[i] = center + radius*eps->eigr[i];
     194             :       }
     195           4 :     } else if (isring) {  /* only supported in complex scalars */
     196             : #if defined(PETSC_USE_COMPLEX)
     197           4 :       PetscCall(RGRingGetParameters(eps->rg,&center,&radius,&vscale,&start_ang,&end_ang,NULL));
     198           4 :       if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
     199           0 :         for (i=0;i<nv;i++) {
     200           0 :           theta = (start_ang*2.0+(end_ang-start_ang)*(PetscRealPart(eps->eigr[i])+1.0))*PETSC_PI;
     201           0 :           eps->eigr[i] = rgscale*center + (rgscale*radius+PetscImaginaryPart(eps->eigr[i]))*PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
     202             :         }
     203             :       } else {
     204          34 :         for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
     205             :       }
     206             : #endif
     207             :     }
     208             :   }
     209          70 :   PetscFunctionReturn(PETSC_SUCCESS);
     210             : }
     211             : 
     212          33 : static PetscErrorCode EPSSetUp_CISS(EPS eps)
     213             : {
     214          33 :   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
     215          33 :   SlepcContourData contour;
     216          33 :   PetscBool        istrivial,isring,isellipse,isinterval,flg;
     217          33 :   PetscReal        c,d;
     218          33 :   PetscInt         nsplit;
     219          33 :   PetscRandom      rand;
     220          33 :   PetscObjectId    id;
     221          33 :   PetscObjectState state;
     222          33 :   Mat              A[2],Psplit[2];
     223          33 :   Vec              v0;
     224             : 
     225          33 :   PetscFunctionBegin;
     226          33 :   if (eps->ncv==PETSC_DEFAULT) {
     227          31 :     eps->ncv = ctx->L_max*ctx->M;
     228          31 :     if (eps->ncv>eps->n) {
     229          19 :       eps->ncv = eps->n;
     230          19 :       ctx->L_max = eps->ncv/ctx->M;
     231          19 :       PetscCheck(ctx->L_max,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot adjust solver parameters, try setting a smaller value of M (moment size)");
     232             :     }
     233             :   } else {
     234           2 :     PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
     235           2 :     ctx->L_max = eps->ncv/ctx->M;
     236           2 :     if (!ctx->L_max) {
     237           0 :       ctx->L_max = 1;
     238           0 :       eps->ncv = ctx->L_max*ctx->M;
     239             :     }
     240             :   }
     241          33 :   ctx->L = PetscMin(ctx->L,ctx->L_max);
     242          33 :   if (eps->max_it==PETSC_DEFAULT) eps->max_it = 5;
     243          33 :   if (eps->mpd==PETSC_DEFAULT) eps->mpd = eps->ncv;
     244          33 :   if (!eps->which) eps->which = EPS_ALL;
     245          33 :   PetscCheck(eps->which==EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only computing all eigenvalues");
     246          33 :   EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);
     247             : 
     248             :   /* check region */
     249          33 :   PetscCall(RGIsTrivial(eps->rg,&istrivial));
     250          33 :   PetscCheck(!istrivial,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
     251          33 :   PetscCall(RGGetComplement(eps->rg,&flg));
     252          33 :   PetscCheck(!flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
     253          33 :   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse));
     254          33 :   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring));
     255          33 :   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval));
     256          33 :   PetscCheck(isellipse || isring || isinterval,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Currently only implemented for interval, elliptic or ring regions");
     257             : 
     258             :   /* if the region has changed, then reset contour data */
     259          33 :   PetscCall(PetscObjectGetId((PetscObject)eps->rg,&id));
     260          33 :   PetscCall(PetscObjectStateGet((PetscObject)eps->rg,&state));
     261          33 :   if (ctx->rgid && (id != ctx->rgid || state != ctx->rgstate)) {
     262           0 :     PetscCall(SlepcContourDataDestroy(&ctx->contour));
     263           0 :     PetscCall(PetscInfo(eps,"Resetting the contour data structure due to a change of region\n"));
     264           0 :     ctx->rgid = id; ctx->rgstate = state;
     265             :   }
     266             : 
     267             : #if !defined(PETSC_USE_COMPLEX)
     268             :   PetscCheck(!isring,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Ring region only supported for complex scalars");
     269             : #endif
     270          33 :   if (isinterval) {
     271          10 :     PetscCall(RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d));
     272             : #if !defined(PETSC_USE_COMPLEX)
     273             :     PetscCheck(c==d && c==0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"In real scalars, endpoints of the imaginary axis must be both zero");
     274             : #endif
     275          10 :     if (!ctx->quad && c==d) ctx->quad = EPS_CISS_QUADRULE_CHEBYSHEV;
     276             :   }
     277          33 :   if (!ctx->quad) ctx->quad = EPS_CISS_QUADRULE_TRAPEZOIDAL;
     278             : 
     279             :   /* create contour data structure */
     280          33 :   if (!ctx->contour) {
     281           0 :     PetscCall(RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj));
     282           0 :     PetscCall(SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour));
     283             :   }
     284             : 
     285          33 :   PetscCall(EPSAllocateSolution(eps,0));
     286          33 :   PetscCall(BVGetRandomContext(eps->V,&rand));  /* make sure the random context is available when duplicating */
     287          33 :   if (ctx->weight) PetscCall(PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma));
     288          33 :   PetscCall(PetscMalloc4(ctx->N,&ctx->weight,ctx->N+1,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma));
     289             : 
     290             :   /* allocate basis vectors */
     291          33 :   PetscCall(BVDestroy(&ctx->S));
     292          33 :   PetscCall(BVDuplicateResize(eps->V,ctx->L*ctx->M,&ctx->S));
     293          33 :   PetscCall(BVDestroy(&ctx->V));
     294          33 :   PetscCall(BVDuplicateResize(eps->V,ctx->L,&ctx->V));
     295             : 
     296          33 :   PetscCall(STGetMatrix(eps->st,0,&A[0]));
     297          33 :   PetscCall(MatIsShell(A[0],&flg));
     298          33 :   PetscCheck(!flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Matrix type shell is not supported in this solver");
     299          33 :   if (eps->isgeneralized) PetscCall(STGetMatrix(eps->st,1,&A[1]));
     300             : 
     301          33 :   if (!ctx->usest_set) ctx->usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
     302          33 :   PetscCheck(!ctx->usest || ctx->npart==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The usest flag is not supported when partitions > 1");
     303             : 
     304             :   /* check if a user-defined split preconditioner has been set */
     305          33 :   PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
     306          33 :   if (nsplit) {
     307           4 :     PetscCall(STGetSplitPreconditionerTerm(eps->st,0,&Psplit[0]));
     308           4 :     if (eps->isgeneralized) PetscCall(STGetSplitPreconditionerTerm(eps->st,1,&Psplit[1]));
     309             :   }
     310             : 
     311          33 :   contour = ctx->contour;
     312          82 :   PetscCall(SlepcContourRedundantMat(contour,eps->isgeneralized?2:1,A,nsplit?Psplit:NULL));
     313          33 :   if (contour->pA) {
     314           8 :     PetscCall(BVGetColumn(ctx->V,0,&v0));
     315           8 :     PetscCall(SlepcContourScatterCreate(contour,v0));
     316           8 :     PetscCall(BVRestoreColumn(ctx->V,0,&v0));
     317           8 :     PetscCall(BVDestroy(&ctx->pV));
     318           8 :     PetscCall(BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->pV));
     319           8 :     PetscCall(BVSetSizesFromVec(ctx->pV,contour->xsub,eps->n));
     320           8 :     PetscCall(BVSetFromOptions(ctx->pV));
     321           8 :     PetscCall(BVResize(ctx->pV,ctx->L,PETSC_FALSE));
     322             :   }
     323             : 
     324          33 :   EPSCheckDefinite(eps);
     325          33 :   EPSCheckSinvertCondition(eps,ctx->usest," (with the usest flag set)");
     326             : 
     327          33 :   PetscCall(BVDestroy(&ctx->Y));
     328          33 :   if (contour->pA) {
     329           8 :     PetscCall(BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->Y));
     330           8 :     PetscCall(BVSetSizesFromVec(ctx->Y,contour->xsub,eps->n));
     331           8 :     PetscCall(BVSetFromOptions(ctx->Y));
     332           8 :     PetscCall(BVResize(ctx->Y,contour->npoints*ctx->L,PETSC_FALSE));
     333          25 :   } else PetscCall(BVDuplicateResize(eps->V,contour->npoints*ctx->L,&ctx->Y));
     334             : 
     335          33 :   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(DSSetType(eps->ds,DSGNHEP));
     336          27 :   else if (eps->isgeneralized) {
     337          11 :     if (eps->ishermitian && eps->ispositive) PetscCall(DSSetType(eps->ds,DSGHEP));
     338           2 :     else PetscCall(DSSetType(eps->ds,DSGNHEP));
     339             :   } else {
     340          16 :     if (eps->ishermitian) PetscCall(DSSetType(eps->ds,DSHEP));
     341           9 :     else PetscCall(DSSetType(eps->ds,DSNHEP));
     342             :   }
     343          33 :   PetscCall(DSAllocate(eps->ds,eps->ncv));
     344             : 
     345             : #if !defined(PETSC_USE_COMPLEX)
     346             :   PetscCall(EPSSetWorkVecs(eps,3));
     347             :   if (!eps->ishermitian) PetscCall(PetscInfo(eps,"Warning: complex eigenvalues are not calculated exactly without --with-scalar-type=complex in PETSc\n"));
     348             : #else
     349          33 :   PetscCall(EPSSetWorkVecs(eps,2));
     350             : #endif
     351          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     352             : }
     353             : 
     354          33 : static PetscErrorCode EPSSetUpSort_CISS(EPS eps)
     355             : {
     356          33 :   SlepcSC        sc;
     357             : 
     358          33 :   PetscFunctionBegin;
     359             :   /* fill sorting criterion context */
     360          33 :   eps->sc->comparison    = SlepcCompareSmallestReal;
     361          33 :   eps->sc->comparisonctx = NULL;
     362          33 :   eps->sc->map           = NULL;
     363          33 :   eps->sc->mapobj        = NULL;
     364             : 
     365             :   /* fill sorting criterion for DS */
     366          33 :   PetscCall(DSGetSlepcSC(eps->ds,&sc));
     367          33 :   sc->comparison    = SlepcCompareLargestMagnitude;
     368          33 :   sc->comparisonctx = NULL;
     369          33 :   sc->map           = NULL;
     370          33 :   sc->mapobj        = NULL;
     371          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     372             : }
     373             : 
     374          33 : static PetscErrorCode EPSSolve_CISS(EPS eps)
     375             : {
     376          33 :   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
     377          33 :   SlepcContourData contour = ctx->contour;
     378          33 :   Mat              A,B,X,M,pA,pB,T,J,Pa=NULL,Pb=NULL;
     379          33 :   BV               V;
     380          33 :   PetscInt         i,j,ld,nmat,L_add=0,nv=0,L_base=ctx->L,inner,nlocal,*inside,nsplit;
     381          33 :   PetscScalar      *Mu,*H0,*H1=NULL,*rr,*temp;
     382          33 :   PetscReal        error,max_error,norm;
     383          33 :   PetscBool        *fl1;
     384          33 :   Vec              si,si1=NULL,w[3];
     385          33 :   PetscRandom      rand;
     386             : #if defined(PETSC_USE_COMPLEX)
     387          33 :   PetscBool        isellipse;
     388          33 :   PetscReal        est_eig,eta;
     389             : #else
     390             :   PetscReal        normi;
     391             : #endif
     392             : 
     393          33 :   PetscFunctionBegin;
     394          33 :   w[0] = eps->work[0];
     395             : #if defined(PETSC_USE_COMPLEX)
     396          33 :   w[1] = NULL;
     397             : #else
     398             :   w[1] = eps->work[2];
     399             : #endif
     400          33 :   w[2] = eps->work[1];
     401          33 :   PetscCall(VecGetLocalSize(w[0],&nlocal));
     402          33 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     403          55 :   PetscCall(RGComputeQuadrature(eps->rg,ctx->quad==EPS_CISS_QUADRULE_CHEBYSHEV?RG_QUADRULE_CHEBYSHEV:RG_QUADRULE_TRAPEZOIDAL,ctx->N,ctx->omega,ctx->pp,ctx->weight));
     404          33 :   PetscCall(STGetNumMatrices(eps->st,&nmat));
     405          33 :   PetscCall(STGetMatrix(eps->st,0,&A));
     406          33 :   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
     407          20 :   else B = NULL;
     408          33 :   J = (contour->pA && nmat>1)? contour->pA[1]: B;
     409          33 :   V = contour->pA? ctx->pV: ctx->V;
     410          33 :   if (!ctx->usest) {
     411          18 :     T = contour->pA? contour->pA[0]: A;
     412          18 :     PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
     413          18 :     if (nsplit) {
     414           3 :       if (contour->pA) {
     415           2 :         Pa = contour->pP[0];
     416           2 :         if (nsplit>1) Pb = contour->pP[1];
     417             :       } else {
     418           1 :         PetscCall(STGetSplitPreconditionerTerm(eps->st,0,&Pa));
     419           1 :         if (nsplit>1) PetscCall(STGetSplitPreconditionerTerm(eps->st,1,&Pb));
     420             :       }
     421             :     }
     422          18 :     PetscCall(EPSCISSSetUp(eps,T,J,Pa,Pb));
     423             :   }
     424          33 :   PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
     425          33 :   PetscCall(BVSetRandomSign(ctx->V));
     426          33 :   PetscCall(BVGetRandomContext(ctx->V,&rand));
     427             : 
     428          33 :   if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
     429          33 :   PetscCall(EPSCISSSolve(eps,J,V,0,ctx->L));
     430             : #if defined(PETSC_USE_COMPLEX)
     431          33 :   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse));
     432          33 :   if (isellipse) {
     433          20 :     PetscCall(BVTraceQuadrature(ctx->Y,ctx->V,ctx->L,ctx->L,ctx->weight,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj,&est_eig));
     434          20 :     PetscCall(PetscInfo(eps,"Estimated eigenvalue count: %f\n",(double)est_eig));
     435          20 :     eta = PetscPowReal(10.0,-PetscLog10Real(eps->tol)/ctx->N);
     436          20 :     L_add = PetscMax(0,(PetscInt)PetscCeilReal((est_eig*eta)/ctx->M)-ctx->L);
     437          20 :     if (L_add>ctx->L_max-ctx->L) {
     438           1 :       PetscCall(PetscInfo(eps,"Number of eigenvalues inside the contour path may be too large\n"));
     439           1 :       L_add = ctx->L_max-ctx->L;
     440             :     }
     441             :   }
     442             : #endif
     443          20 :   if (L_add>0) {
     444           2 :     PetscCall(PetscInfo(eps,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by Estimate #Eig\n",ctx->L,ctx->L+L_add));
     445           2 :     PetscCall(BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints));
     446           2 :     PetscCall(BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add));
     447           2 :     PetscCall(BVSetRandomSign(ctx->V));
     448           2 :     if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
     449           2 :     ctx->L += L_add;
     450           2 :     PetscCall(EPSCISSSolve(eps,J,V,ctx->L-L_add,ctx->L));
     451             :   }
     452          33 :   PetscCall(PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0));
     453          33 :   for (i=0;i<ctx->refine_blocksize;i++) {
     454           1 :     PetscCall(BVDotQuadrature(ctx->Y,(contour->pA)?ctx->pV:ctx->V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj));
     455           1 :     PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
     456           1 :     PetscCall(PetscLogEventBegin(EPS_CISS_SVD,eps,0,0,0));
     457           1 :     PetscCall(SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv));
     458           1 :     PetscCall(PetscLogEventEnd(EPS_CISS_SVD,eps,0,0,0));
     459           1 :     if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
     460           0 :     L_add = L_base;
     461           0 :     if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
     462           0 :     PetscCall(PetscInfo(eps,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by SVD(H0)\n",ctx->L,ctx->L+L_add));
     463           0 :     PetscCall(BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints));
     464           0 :     PetscCall(BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add));
     465           0 :     PetscCall(BVSetRandomSign(ctx->V));
     466           0 :     if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
     467           0 :     ctx->L += L_add;
     468           0 :     PetscCall(EPSCISSSolve(eps,J,V,ctx->L-L_add,ctx->L));
     469           0 :     if (L_add) {
     470           0 :       PetscCall(PetscFree2(Mu,H0));
     471           0 :       PetscCall(PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0));
     472             :     }
     473             :   }
     474          33 :   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1));
     475             : 
     476          68 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     477          35 :     eps->its++;
     478          35 :     for (inner=0;inner<=ctx->refine_inner;inner++) {
     479          35 :       if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
     480           7 :         PetscCall(BVDotQuadrature(ctx->Y,(contour->pA)?ctx->pV:ctx->V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj));
     481           7 :         PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
     482           7 :         PetscCall(PetscLogEventBegin(EPS_CISS_SVD,eps,0,0,0));
     483           7 :         PetscCall(SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv));
     484           7 :         PetscCall(PetscLogEventEnd(EPS_CISS_SVD,eps,0,0,0));
     485             :         break;
     486             :       } else {
     487          28 :         PetscCall(BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj));
     488          28 :         PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
     489          28 :         PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
     490          28 :         PetscCall(BVCopy(ctx->S,ctx->V));
     491          28 :         PetscCall(BVSVDAndRank(ctx->S,ctx->M,ctx->L,ctx->delta,BV_SVD_METHOD_REFINE,H0,ctx->sigma,&nv));
     492          28 :         if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
     493           0 :           if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
     494           0 :           PetscCall(EPSCISSSolve(eps,J,V,0,ctx->L));
     495             :         } else break;
     496             :       }
     497             :     }
     498          35 :     eps->nconv = 0;
     499          35 :     if (nv == 0) eps->reason = EPS_CONVERGED_TOL;
     500             :     else {
     501          35 :       PetscCall(DSSetDimensions(eps->ds,nv,0,0));
     502          35 :       PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
     503             : 
     504          35 :       if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
     505           7 :         PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
     506           7 :         PetscCall(CISS_BlockHankel(Mu,1,ctx->L,ctx->M,H1));
     507           7 :         PetscCall(DSGetArray(eps->ds,DS_MAT_A,&temp));
     508          58 :         for (j=0;j<nv;j++) {
     509         438 :           for (i=0;i<nv;i++) {
     510         387 :             temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
     511             :           }
     512             :         }
     513           7 :         PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&temp));
     514           7 :         PetscCall(DSGetArray(eps->ds,DS_MAT_B,&temp));
     515          58 :         for (j=0;j<nv;j++) {
     516         438 :           for (i=0;i<nv;i++) {
     517         387 :             temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
     518             :           }
     519             :         }
     520           7 :         PetscCall(DSRestoreArray(eps->ds,DS_MAT_B,&temp));
     521             :       } else {
     522          28 :         PetscCall(BVSetActiveColumns(ctx->S,0,nv));
     523          28 :         PetscCall(DSGetMat(eps->ds,DS_MAT_A,&pA));
     524          28 :         PetscCall(MatZeroEntries(pA));
     525          28 :         PetscCall(BVMatProject(ctx->S,A,ctx->S,pA));
     526          28 :         PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&pA));
     527          28 :         if (B) {
     528          11 :           PetscCall(DSGetMat(eps->ds,DS_MAT_B,&pB));
     529          11 :           PetscCall(MatZeroEntries(pB));
     530          11 :           PetscCall(BVMatProject(ctx->S,B,ctx->S,pB));
     531          11 :           PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&pB));
     532             :         }
     533             :       }
     534             : 
     535          35 :       PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     536          35 :       PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     537             : 
     538          35 :       PetscCall(PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr));
     539          35 :       PetscCall(rescale_eig(eps,nv));
     540          35 :       PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
     541          35 :       PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
     542          35 :       PetscCall(SlepcCISS_isGhost(X,nv,ctx->sigma,ctx->spurious_threshold,fl1));
     543          35 :       PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
     544          35 :       PetscCall(RGCheckInside(eps->rg,nv,eps->eigr,eps->eigi,inside));
     545        1153 :       for (i=0;i<nv;i++) {
     546        1118 :         if (fl1[i] && inside[i]>=0) {
     547         269 :           rr[i] = 1.0;
     548         269 :           eps->nconv++;
     549         849 :         } else rr[i] = 0.0;
     550             :       }
     551          35 :       PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,rr,NULL,&eps->nconv));
     552          35 :       PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     553          35 :       PetscCall(rescale_eig(eps,nv));
     554          35 :       PetscCall(PetscFree3(fl1,inside,rr));
     555          35 :       PetscCall(BVSetActiveColumns(eps->V,0,nv));
     556          35 :       if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
     557           7 :         PetscCall(BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj));
     558           7 :         PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
     559           7 :         PetscCall(BVCopy(ctx->S,ctx->V));
     560           7 :         PetscCall(BVSetActiveColumns(ctx->S,0,nv));
     561             :       }
     562          35 :       PetscCall(BVCopy(ctx->S,eps->V));
     563             : 
     564          35 :       PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
     565          35 :       PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
     566          35 :       PetscCall(BVMultInPlace(ctx->S,X,0,eps->nconv));
     567          35 :       if (eps->ishermitian) PetscCall(BVMultInPlace(eps->V,X,0,eps->nconv));
     568          35 :       PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
     569             :       max_error = 0.0;
     570         304 :       for (i=0;i<eps->nconv;i++) {
     571         269 :         PetscCall(BVGetColumn(ctx->S,i,&si));
     572             : #if !defined(PETSC_USE_COMPLEX)
     573             :         if (eps->eigi[i]!=0.0) PetscCall(BVGetColumn(ctx->S,i+1,&si1));
     574             : #endif
     575         269 :         PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,eps->eigr[i],eps->eigi[i],si,si1,w,&error));
     576         269 :         if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {  /* vector is not normalized */
     577          35 :           PetscCall(VecNorm(si,NORM_2,&norm));
     578             : #if !defined(PETSC_USE_COMPLEX)
     579             :           if (eps->eigi[i]!=0.0) {
     580             :             PetscCall(VecNorm(si1,NORM_2,&normi));
     581             :             norm = SlepcAbsEigenvalue(norm,normi);
     582             :           }
     583             : #endif
     584          35 :           error /= norm;
     585             :         }
     586         269 :         PetscCall((*eps->converged)(eps,eps->eigr[i],eps->eigi[i],error,&error,eps->convergedctx));
     587         269 :         PetscCall(BVRestoreColumn(ctx->S,i,&si));
     588             : #if !defined(PETSC_USE_COMPLEX)
     589             :         if (eps->eigi[i]!=0.0) {
     590             :           PetscCall(BVRestoreColumn(ctx->S,i+1,&si1));
     591             :           i++;
     592             :         }
     593             : #endif
     594         356 :         max_error = PetscMax(max_error,error);
     595             :       }
     596             : 
     597          35 :       if (max_error <= eps->tol) eps->reason = EPS_CONVERGED_TOL;
     598           2 :       else if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
     599             :       else {
     600           2 :         if (eps->nconv > ctx->L) nv = eps->nconv;
     601           0 :         else if (ctx->L > nv) nv = ctx->L;
     602           2 :         nv = PetscMin(nv,ctx->L*ctx->M);
     603           2 :         PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M));
     604           2 :         PetscCall(MatSetRandom(M,rand));
     605           2 :         PetscCall(BVSetActiveColumns(ctx->S,0,nv));
     606           2 :         PetscCall(BVMultInPlace(ctx->S,M,0,ctx->L));
     607           2 :         PetscCall(MatDestroy(&M));
     608           2 :         PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
     609           2 :         PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
     610           2 :         PetscCall(BVCopy(ctx->S,ctx->V));
     611           2 :         if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
     612          70 :         PetscCall(EPSCISSSolve(eps,J,V,0,ctx->L));
     613             :       }
     614             :     }
     615             :   }
     616          33 :   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(PetscFree(H1));
     617          33 :   PetscCall(PetscFree2(Mu,H0));
     618          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     619             : }
     620             : 
     621          33 : static PetscErrorCode EPSComputeVectors_CISS(EPS eps)
     622             : {
     623          33 :   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
     624          33 :   PetscInt       n;
     625          33 :   Mat            Z,B=NULL;
     626             : 
     627          33 :   PetscFunctionBegin;
     628          33 :   if (eps->ishermitian) {
     629          18 :     if (eps->isgeneralized && !eps->ispositive) PetscCall(EPSComputeVectors_Indefinite(eps));
     630          18 :     else PetscCall(EPSComputeVectors_Hermitian(eps));
     631          18 :     if (eps->isgeneralized && eps->ispositive && ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
     632             :       /* normalize to have unit B-norm */
     633           2 :       PetscCall(STGetMatrix(eps->st,1,&B));
     634           2 :       PetscCall(BVSetMatrix(eps->V,B,PETSC_FALSE));
     635           2 :       PetscCall(BVNormalize(eps->V,NULL));
     636           2 :       PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
     637             :     }
     638          18 :     PetscFunctionReturn(PETSC_SUCCESS);
     639             :   }
     640          15 :   PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
     641          15 :   PetscCall(BVSetActiveColumns(eps->V,0,n));
     642             : 
     643             :   /* right eigenvectors */
     644          15 :   PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
     645             : 
     646             :   /* V = V * Z */
     647          15 :   PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
     648          15 :   PetscCall(BVMultInPlace(eps->V,Z,0,n));
     649          15 :   PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
     650          15 :   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
     651             : 
     652             :   /* normalize */
     653          15 :   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(BVNormalize(eps->V,NULL));
     654          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     655             : }
     656             : 
     657          15 : static PetscErrorCode EPSCISSSetSizes_CISS(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
     658             : {
     659          15 :   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
     660          15 :   PetscInt       oN,oL,oM,oLmax,onpart;
     661          15 :   PetscMPIInt    size;
     662             : 
     663          15 :   PetscFunctionBegin;
     664          15 :   oN = ctx->N;
     665          15 :   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
     666           0 :     if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
     667             :   } else {
     668          15 :     PetscCheck(ip>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
     669          15 :     PetscCheck(ip%2==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
     670          15 :     if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
     671             :   }
     672          15 :   oL = ctx->L;
     673          15 :   if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
     674           0 :     ctx->L = 16;
     675             :   } else {
     676          15 :     PetscCheck(bs>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
     677          15 :     ctx->L = bs;
     678             :   }
     679          15 :   oM = ctx->M;
     680          15 :   if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
     681           0 :     ctx->M = ctx->N/4;
     682             :   } else {
     683          15 :     PetscCheck(ms>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
     684          15 :     PetscCheck(ms<=ctx->N,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be less than or equal to the number of integration points");
     685          15 :     ctx->M = ms;
     686             :   }
     687          15 :   onpart = ctx->npart;
     688          15 :   if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
     689           0 :     ctx->npart = 1;
     690             :   } else {
     691          15 :     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
     692          15 :     PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
     693          15 :     ctx->npart = npart;
     694             :   }
     695          15 :   oLmax = ctx->L_max;
     696          15 :   if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
     697           0 :     ctx->L_max = 64;
     698             :   } else {
     699          15 :     PetscCheck(bsmax>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
     700          15 :     ctx->L_max = PetscMax(bsmax,ctx->L);
     701             :   }
     702          15 :   if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) {
     703          14 :     PetscCall(SlepcContourDataDestroy(&ctx->contour));
     704          14 :     PetscCall(PetscInfo(eps,"Resetting the contour data structure due to a change of parameters\n"));
     705          14 :     eps->state = EPS_STATE_INITIAL;
     706             :   }
     707          15 :   ctx->isreal = realmats;
     708          15 :   if (oL != ctx->L || oM != ctx->M || oLmax != ctx->L_max) eps->state = EPS_STATE_INITIAL;
     709          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     710             : }
     711             : 
     712             : /*@
     713             :    EPSCISSSetSizes - Sets the values of various size parameters in the CISS solver.
     714             : 
     715             :    Logically Collective
     716             : 
     717             :    Input Parameters:
     718             : +  eps   - the eigenproblem solver context
     719             : .  ip    - number of integration points
     720             : .  bs    - block size
     721             : .  ms    - moment size
     722             : .  npart - number of partitions when splitting the communicator
     723             : .  bsmax - max block size
     724             : -  realmats - A and B are real
     725             : 
     726             :    Options Database Keys:
     727             : +  -eps_ciss_integration_points - Sets the number of integration points
     728             : .  -eps_ciss_blocksize - Sets the block size
     729             : .  -eps_ciss_moments - Sets the moment size
     730             : .  -eps_ciss_partitions - Sets the number of partitions
     731             : .  -eps_ciss_maxblocksize - Sets the maximum block size
     732             : -  -eps_ciss_realmats - A and B are real
     733             : 
     734             :    Note:
     735             :    The default number of partitions is 1. This means the internal KSP object is shared
     736             :    among all processes of the EPS communicator. Otherwise, the communicator is split
     737             :    into npart communicators, so that npart KSP solves proceed simultaneously.
     738             : 
     739             :    Level: advanced
     740             : 
     741             : .seealso: EPSCISSGetSizes()
     742             : @*/
     743          15 : PetscErrorCode EPSCISSSetSizes(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
     744             : {
     745          15 :   PetscFunctionBegin;
     746          15 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     747          60 :   PetscValidLogicalCollectiveInt(eps,ip,2);
     748          60 :   PetscValidLogicalCollectiveInt(eps,bs,3);
     749          60 :   PetscValidLogicalCollectiveInt(eps,ms,4);
     750          60 :   PetscValidLogicalCollectiveInt(eps,npart,5);
     751          60 :   PetscValidLogicalCollectiveInt(eps,bsmax,6);
     752          60 :   PetscValidLogicalCollectiveBool(eps,realmats,7);
     753          15 :   PetscTryMethod(eps,"EPSCISSSetSizes_C",(EPS,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(eps,ip,bs,ms,npart,bsmax,realmats));
     754          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     755             : }
     756             : 
     757          33 : static PetscErrorCode EPSCISSGetSizes_CISS(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
     758             : {
     759          33 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
     760             : 
     761          33 :   PetscFunctionBegin;
     762          33 :   if (ip) *ip = ctx->N;
     763          33 :   if (bs) *bs = ctx->L;
     764          33 :   if (ms) *ms = ctx->M;
     765          33 :   if (npart) *npart = ctx->npart;
     766          33 :   if (bsmax) *bsmax = ctx->L_max;
     767          33 :   if (realmats) *realmats = ctx->isreal;
     768          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     769             : }
     770             : 
     771             : /*@
     772             :    EPSCISSGetSizes - Gets the values of various size parameters in the CISS solver.
     773             : 
     774             :    Not Collective
     775             : 
     776             :    Input Parameter:
     777             : .  eps - the eigenproblem solver context
     778             : 
     779             :    Output Parameters:
     780             : +  ip    - number of integration points
     781             : .  bs    - block size
     782             : .  ms    - moment size
     783             : .  npart - number of partitions when splitting the communicator
     784             : .  bsmax - max block size
     785             : -  realmats - A and B are real
     786             : 
     787             :    Level: advanced
     788             : 
     789             : .seealso: EPSCISSSetSizes()
     790             : @*/
     791          33 : PetscErrorCode EPSCISSGetSizes(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
     792             : {
     793          33 :   PetscFunctionBegin;
     794          33 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     795          33 :   PetscUseMethod(eps,"EPSCISSGetSizes_C",(EPS,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(eps,ip,bs,ms,npart,bsmax,realmats));
     796          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     797             : }
     798             : 
     799           2 : static PetscErrorCode EPSCISSSetThreshold_CISS(EPS eps,PetscReal delta,PetscReal spur)
     800             : {
     801           2 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
     802             : 
     803           2 :   PetscFunctionBegin;
     804           2 :   if (delta == (PetscReal)PETSC_DEFAULT) {
     805           0 :     ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
     806             :   } else {
     807           2 :     PetscCheck(delta>0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
     808           2 :     ctx->delta = delta;
     809             :   }
     810           2 :   if (spur == (PetscReal)PETSC_DEFAULT) {
     811           0 :     ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
     812             :   } else {
     813           2 :     PetscCheck(spur>0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
     814           2 :     ctx->spurious_threshold = spur;
     815             :   }
     816           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     817             : }
     818             : 
     819             : /*@
     820             :    EPSCISSSetThreshold - Sets the values of various threshold parameters in
     821             :    the CISS solver.
     822             : 
     823             :    Logically Collective
     824             : 
     825             :    Input Parameters:
     826             : +  eps   - the eigenproblem solver context
     827             : .  delta - threshold for numerical rank
     828             : -  spur  - spurious threshold (to discard spurious eigenpairs)
     829             : 
     830             :    Options Database Keys:
     831             : +  -eps_ciss_delta - Sets the delta
     832             : -  -eps_ciss_spurious_threshold - Sets the spurious threshold
     833             : 
     834             :    Level: advanced
     835             : 
     836             : .seealso: EPSCISSGetThreshold()
     837             : @*/
     838           2 : PetscErrorCode EPSCISSSetThreshold(EPS eps,PetscReal delta,PetscReal spur)
     839             : {
     840           2 :   PetscFunctionBegin;
     841           2 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     842           8 :   PetscValidLogicalCollectiveReal(eps,delta,2);
     843           8 :   PetscValidLogicalCollectiveReal(eps,spur,3);
     844           2 :   PetscTryMethod(eps,"EPSCISSSetThreshold_C",(EPS,PetscReal,PetscReal),(eps,delta,spur));
     845           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     846             : }
     847             : 
     848          33 : static PetscErrorCode EPSCISSGetThreshold_CISS(EPS eps,PetscReal *delta,PetscReal *spur)
     849             : {
     850          33 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
     851             : 
     852          33 :   PetscFunctionBegin;
     853          33 :   if (delta) *delta = ctx->delta;
     854          33 :   if (spur)  *spur = ctx->spurious_threshold;
     855          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     856             : }
     857             : 
     858             : /*@
     859             :    EPSCISSGetThreshold - Gets the values of various threshold parameters
     860             :    in the CISS solver.
     861             : 
     862             :    Not Collective
     863             : 
     864             :    Input Parameter:
     865             : .  eps - the eigenproblem solver context
     866             : 
     867             :    Output Parameters:
     868             : +  delta - threshold for numerical rank
     869             : -  spur  - spurious threshold (to discard spurious eigenpairs)
     870             : 
     871             :    Level: advanced
     872             : 
     873             : .seealso: EPSCISSSetThreshold()
     874             : @*/
     875          33 : PetscErrorCode EPSCISSGetThreshold(EPS eps,PetscReal *delta,PetscReal *spur)
     876             : {
     877          33 :   PetscFunctionBegin;
     878          33 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     879          33 :   PetscUseMethod(eps,"EPSCISSGetThreshold_C",(EPS,PetscReal*,PetscReal*),(eps,delta,spur));
     880          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     881             : }
     882             : 
     883           1 : static PetscErrorCode EPSCISSSetRefinement_CISS(EPS eps,PetscInt inner,PetscInt blsize)
     884             : {
     885           1 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
     886             : 
     887           1 :   PetscFunctionBegin;
     888           1 :   if (inner == PETSC_DEFAULT) {
     889           0 :     ctx->refine_inner = 0;
     890             :   } else {
     891           1 :     PetscCheck(inner>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
     892           1 :     ctx->refine_inner = inner;
     893             :   }
     894           1 :   if (blsize == PETSC_DEFAULT) {
     895           0 :     ctx->refine_blocksize = 0;
     896             :   } else {
     897           1 :     PetscCheck(blsize>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
     898           1 :     ctx->refine_blocksize = blsize;
     899             :   }
     900           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     901             : }
     902             : 
     903             : /*@
     904             :    EPSCISSSetRefinement - Sets the values of various refinement parameters
     905             :    in the CISS solver.
     906             : 
     907             :    Logically Collective
     908             : 
     909             :    Input Parameters:
     910             : +  eps    - the eigenproblem solver context
     911             : .  inner  - number of iterative refinement iterations (inner loop)
     912             : -  blsize - number of iterative refinement iterations (blocksize loop)
     913             : 
     914             :    Options Database Keys:
     915             : +  -eps_ciss_refine_inner - Sets number of inner iterations
     916             : -  -eps_ciss_refine_blocksize - Sets number of blocksize iterations
     917             : 
     918             :    Level: advanced
     919             : 
     920             : .seealso: EPSCISSGetRefinement()
     921             : @*/
     922           1 : PetscErrorCode EPSCISSSetRefinement(EPS eps,PetscInt inner,PetscInt blsize)
     923             : {
     924           1 :   PetscFunctionBegin;
     925           1 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     926           4 :   PetscValidLogicalCollectiveInt(eps,inner,2);
     927           4 :   PetscValidLogicalCollectiveInt(eps,blsize,3);
     928           1 :   PetscTryMethod(eps,"EPSCISSSetRefinement_C",(EPS,PetscInt,PetscInt),(eps,inner,blsize));
     929           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     930             : }
     931             : 
     932          33 : static PetscErrorCode EPSCISSGetRefinement_CISS(EPS eps,PetscInt *inner,PetscInt *blsize)
     933             : {
     934          33 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
     935             : 
     936          33 :   PetscFunctionBegin;
     937          33 :   if (inner)  *inner = ctx->refine_inner;
     938          33 :   if (blsize) *blsize = ctx->refine_blocksize;
     939          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     940             : }
     941             : 
     942             : /*@
     943             :    EPSCISSGetRefinement - Gets the values of various refinement parameters
     944             :    in the CISS solver.
     945             : 
     946             :    Not Collective
     947             : 
     948             :    Input Parameter:
     949             : .  eps - the eigenproblem solver context
     950             : 
     951             :    Output Parameters:
     952             : +  inner  - number of iterative refinement iterations (inner loop)
     953             : -  blsize - number of iterative refinement iterations (blocksize loop)
     954             : 
     955             :    Level: advanced
     956             : 
     957             : .seealso: EPSCISSSetRefinement()
     958             : @*/
     959          33 : PetscErrorCode EPSCISSGetRefinement(EPS eps, PetscInt *inner, PetscInt *blsize)
     960             : {
     961          33 :   PetscFunctionBegin;
     962          33 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     963          33 :   PetscUseMethod(eps,"EPSCISSGetRefinement_C",(EPS,PetscInt*,PetscInt*),(eps,inner,blsize));
     964          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     965             : }
     966             : 
     967          11 : static PetscErrorCode EPSCISSSetUseST_CISS(EPS eps,PetscBool usest)
     968             : {
     969          11 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
     970             : 
     971          11 :   PetscFunctionBegin;
     972          11 :   ctx->usest     = usest;
     973          11 :   ctx->usest_set = PETSC_TRUE;
     974          11 :   eps->state     = EPS_STATE_INITIAL;
     975          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     976             : }
     977             : 
     978             : /*@
     979             :    EPSCISSSetUseST - Sets a flag indicating that the CISS solver will
     980             :    use the ST object for the linear solves.
     981             : 
     982             :    Logically Collective
     983             : 
     984             :    Input Parameters:
     985             : +  eps    - the eigenproblem solver context
     986             : -  usest  - boolean flag to use the ST object or not
     987             : 
     988             :    Options Database Keys:
     989             : .  -eps_ciss_usest <bool> - whether the ST object will be used or not
     990             : 
     991             :    Level: advanced
     992             : 
     993             : .seealso: EPSCISSGetUseST()
     994             : @*/
     995          11 : PetscErrorCode EPSCISSSetUseST(EPS eps,PetscBool usest)
     996             : {
     997          11 :   PetscFunctionBegin;
     998          11 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     999          44 :   PetscValidLogicalCollectiveBool(eps,usest,2);
    1000          11 :   PetscTryMethod(eps,"EPSCISSSetUseST_C",(EPS,PetscBool),(eps,usest));
    1001          11 :   PetscFunctionReturn(PETSC_SUCCESS);
    1002             : }
    1003             : 
    1004          33 : static PetscErrorCode EPSCISSGetUseST_CISS(EPS eps,PetscBool *usest)
    1005             : {
    1006          33 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
    1007             : 
    1008          33 :   PetscFunctionBegin;
    1009          33 :   *usest = ctx->usest;
    1010          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1011             : }
    1012             : 
    1013             : /*@
    1014             :    EPSCISSGetUseST - Gets the flag for using the ST object
    1015             :    in the CISS solver.
    1016             : 
    1017             :    Not Collective
    1018             : 
    1019             :    Input Parameter:
    1020             : .  eps - the eigenproblem solver context
    1021             : 
    1022             :    Output Parameters:
    1023             : .  usest - boolean flag indicating if the ST object is being used
    1024             : 
    1025             :    Level: advanced
    1026             : 
    1027             : .seealso: EPSCISSSetUseST()
    1028             : @*/
    1029          33 : PetscErrorCode EPSCISSGetUseST(EPS eps,PetscBool *usest)
    1030             : {
    1031          33 :   PetscFunctionBegin;
    1032          33 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
    1033          33 :   PetscAssertPointer(usest,2);
    1034          33 :   PetscUseMethod(eps,"EPSCISSGetUseST_C",(EPS,PetscBool*),(eps,usest));
    1035          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1036             : }
    1037             : 
    1038           6 : static PetscErrorCode EPSCISSSetQuadRule_CISS(EPS eps,EPSCISSQuadRule quad)
    1039             : {
    1040           6 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
    1041             : 
    1042           6 :   PetscFunctionBegin;
    1043           6 :   if (ctx->quad != quad) {
    1044           6 :     ctx->quad  = quad;
    1045           6 :     eps->state = EPS_STATE_INITIAL;
    1046             :   }
    1047           6 :   PetscFunctionReturn(PETSC_SUCCESS);
    1048             : }
    1049             : 
    1050             : /*@
    1051             :    EPSCISSSetQuadRule - Sets the quadrature rule used in the CISS solver.
    1052             : 
    1053             :    Logically Collective
    1054             : 
    1055             :    Input Parameters:
    1056             : +  eps  - the eigenproblem solver context
    1057             : -  quad - the quadrature rule
    1058             : 
    1059             :    Options Database Key:
    1060             : .  -eps_ciss_quadrule - Sets the quadrature rule (either 'trapezoidal' or
    1061             :                            'chebyshev')
    1062             : 
    1063             :    Notes:
    1064             :    By default, the trapezoidal rule is used (EPS_CISS_QUADRULE_TRAPEZOIDAL).
    1065             : 
    1066             :    If the 'chebyshev' option is specified (EPS_CISS_QUADRULE_CHEBYSHEV), then
    1067             :    Chebyshev points are used as quadrature points.
    1068             : 
    1069             :    Level: advanced
    1070             : 
    1071             : .seealso: EPSCISSGetQuadRule(), EPSCISSQuadRule
    1072             : @*/
    1073           6 : PetscErrorCode EPSCISSSetQuadRule(EPS eps,EPSCISSQuadRule quad)
    1074             : {
    1075           6 :   PetscFunctionBegin;
    1076           6 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
    1077          24 :   PetscValidLogicalCollectiveEnum(eps,quad,2);
    1078           6 :   PetscTryMethod(eps,"EPSCISSSetQuadRule_C",(EPS,EPSCISSQuadRule),(eps,quad));
    1079           6 :   PetscFunctionReturn(PETSC_SUCCESS);
    1080             : }
    1081             : 
    1082           1 : static PetscErrorCode EPSCISSGetQuadRule_CISS(EPS eps,EPSCISSQuadRule *quad)
    1083             : {
    1084           1 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
    1085             : 
    1086           1 :   PetscFunctionBegin;
    1087           1 :   *quad = ctx->quad;
    1088           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1089             : }
    1090             : 
    1091             : /*@
    1092             :    EPSCISSGetQuadRule - Gets the quadrature rule used in the CISS solver.
    1093             : 
    1094             :    Not Collective
    1095             : 
    1096             :    Input Parameter:
    1097             : .  eps - the eigenproblem solver context
    1098             : 
    1099             :    Output Parameters:
    1100             : .  quad - quadrature rule
    1101             : 
    1102             :    Level: advanced
    1103             : 
    1104             : .seealso: EPSCISSSetQuadRule() EPSCISSQuadRule
    1105             : @*/
    1106           1 : PetscErrorCode EPSCISSGetQuadRule(EPS eps,EPSCISSQuadRule *quad)
    1107             : {
    1108           1 :   PetscFunctionBegin;
    1109           1 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
    1110           1 :   PetscAssertPointer(quad,2);
    1111           1 :   PetscUseMethod(eps,"EPSCISSGetQuadRule_C",(EPS,EPSCISSQuadRule*),(eps,quad));
    1112           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1113             : }
    1114             : 
    1115           8 : static PetscErrorCode EPSCISSSetExtraction_CISS(EPS eps,EPSCISSExtraction extraction)
    1116             : {
    1117           8 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
    1118             : 
    1119           8 :   PetscFunctionBegin;
    1120           8 :   if (ctx->extraction != extraction) {
    1121           6 :     ctx->extraction = extraction;
    1122           6 :     eps->state      = EPS_STATE_INITIAL;
    1123             :   }
    1124           8 :   PetscFunctionReturn(PETSC_SUCCESS);
    1125             : }
    1126             : 
    1127             : /*@
    1128             :    EPSCISSSetExtraction - Sets the extraction technique used in the CISS solver.
    1129             : 
    1130             :    Logically Collective
    1131             : 
    1132             :    Input Parameters:
    1133             : +  eps        - the eigenproblem solver context
    1134             : -  extraction - the extraction technique
    1135             : 
    1136             :    Options Database Key:
    1137             : .  -eps_ciss_extraction - Sets the extraction technique (either 'ritz' or
    1138             :                            'hankel')
    1139             : 
    1140             :    Notes:
    1141             :    By default, the Rayleigh-Ritz extraction is used (EPS_CISS_EXTRACTION_RITZ).
    1142             : 
    1143             :    If the 'hankel' option is specified (EPS_CISS_EXTRACTION_HANKEL), then
    1144             :    the Block Hankel method is used for extracting eigenpairs.
    1145             : 
    1146             :    Level: advanced
    1147             : 
    1148             : .seealso: EPSCISSGetExtraction(), EPSCISSExtraction
    1149             : @*/
    1150           8 : PetscErrorCode EPSCISSSetExtraction(EPS eps,EPSCISSExtraction extraction)
    1151             : {
    1152           8 :   PetscFunctionBegin;
    1153           8 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
    1154          32 :   PetscValidLogicalCollectiveEnum(eps,extraction,2);
    1155           8 :   PetscTryMethod(eps,"EPSCISSSetExtraction_C",(EPS,EPSCISSExtraction),(eps,extraction));
    1156           8 :   PetscFunctionReturn(PETSC_SUCCESS);
    1157             : }
    1158             : 
    1159           1 : static PetscErrorCode EPSCISSGetExtraction_CISS(EPS eps,EPSCISSExtraction *extraction)
    1160             : {
    1161           1 :   EPS_CISS *ctx = (EPS_CISS*)eps->data;
    1162             : 
    1163           1 :   PetscFunctionBegin;
    1164           1 :   *extraction = ctx->extraction;
    1165           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1166             : }
    1167             : 
    1168             : /*@
    1169             :    EPSCISSGetExtraction - Gets the extraction technique used in the CISS solver.
    1170             : 
    1171             :    Not Collective
    1172             : 
    1173             :    Input Parameter:
    1174             : .  eps - the eigenproblem solver context
    1175             : 
    1176             :    Output Parameters:
    1177             : .  extraction - extraction technique
    1178             : 
    1179             :    Level: advanced
    1180             : 
    1181             : .seealso: EPSCISSSetExtraction() EPSCISSExtraction
    1182             : @*/
    1183           1 : PetscErrorCode EPSCISSGetExtraction(EPS eps,EPSCISSExtraction *extraction)
    1184             : {
    1185           1 :   PetscFunctionBegin;
    1186           1 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
    1187           1 :   PetscAssertPointer(extraction,2);
    1188           1 :   PetscUseMethod(eps,"EPSCISSGetExtraction_C",(EPS,EPSCISSExtraction*),(eps,extraction));
    1189           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1190             : }
    1191             : 
    1192          33 : static PetscErrorCode EPSCISSGetKSPs_CISS(EPS eps,PetscInt *nsolve,KSP **ksp)
    1193             : {
    1194          33 :   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
    1195          33 :   SlepcContourData contour;
    1196          33 :   PetscInt         i,nsplit;
    1197          33 :   PC               pc;
    1198          33 :   MPI_Comm         child;
    1199             : 
    1200          33 :   PetscFunctionBegin;
    1201          33 :   if (!ctx->contour) {  /* initialize contour data structure first */
    1202          33 :     PetscCall(RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj));
    1203          33 :     PetscCall(SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour));
    1204             :   }
    1205          33 :   contour = ctx->contour;
    1206          33 :   if (!contour->ksp) {
    1207          33 :     PetscCall(PetscMalloc1(contour->npoints,&contour->ksp));
    1208          33 :     PetscCall(EPSGetST(eps,&eps->st));
    1209          33 :     PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
    1210          33 :     PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
    1211         841 :     for (i=0;i<contour->npoints;i++) {
    1212         808 :       PetscCall(KSPCreate(child,&contour->ksp[i]));
    1213         808 :       PetscCall(PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)eps,1));
    1214         808 :       PetscCall(KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)eps)->prefix));
    1215         808 :       PetscCall(KSPAppendOptionsPrefix(contour->ksp[i],"eps_ciss_"));
    1216         808 :       PetscCall(PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)eps)->options));
    1217         808 :       PetscCall(KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE));
    1218        1224 :       PetscCall(KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(eps->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
    1219         808 :       PetscCall(KSPGetPC(contour->ksp[i],&pc));
    1220         808 :       if (nsplit) {
    1221          96 :         PetscCall(KSPSetType(contour->ksp[i],KSPBCGS));
    1222          96 :         PetscCall(PCSetType(pc,PCBJACOBI));
    1223             :       } else {
    1224         712 :         PetscCall(KSPSetType(contour->ksp[i],KSPPREONLY));
    1225         808 :         PetscCall(PCSetType(pc,PCLU));
    1226             :       }
    1227             :     }
    1228             :   }
    1229          33 :   if (nsolve) *nsolve = contour->npoints;
    1230          33 :   if (ksp)    *ksp    = contour->ksp;
    1231          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1232             : }
    1233             : 
    1234             : /*@C
    1235             :    EPSCISSGetKSPs - Retrieve the array of linear solver objects associated with
    1236             :    the CISS solver.
    1237             : 
    1238             :    Not Collective
    1239             : 
    1240             :    Input Parameter:
    1241             : .  eps - the eigenproblem solver solver
    1242             : 
    1243             :    Output Parameters:
    1244             : +  nsolve - number of solver objects
    1245             : -  ksp - array of linear solver object
    1246             : 
    1247             :    Notes:
    1248             :    The number of KSP solvers is equal to the number of integration points divided by
    1249             :    the number of partitions. This value is halved in the case of real matrices with
    1250             :    a region centered at the real axis.
    1251             : 
    1252             :    Level: advanced
    1253             : 
    1254             : .seealso: EPSCISSSetSizes()
    1255             : @*/
    1256          33 : PetscErrorCode EPSCISSGetKSPs(EPS eps,PetscInt *nsolve,KSP **ksp)
    1257             : {
    1258          33 :   PetscFunctionBegin;
    1259          33 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
    1260          33 :   PetscUseMethod(eps,"EPSCISSGetKSPs_C",(EPS,PetscInt*,KSP**),(eps,nsolve,ksp));
    1261          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1262             : }
    1263             : 
    1264          33 : static PetscErrorCode EPSReset_CISS(EPS eps)
    1265             : {
    1266          33 :   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
    1267             : 
    1268          33 :   PetscFunctionBegin;
    1269          33 :   PetscCall(BVDestroy(&ctx->S));
    1270          33 :   PetscCall(BVDestroy(&ctx->V));
    1271          33 :   PetscCall(BVDestroy(&ctx->Y));
    1272          33 :   if (!ctx->usest) PetscCall(SlepcContourDataReset(ctx->contour));
    1273          33 :   PetscCall(BVDestroy(&ctx->pV));
    1274          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1275             : }
    1276             : 
    1277          33 : static PetscErrorCode EPSSetFromOptions_CISS(EPS eps,PetscOptionItems *PetscOptionsObject)
    1278             : {
    1279          33 :   PetscReal         r3,r4;
    1280          33 :   PetscInt          i,i1,i2,i3,i4,i5,i6,i7;
    1281          33 :   PetscBool         b1,b2,flg,flg2,flg3,flg4,flg5,flg6;
    1282          33 :   EPS_CISS          *ctx = (EPS_CISS*)eps->data;
    1283          33 :   EPSCISSQuadRule   quad;
    1284          33 :   EPSCISSExtraction extraction;
    1285             : 
    1286          33 :   PetscFunctionBegin;
    1287          33 :   PetscOptionsHeadBegin(PetscOptionsObject,"EPS CISS Options");
    1288             : 
    1289          33 :     PetscCall(EPSCISSGetSizes(eps,&i1,&i2,&i3,&i4,&i5,&b1));
    1290          33 :     PetscCall(PetscOptionsInt("-eps_ciss_integration_points","Number of integration points","EPSCISSSetSizes",i1,&i1,&flg));
    1291          33 :     PetscCall(PetscOptionsInt("-eps_ciss_blocksize","Block size","EPSCISSSetSizes",i2,&i2,&flg2));
    1292          33 :     PetscCall(PetscOptionsInt("-eps_ciss_moments","Moment size","EPSCISSSetSizes",i3,&i3,&flg3));
    1293          33 :     PetscCall(PetscOptionsInt("-eps_ciss_partitions","Number of partitions","EPSCISSSetSizes",i4,&i4,&flg4));
    1294          33 :     PetscCall(PetscOptionsInt("-eps_ciss_maxblocksize","Maximum block size","EPSCISSSetSizes",i5,&i5,&flg5));
    1295          33 :     PetscCall(PetscOptionsBool("-eps_ciss_realmats","True if A and B are real","EPSCISSSetSizes",b1,&b1,&flg6));
    1296          33 :     if (flg || flg2 || flg3 || flg4 || flg5 || flg6) PetscCall(EPSCISSSetSizes(eps,i1,i2,i3,i4,i5,b1));
    1297             : 
    1298          33 :     PetscCall(EPSCISSGetThreshold(eps,&r3,&r4));
    1299          33 :     PetscCall(PetscOptionsReal("-eps_ciss_delta","Threshold for numerical rank","EPSCISSSetThreshold",r3,&r3,&flg));
    1300          33 :     PetscCall(PetscOptionsReal("-eps_ciss_spurious_threshold","Threshold for the spurious eigenpairs","EPSCISSSetThreshold",r4,&r4,&flg2));
    1301          33 :     if (flg || flg2) PetscCall(EPSCISSSetThreshold(eps,r3,r4));
    1302             : 
    1303          33 :     PetscCall(EPSCISSGetRefinement(eps,&i6,&i7));
    1304          33 :     PetscCall(PetscOptionsInt("-eps_ciss_refine_inner","Number of inner iterative refinement iterations","EPSCISSSetRefinement",i6,&i6,&flg));
    1305          33 :     PetscCall(PetscOptionsInt("-eps_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","EPSCISSSetRefinement",i7,&i7,&flg2));
    1306          33 :     if (flg || flg2) PetscCall(EPSCISSSetRefinement(eps,i6,i7));
    1307             : 
    1308          33 :     PetscCall(EPSCISSGetUseST(eps,&b2));
    1309          33 :     PetscCall(PetscOptionsBool("-eps_ciss_usest","Use ST for linear solves","EPSCISSSetUseST",b2,&b2,&flg));
    1310          33 :     if (flg) PetscCall(EPSCISSSetUseST(eps,b2));
    1311             : 
    1312          33 :     PetscCall(PetscOptionsEnum("-eps_ciss_quadrule","Quadrature rule","EPSCISSSetQuadRule",EPSCISSQuadRules,(PetscEnum)ctx->quad,(PetscEnum*)&quad,&flg));
    1313          33 :     if (flg) PetscCall(EPSCISSSetQuadRule(eps,quad));
    1314             : 
    1315          33 :     PetscCall(PetscOptionsEnum("-eps_ciss_extraction","Extraction technique","EPSCISSSetExtraction",EPSCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg));
    1316          33 :     if (flg) PetscCall(EPSCISSSetExtraction(eps,extraction));
    1317             : 
    1318          33 :   PetscOptionsHeadEnd();
    1319             : 
    1320          33 :   if (!eps->rg) PetscCall(EPSGetRG(eps,&eps->rg));
    1321          33 :   PetscCall(RGSetFromOptions(eps->rg)); /* this is necessary here to set useconj */
    1322          33 :   if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
    1323          33 :   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
    1324         841 :   for (i=0;i<ctx->contour->npoints;i++) PetscCall(KSPSetFromOptions(ctx->contour->ksp[i]));
    1325          33 :   PetscCall(PetscSubcommSetFromOptions(ctx->contour->subcomm));
    1326          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1327             : }
    1328             : 
    1329          33 : static PetscErrorCode EPSDestroy_CISS(EPS eps)
    1330             : {
    1331          33 :   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
    1332             : 
    1333          33 :   PetscFunctionBegin;
    1334          33 :   PetscCall(SlepcContourDataDestroy(&ctx->contour));
    1335          33 :   PetscCall(PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma));
    1336          33 :   PetscCall(PetscFree(eps->data));
    1337          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",NULL));
    1338          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",NULL));
    1339          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",NULL));
    1340          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",NULL));
    1341          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",NULL));
    1342          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",NULL));
    1343          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",NULL));
    1344          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",NULL));
    1345          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",NULL));
    1346          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",NULL));
    1347          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",NULL));
    1348          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",NULL));
    1349          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",NULL));
    1350          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1351             : }
    1352             : 
    1353           0 : static PetscErrorCode EPSView_CISS(EPS eps,PetscViewer viewer)
    1354             : {
    1355           0 :   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
    1356           0 :   PetscBool      isascii;
    1357           0 :   PetscViewer    sviewer;
    1358             : 
    1359           0 :   PetscFunctionBegin;
    1360           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    1361           0 :   if (isascii) {
    1362           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  sizes { integration points: %" PetscInt_FMT ", block size: %" PetscInt_FMT ", moment size: %" PetscInt_FMT ", partitions: %" PetscInt_FMT ", maximum block size: %" PetscInt_FMT " }\n",ctx->N,ctx->L,ctx->M,ctx->npart,ctx->L_max));
    1363           0 :     if (ctx->isreal) PetscCall(PetscViewerASCIIPrintf(viewer,"  exploiting symmetry of integration points\n"));
    1364           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold));
    1365           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  iterative refinement { inner: %" PetscInt_FMT ", blocksize: %" PetscInt_FMT " }\n",ctx->refine_inner, ctx->refine_blocksize));
    1366           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  extraction: %s\n",EPSCISSExtractions[ctx->extraction]));
    1367           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  quadrature rule: %s\n",EPSCISSQuadRules[ctx->quad]));
    1368           0 :     if (ctx->usest) PetscCall(PetscViewerASCIIPrintf(viewer,"  using ST for linear solves\n"));
    1369             :     else {
    1370           0 :       if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
    1371           0 :       PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
    1372           0 :       PetscCall(PetscViewerASCIIPushTab(viewer));
    1373           0 :       if (ctx->npart>1 && ctx->contour->subcomm) {
    1374           0 :         PetscCall(PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
    1375           0 :         if (!ctx->contour->subcomm->color) PetscCall(KSPView(ctx->contour->ksp[0],sviewer));
    1376           0 :         PetscCall(PetscViewerFlush(sviewer));
    1377           0 :         PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
    1378             :         /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
    1379           0 :         PetscCall(PetscViewerASCIIPopSynchronized(viewer));
    1380           0 :       } else PetscCall(KSPView(ctx->contour->ksp[0],viewer));
    1381           0 :       PetscCall(PetscViewerASCIIPopTab(viewer));
    1382             :     }
    1383             :   }
    1384           0 :   PetscFunctionReturn(PETSC_SUCCESS);
    1385             : }
    1386             : 
    1387          66 : static PetscErrorCode EPSSetDefaultST_CISS(EPS eps)
    1388             : {
    1389          66 :   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
    1390          66 :   PetscBool      usest = ctx->usest;
    1391          66 :   KSP            ksp;
    1392          66 :   PC             pc;
    1393             : 
    1394          66 :   PetscFunctionBegin;
    1395          66 :   if (!((PetscObject)eps->st)->type_name) {
    1396          28 :     if (!ctx->usest_set) usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
    1397          28 :     if (usest) PetscCall(STSetType(eps->st,STSINVERT));
    1398             :     else {
    1399             :       /* we are not going to use ST, so avoid factorizing the matrix */
    1400          15 :       PetscCall(STSetType(eps->st,STSHIFT));
    1401          15 :       if (eps->isgeneralized) {
    1402           3 :         PetscCall(STGetKSP(eps->st,&ksp));
    1403           3 :         PetscCall(KSPGetPC(ksp,&pc));
    1404           3 :         PetscCall(PCSetType(pc,PCNONE));
    1405             :       }
    1406             :     }
    1407             :   }
    1408          66 :   PetscFunctionReturn(PETSC_SUCCESS);
    1409             : }
    1410             : 
    1411          33 : SLEPC_EXTERN PetscErrorCode EPSCreate_CISS(EPS eps)
    1412             : {
    1413          33 :   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
    1414             : 
    1415          33 :   PetscFunctionBegin;
    1416          33 :   PetscCall(PetscNew(&ctx));
    1417          33 :   eps->data = ctx;
    1418             : 
    1419          33 :   eps->useds = PETSC_TRUE;
    1420          33 :   eps->categ = EPS_CATEGORY_CONTOUR;
    1421             : 
    1422          33 :   eps->ops->solve          = EPSSolve_CISS;
    1423          33 :   eps->ops->setup          = EPSSetUp_CISS;
    1424          33 :   eps->ops->setupsort      = EPSSetUpSort_CISS;
    1425          33 :   eps->ops->setfromoptions = EPSSetFromOptions_CISS;
    1426          33 :   eps->ops->destroy        = EPSDestroy_CISS;
    1427          33 :   eps->ops->reset          = EPSReset_CISS;
    1428          33 :   eps->ops->view           = EPSView_CISS;
    1429          33 :   eps->ops->computevectors = EPSComputeVectors_CISS;
    1430          33 :   eps->ops->setdefaultst   = EPSSetDefaultST_CISS;
    1431             : 
    1432          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",EPSCISSSetSizes_CISS));
    1433          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",EPSCISSGetSizes_CISS));
    1434          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",EPSCISSSetThreshold_CISS));
    1435          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",EPSCISSGetThreshold_CISS));
    1436          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",EPSCISSSetRefinement_CISS));
    1437          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",EPSCISSGetRefinement_CISS));
    1438          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",EPSCISSSetUseST_CISS));
    1439          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",EPSCISSGetUseST_CISS));
    1440          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",EPSCISSSetQuadRule_CISS));
    1441          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",EPSCISSGetQuadRule_CISS));
    1442          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",EPSCISSSetExtraction_CISS));
    1443          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",EPSCISSGetExtraction_CISS));
    1444          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",EPSCISSGetKSPs_CISS));
    1445             : 
    1446             :   /* set default values of parameters */
    1447          33 :   ctx->N                  = 32;
    1448          33 :   ctx->L                  = 16;
    1449          33 :   ctx->M                  = ctx->N/4;
    1450          33 :   ctx->delta              = SLEPC_DEFAULT_TOL*1e-4;
    1451          33 :   ctx->L_max              = 64;
    1452          33 :   ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
    1453          33 :   ctx->usest              = PETSC_TRUE;
    1454          33 :   ctx->usest_set          = PETSC_FALSE;
    1455          33 :   ctx->isreal             = PETSC_FALSE;
    1456          33 :   ctx->refine_inner       = 0;
    1457          33 :   ctx->refine_blocksize   = 0;
    1458          33 :   ctx->npart              = 1;
    1459          33 :   ctx->quad               = (EPSCISSQuadRule)0;
    1460          33 :   ctx->extraction         = EPS_CISS_EXTRACTION_RITZ;
    1461          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1462             : }

Generated by: LCOV version 1.14