LCOV - code coverage report
Current view: top level - eps/impls/ciss - ciss.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 779 871 89.4 %
Date: 2021-08-02 00:32:28 Functions: 39 40 97.5 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14