Actual source code: ciss.c

slepc-3.13.4 2020-09-02
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc eigensolver: "ciss"

 13:    Method: Contour Integral Spectral Slicing

 15:    Algorithm:

 17:        Contour integral based on Sakurai-Sugiura method to construct a
 18:        subspace, with various eigenpair extractions (Rayleigh-Ritz,
 19:        explicit moment).

 21:    Based on code contributed by Y. Maeda, T. Sakurai.

 23:    References:

 25:        [1] T. Sakurai and H. Sugiura, "A projection method for generalized
 26:            eigenvalue problems", J. Comput. Appl. Math. 159:119-128, 2003.

 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: */

 33:  #include <slepc/private/epsimpl.h>
 34:  #include <slepcblaslapack.h>

 36: typedef struct {
 37:   /* parameters */
 38:   PetscInt          N;          /* number of integration points (32) */
 39:   PetscInt          L;          /* block size (16) */
 40:   PetscInt          M;          /* moment degree (N/4 = 4) */
 41:   PetscReal         delta;      /* threshold of singular value (1e-12) */
 42:   PetscInt          L_max;      /* maximum number of columns of the source matrix V */
 43:   PetscReal         spurious_threshold; /* discard spurious eigenpairs */
 44:   PetscBool         isreal;     /* A and B are real */
 45:   PetscInt          npart;      /* number of partitions */
 46:   PetscInt          refine_inner;
 47:   PetscInt          refine_blocksize;
 48:   /* private data */
 49:   PetscReal         *sigma;     /* threshold for numerical rank */
 50:   PetscInt          subcomm_id;
 51:   PetscInt          num_solve_point;
 52:   PetscScalar       *weight;
 53:   PetscScalar       *omega;
 54:   PetscScalar       *pp;
 55:   BV                V;
 56:   BV                S;
 57:   BV                pV;
 58:   BV                Y;
 59:   Vec               xsub;
 60:   Vec               xdup;
 61:   KSP               *ksp;       /* ksp array for storing factorizations at integration points */
 62:   PetscBool         useconj;
 63:   PetscReal         est_eig;
 64:   VecScatter        scatterin;
 65:   Mat               pA,pB;
 66:   PetscSubcomm      subcomm;
 67:   PetscBool         usest;
 68:   PetscBool         usest_set;  /* whether the user set the usest flag or not */
 69:   EPSCISSQuadRule   quad;
 70:   EPSCISSExtraction extraction;
 71: } EPS_CISS;

 73: /* destroy KSP objects when the number of solve points changes */
 74: PETSC_STATIC_INLINE PetscErrorCode EPSCISSResetSolvers(EPS eps)
 75: {
 77:   PetscInt       i;
 78:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

 81:   if (ctx->ksp) {
 82:     for (i=0;i<ctx->num_solve_point;i++) {
 83:       KSPDestroy(&ctx->ksp[i]);
 84:     }
 85:     PetscFree(ctx->ksp);
 86:   }
 87:   return(0);
 88: }

 90: /* clean PetscSubcomm object when the number of partitions changes */
 91: PETSC_STATIC_INLINE PetscErrorCode EPSCISSResetSubcomm(EPS eps)
 92: {
 94:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

 97:   EPSCISSResetSolvers(eps);
 98:   PetscSubcommDestroy(&ctx->subcomm);
 99:   return(0);
100: }

102: /* determine whether half of integration points can be avoided (use its conjugates);
103:    depends on isreal and the center of the region */
104: PETSC_STATIC_INLINE PetscErrorCode EPSCISSSetUseConj(EPS eps,PetscBool *useconj)
105: {
107:   PetscScalar    center;
108:   PetscReal      c,d;
109:   PetscBool      isellipse,isinterval;
110: #if defined(PETSC_USE_COMPLEX)
111:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
112: #endif

115:   *useconj = PETSC_FALSE;
116:   PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
117:   if (isellipse) {
118:     RGEllipseGetParameters(eps->rg,&center,NULL,NULL);
119: #if defined(PETSC_USE_COMPLEX)
120:     *useconj = (ctx->isreal && PetscImaginaryPart(center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
121: #endif
122:   } else {
123:     PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
124:     if (isinterval) {
125:       RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
126: #if defined(PETSC_USE_COMPLEX)
127:       *useconj = (ctx->isreal && c==d)? PETSC_TRUE: PETSC_FALSE;
128: #endif
129:     }
130:   }
131:   return(0);
132: }

134: /* create PetscSubcomm object and determine num_solve_point (depends on npart, N, useconj) */
135: PETSC_STATIC_INLINE PetscErrorCode EPSCISSSetUpSubComm(EPS eps,PetscInt *num_solve_point)
136: {
138:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
139:   PetscInt       N = ctx->N;

142:   PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subcomm);
143:   PetscSubcommSetNumber(ctx->subcomm,ctx->npart);
144:   PetscSubcommSetType(ctx->subcomm,PETSC_SUBCOMM_INTERLACED);
145:   PetscLogObjectMemory((PetscObject)eps,sizeof(PetscSubcomm));
146:   ctx->subcomm_id = ctx->subcomm->color;
147:   EPSCISSSetUseConj(eps,&ctx->useconj);
148:   if (ctx->useconj) N = N/2;
149:   *num_solve_point = N / ctx->npart;
150:   if (N%ctx->npart > ctx->subcomm_id) (*num_solve_point)++;
151:   return(0);
152: }

154: static PetscErrorCode CISSRedundantMat(EPS eps)
155: {
157:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
158:   Mat            A,B;
159:   PetscInt       nmat;

162:   STGetNumMatrices(eps->st,&nmat);
163:   if (ctx->subcomm->n != 1) {
164:     STGetMatrix(eps->st,0,&A);
165:     MatDestroy(&ctx->pA);
166:     MatCreateRedundantMatrix(A,ctx->subcomm->n,PetscSubcommChild(ctx->subcomm),MAT_INITIAL_MATRIX,&ctx->pA);
167:     if (nmat>1) {
168:       STGetMatrix(eps->st,1,&B);
169:       MatDestroy(&ctx->pB);
170:       MatCreateRedundantMatrix(B,ctx->subcomm->n,PetscSubcommChild(ctx->subcomm),MAT_INITIAL_MATRIX,&ctx->pB);
171:     } else ctx->pB = NULL;
172:   } else {
173:     ctx->pA = NULL;
174:     ctx->pB = NULL;
175:   }
176:   return(0);
177: }

179: static PetscErrorCode CISSScatterVec(EPS eps)
180: {
182:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
183:   IS             is1,is2;
184:   Vec            v0;
185:   PetscInt       i,j,k,mstart,mend,mlocal;
186:   PetscInt       *idx1,*idx2,mloc_sub;

189:   VecDestroy(&ctx->xsub);
190:   MatCreateVecs(ctx->pA,&ctx->xsub,NULL);

192:   VecDestroy(&ctx->xdup);
193:   MatGetLocalSize(ctx->pA,&mloc_sub,NULL);
194:   VecCreateMPI(PetscSubcommContiguousParent(ctx->subcomm),mloc_sub,PETSC_DECIDE,&ctx->xdup);

196:   VecScatterDestroy(&ctx->scatterin);
197:   BVGetColumn(ctx->V,0,&v0);
198:   VecGetOwnershipRange(v0,&mstart,&mend);
199:   mlocal = mend - mstart;
200:   PetscMalloc2(ctx->subcomm->n*mlocal,&idx1,ctx->subcomm->n*mlocal,&idx2);
201:   j = 0;
202:   for (k=0;k<ctx->subcomm->n;k++) {
203:     for (i=mstart;i<mend;i++) {
204:       idx1[j]   = i;
205:       idx2[j++] = i + eps->n*k;
206:     }
207:   }
208:   ISCreateGeneral(PetscObjectComm((PetscObject)eps),ctx->subcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);
209:   ISCreateGeneral(PetscObjectComm((PetscObject)eps),ctx->subcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);
210:   VecScatterCreate(v0,is1,ctx->xdup,is2,&ctx->scatterin);
211:   ISDestroy(&is1);
212:   ISDestroy(&is2);
213:   PetscFree2(idx1,idx2);
214:   BVRestoreColumn(ctx->V,0,&v0);
215:   return(0);
216: }

218: static PetscErrorCode SetPathParameter(EPS eps)
219: {
221:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
222:   PetscInt       i,j;
223:   PetscScalar    center=0.0,tmp,tmp2,*omegai;
224:   PetscReal      theta,radius=1.0,vscale,a,b,c,d,max_w=0.0,rgscale;
225: #if defined(PETSC_USE_COMPLEX)
226:   PetscReal      start_ang,end_ang;
227: #endif
228:   PetscBool      isring=PETSC_FALSE,isellipse=PETSC_FALSE,isinterval=PETSC_FALSE;

231:   PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
232:   PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
233:   PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
234:   RGGetScale(eps->rg,&rgscale);
235:   PetscMalloc1(ctx->N+1l,&omegai);
236:   RGComputeContour(eps->rg,ctx->N,ctx->omega,omegai);
237:   if (isellipse) {
238:     RGEllipseGetParameters(eps->rg,&center,&radius,&vscale);
239:     for (i=0;i<ctx->N;i++) {
240: #if defined(PETSC_USE_COMPLEX)
241:       theta = 2.0*PETSC_PI*(i+0.5)/ctx->N;
242:       ctx->pp[i] = PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
243:       ctx->weight[i] = rgscale*radius*(PetscCMPLX(vscale*PetscCosReal(theta),PetscSinReal(theta)))/(PetscReal)ctx->N;
244: #else
245:       theta = (PETSC_PI/ctx->N)*(i+0.5);
246:       ctx->pp[i] = PetscCosReal(theta);
247:       ctx->weight[i] = PetscCosReal((ctx->N-1)*theta)/ctx->N;
248:       ctx->omega[i] = rgscale*(center + radius*ctx->pp[i]);
249: #endif
250:     }
251:   } else if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
252:     for (i=0;i<ctx->N;i++) {
253:       theta = (PETSC_PI/ctx->N)*(i+0.5);
254:       ctx->pp[i] = PetscCosReal(theta);
255:       ctx->weight[i] = PetscCosReal((ctx->N-1)*theta)/ctx->N;
256:     }
257:     if (isinterval) {
258:       RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d);
259:       if ((c!=d || c!=0.0) && (a!=b || a!=0.0)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Endpoints of the imaginary axis or the real axis must be both zero");
260:       for (i=0;i<ctx->N;i++) {
261:         if (c==d) ctx->omega[i] = ((b-a)*(ctx->pp[i]+1.0)/2.0+a)*rgscale;
262:         if (a==b) {
263: #if defined(PETSC_USE_COMPLEX)
264:           ctx->omega[i] = ((d-c)*(ctx->pp[i]+1.0)/2.0+c)*rgscale*PETSC_i;
265: #else
266:           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
267: #endif
268:         }
269:       }
270:     }
271:     if (isring) {  /* only supported in complex scalars */
272: #if defined(PETSC_USE_COMPLEX)
273:       RGRingGetParameters(eps->rg,&center,&radius,&vscale,&start_ang,&end_ang,NULL);
274:       for (i=0;i<ctx->N;i++) {
275:         theta = (start_ang*2.0+(end_ang-start_ang)*(PetscRealPart(ctx->pp[i])+1.0))*PETSC_PI;
276:         ctx->omega[i] = rgscale*(center + radius*PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta)));
277:       }
278: #endif
279:     }
280:   } else {
281:     if (isinterval) {
282:       RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d);
283:       center = rgscale*((b+a)/2.0+(d+c)/2.0*PETSC_PI);
284:       radius = PetscSqrtReal(PetscPowRealInt(rgscale*(b-a)/2.0,2)+PetscPowRealInt(rgscale*(d-c)/2.0,2));
285:     } else if (isring) {
286:       RGRingGetParameters(eps->rg,&center,&radius,NULL,NULL,NULL,NULL);
287:       center *= rgscale;
288:       radius *= rgscale;
289:     }
290:     for (i=0;i<ctx->N;i++) {
291:       ctx->pp[i] = (ctx->omega[i]-center)/radius;
292:       tmp = 1; tmp2 = 1;
293:       for (j=0;j<ctx->N;j++) {
294:         tmp *= ctx->omega[j];
295:         if (i != j) tmp2 *= ctx->omega[j]-ctx->omega[i];
296:       }
297:       ctx->weight[i] = tmp/tmp2;
298:       max_w = PetscMax(PetscAbsScalar(ctx->weight[i]),max_w);
299:     }
300:     for (i=0;i<ctx->N;i++) ctx->weight[i] /= (PetscScalar)max_w;
301:   }
302:   PetscFree(omegai);
303:   return(0);
304: }

306: static PetscErrorCode CISSVecSetRandom(BV V,PetscInt i0,PetscInt i1)
307: {
309:   PetscInt       i,j,nlocal;
310:   PetscScalar    *vdata;
311:   Vec            x;

314:   BVGetSizes(V,&nlocal,NULL,NULL);
315:   for (i=i0;i<i1;i++) {
316:     BVSetRandomColumn(V,i);
317:     BVGetColumn(V,i,&x);
318:     VecGetArray(x,&vdata);
319:     for (j=0;j<nlocal;j++) {
320:       vdata[j] = PetscRealPart(vdata[j]);
321:       if (PetscRealPart(vdata[j]) < 0.5) vdata[j] = -1.0;
322:       else vdata[j] = 1.0;
323:     }
324:     VecRestoreArray(x,&vdata);
325:     BVRestoreColumn(V,i,&x);
326:   }
327:   return(0);
328: }

330: static PetscErrorCode VecScatterVecs(EPS eps,BV Vin,PetscInt n)
331: {
332:   PetscErrorCode    ierr;
333:   EPS_CISS          *ctx = (EPS_CISS*)eps->data;
334:   PetscInt          i;
335:   Vec               vi,pvi;
336:   const PetscScalar *array;

339:   for (i=0;i<n;i++) {
340:     BVGetColumn(Vin,i,&vi);
341:     VecScatterBegin(ctx->scatterin,vi,ctx->xdup,INSERT_VALUES,SCATTER_FORWARD);
342:     VecScatterEnd(ctx->scatterin,vi,ctx->xdup,INSERT_VALUES,SCATTER_FORWARD);
343:     BVRestoreColumn(Vin,i,&vi);
344:     VecGetArrayRead(ctx->xdup,&array);
345:     VecPlaceArray(ctx->xsub,array);
346:     BVGetColumn(ctx->pV,i,&pvi);
347:     VecCopy(ctx->xsub,pvi);
348:     BVRestoreColumn(ctx->pV,i,&pvi);
349:     VecResetArray(ctx->xsub);
350:     VecRestoreArrayRead(ctx->xdup,&array);
351:   }
352:   return(0);
353: }

355: static PetscErrorCode SolveLinearSystem(EPS eps,Mat A,Mat B,BV V,PetscInt L_start,PetscInt L_end,PetscBool initksp)
356: {
358:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
359:   PetscInt       i,j,p_id;
360:   Mat            Fz,kspMat;
361:   Vec            Bvj,vj,yj;
362:   KSP            ksp;

365:   if (!ctx->ksp) { EPSCISSGetKSPs(eps,&ctx->num_solve_point,&ctx->ksp); }
366:   BVCreateVec(V,&Bvj);
367:   if (ctx->usest) {
368:     MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Fz);
369:   }
370:   for (i=0;i<ctx->num_solve_point;i++) {
371:     p_id = i*ctx->subcomm->n + ctx->subcomm_id;
372:     if (!ctx->usest && initksp) {
373:       MatDuplicate(A,MAT_COPY_VALUES,&kspMat);
374:       if (B) {
375:         MatAXPY(kspMat,-ctx->omega[p_id],B,DIFFERENT_NONZERO_PATTERN);
376:       } else {
377:         MatShift(kspMat,-ctx->omega[p_id]);
378:       }
379:       KSPSetOperators(ctx->ksp[i],kspMat,kspMat);
380:       MatDestroy(&kspMat);
381:     } else if (ctx->usest) {
382:       STSetShift(eps->st,ctx->omega[p_id]);
383:       STGetKSP(eps->st,&ksp);
384:     }
385:     for (j=L_start;j<L_end;j++) {
386:       BVGetColumn(V,j,&vj);
387:       BVGetColumn(ctx->Y,i*ctx->L_max+j,&yj);
388:       if (B) {
389:         MatMult(B,vj,Bvj);
390:         if (ctx->usest) {
391:           KSPSolve(ksp,Bvj,yj);
392:         } else {
393:           KSPSolve(ctx->ksp[i],Bvj,yj);
394:         }
395:       } else {
396:         if (ctx->usest) {
397:           KSPSolve(ksp,vj,yj);
398:         } else {
399:           KSPSolve(ctx->ksp[i],vj,yj);
400:         }
401:       }
402:       BVRestoreColumn(V,j,&vj);
403:       BVRestoreColumn(ctx->Y,i*ctx->L_max+j,&yj);
404:     }
405:     if (ctx->usest && i<ctx->num_solve_point-1) { KSPReset(ksp); }
406:   }
407:   if (ctx->usest) { MatDestroy(&Fz); }
408:   VecDestroy(&Bvj);
409:   return(0);
410: }

412: #if defined(PETSC_USE_COMPLEX)
413: static PetscErrorCode EstimateNumberEigs(EPS eps,PetscInt *L_add)
414: {
416:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
417:   PetscInt       i,j,p_id;
418:   PetscScalar    tmp,m = 1,sum = 0.0;
419:   PetscReal      eta;
420:   Vec            v,vtemp,vj;

423:   BVCreateVec(ctx->Y,&v);
424:   BVCreateVec(ctx->V,&vtemp);
425:   for (j=0;j<ctx->L;j++) {
426:     VecSet(v,0);
427:     for (i=0;i<ctx->num_solve_point; i++) {
428:       p_id = i*ctx->subcomm->n + ctx->subcomm_id;
429:       BVSetActiveColumns(ctx->Y,i*ctx->L_max+j,i*ctx->L_max+j+1);
430:       BVMultVec(ctx->Y,ctx->weight[p_id],1,v,&m);
431:     }
432:     BVGetColumn(ctx->V,j,&vj);
433:     if (ctx->pA) {
434:       VecSet(vtemp,0);
435:       VecScatterBegin(ctx->scatterin,v,vtemp,ADD_VALUES,SCATTER_REVERSE);
436:       VecScatterEnd(ctx->scatterin,v,vtemp,ADD_VALUES,SCATTER_REVERSE);
437:       VecDot(vj,vtemp,&tmp);
438:     } else {
439:       VecDot(vj,v,&tmp);
440:     }
441:     BVRestoreColumn(ctx->V,j,&vj);
442:     if (ctx->useconj) sum += PetscRealPart(tmp)*2;
443:     else sum += tmp;
444:   }
445:   ctx->est_eig = PetscAbsScalar(sum/(PetscReal)ctx->L);
446:   eta = PetscPowReal(10.0,-PetscLog10Real(eps->tol)/ctx->N);
447:   PetscInfo1(eps,"Estimation_#Eig %f\n",(double)ctx->est_eig);
448:   *L_add = (PetscInt)PetscCeilReal((ctx->est_eig*eta)/ctx->M) - ctx->L;
449:   if (*L_add < 0) *L_add = 0;
450:   if (*L_add>ctx->L_max-ctx->L) {
451:     PetscInfo(eps,"Number of eigenvalues around the contour path may be too large\n");
452:     *L_add = ctx->L_max-ctx->L;
453:   }
454:   VecDestroy(&v);
455:   VecDestroy(&vtemp);
456:   return(0);
457: }
458: #endif

460: static PetscErrorCode CalcMu(EPS eps,PetscScalar *Mu)
461: {
463:   PetscMPIInt    sub_size,len;
464:   PetscInt       i,j,k,s;
465:   PetscScalar    *m,*temp,*temp2,*ppk,alp;
466:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
467:   Mat            M;

470:   MPI_Comm_size(PetscSubcommChild(ctx->subcomm),&sub_size);
471:   PetscMalloc3(ctx->num_solve_point*ctx->L*(ctx->L+1),&temp,2*ctx->M*ctx->L*ctx->L,&temp2,ctx->num_solve_point,&ppk);
472:   MatCreateSeqDense(PETSC_COMM_SELF,ctx->L,ctx->L_max*ctx->num_solve_point,NULL,&M);
473:   for (i=0;i<2*ctx->M*ctx->L*ctx->L;i++) temp2[i] = 0;
474:   BVSetActiveColumns(ctx->Y,0,ctx->L_max*ctx->num_solve_point);
475:   if (ctx->pA) {
476:     BVSetActiveColumns(ctx->pV,0,ctx->L);
477:     BVDot(ctx->Y,ctx->pV,M);
478:   } else {
479:     BVSetActiveColumns(ctx->V,0,ctx->L);
480:     BVDot(ctx->Y,ctx->V,M);
481:   }
482:   MatDenseGetArray(M,&m);
483:   for (i=0;i<ctx->num_solve_point;i++) {
484:     for (j=0;j<ctx->L;j++) {
485:       for (k=0;k<ctx->L;k++) {
486:         temp[k+j*ctx->L+i*ctx->L*ctx->L]=m[k+j*ctx->L+i*ctx->L*ctx->L_max];
487:       }
488:     }
489:   }
490:   MatDenseRestoreArray(M,&m);
491:   for (i=0;i<ctx->num_solve_point;i++) ppk[i] = 1;
492:   for (k=0;k<2*ctx->M;k++) {
493:     for (j=0;j<ctx->L;j++) {
494:       for (i=0;i<ctx->num_solve_point;i++) {
495:         alp = ppk[i]*ctx->weight[i*ctx->subcomm->n + ctx->subcomm_id];
496:         for (s=0;s<ctx->L;s++) {
497:           if (ctx->useconj) temp2[s+(j+k*ctx->L)*ctx->L] += PetscRealPart(alp*temp[s+(j+i*ctx->L)*ctx->L])*2;
498:           else temp2[s+(j+k*ctx->L)*ctx->L] += alp*temp[s+(j+i*ctx->L)*ctx->L];
499:         }
500:       }
501:     }
502:     for (i=0;i<ctx->num_solve_point;i++)
503:       ppk[i] *= ctx->pp[i*ctx->subcomm->n + ctx->subcomm_id];
504:   }
505:   for (i=0;i<2*ctx->M*ctx->L*ctx->L;i++) temp2[i] /= sub_size;
506:   PetscMPIIntCast(2*ctx->M*ctx->L*ctx->L,&len);
507:   MPI_Allreduce(temp2,Mu,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)eps));
508:   PetscFree3(temp,temp2,ppk);
509:   MatDestroy(&M);
510:   return(0);
511: }

513: static PetscErrorCode BlockHankel(EPS eps,PetscScalar *Mu,PetscInt s,PetscScalar *H)
514: {
515:   EPS_CISS *ctx = (EPS_CISS*)eps->data;
516:   PetscInt i,j,k,L=ctx->L,M=ctx->M;

519:   for (k=0;k<L*M;k++)
520:     for (j=0;j<M;j++)
521:       for (i=0;i<L;i++)
522:         H[j*L+i+k*L*M] = Mu[i+k*L+(j+s)*L*L];
523:   return(0);
524: }

526: static PetscErrorCode SVD_H0(EPS eps,PetscScalar *S,PetscInt *K)
527: {
529:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
530:   PetscInt       i,ml=ctx->L*ctx->M;
531:   PetscBLASInt   m,n,lda,ldu,ldvt,lwork,info;
532:   PetscScalar    *work;
533: #if defined(PETSC_USE_COMPLEX)
534:   PetscReal      *rwork;
535: #endif

538:   PetscMalloc1(5*ml,&work);
539: #if defined(PETSC_USE_COMPLEX)
540:   PetscMalloc1(5*ml,&rwork);
541: #endif
542:   PetscBLASIntCast(ml,&m);
543:   n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m;
544:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
545: #if defined(PETSC_USE_COMPLEX)
546:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
547: #else
548:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
549: #endif
550:   SlepcCheckLapackInfo("gesvd",info);
551:   PetscFPTrapPop();
552:   (*K) = 0;
553:   for (i=0;i<ml;i++) {
554:     if (ctx->sigma[i]/PetscMax(ctx->sigma[0],1)>ctx->delta) (*K)++;
555:   }
556:   PetscFree(work);
557: #if defined(PETSC_USE_COMPLEX)
558:   PetscFree(rwork);
559: #endif
560:   return(0);
561: }

563: static PetscErrorCode ConstructS(EPS eps)
564: {
566:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
567:   PetscInt       i,j,k,vec_local_size,p_id;
568:   Vec            v,sj;
569:   PetscScalar    *ppk, *v_data, m = 1;

572:   BVGetSizes(ctx->Y,&vec_local_size,NULL,NULL);
573:   PetscMalloc1(ctx->num_solve_point,&ppk);
574:   for (i=0;i<ctx->num_solve_point;i++) ppk[i] = 1;
575:   BVCreateVec(ctx->Y,&v);
576:   for (k=0;k<ctx->M;k++) {
577:     for (j=0;j<ctx->L;j++) {
578:       VecSet(v,0);
579:       for (i=0;i<ctx->num_solve_point;i++) {
580:         p_id = i*ctx->subcomm->n + ctx->subcomm_id;
581:         BVSetActiveColumns(ctx->Y,i*ctx->L_max+j,i*ctx->L_max+j+1);
582:         BVMultVec(ctx->Y,ppk[i]*ctx->weight[p_id],1.0,v,&m);
583:       }
584:       if (ctx->useconj) {
585:         VecGetArray(v,&v_data);
586:         for (i=0;i<vec_local_size;i++) v_data[i] = PetscRealPart(v_data[i])*2;
587:         VecRestoreArray(v,&v_data);
588:       }
589:       BVGetColumn(ctx->S,k*ctx->L+j,&sj);
590:       if (ctx->pA) {
591:         VecSet(sj,0);
592:         VecScatterBegin(ctx->scatterin,v,sj,ADD_VALUES,SCATTER_REVERSE);
593:         VecScatterEnd(ctx->scatterin,v,sj,ADD_VALUES,SCATTER_REVERSE);
594:       } else {
595:         VecCopy(v,sj);
596:       }
597:       BVRestoreColumn(ctx->S,k*ctx->L+j,&sj);
598:     }
599:     for (i=0;i<ctx->num_solve_point;i++) {
600:       p_id = i*ctx->subcomm->n + ctx->subcomm_id;
601:       ppk[i] *= ctx->pp[p_id];
602:     }
603:   }
604:   PetscFree(ppk);
605:   VecDestroy(&v);
606:   return(0);
607: }

609: static PetscErrorCode SVD_S(BV S,PetscInt ml,PetscReal delta,PetscReal *sigma,PetscInt *K)
610: {
612:   PetscInt       i,j,k,local_size;
613:   PetscMPIInt    len;
614:   PetscScalar    *work,*temp,*B,*tempB,*s_data,*Q1,*Q2,*temp2,alpha=1,beta=0;
615:   PetscBLASInt   l,m,n,lda,ldu,ldvt,lwork,info,ldb,ldc;
616: #if defined(PETSC_USE_COMPLEX)
617:   PetscReal      *rwork;
618: #endif

621:   BVGetSizes(S,&local_size,NULL,NULL);
622:   BVGetArray(S,&s_data);
623:   PetscMalloc7(ml*ml,&temp,ml*ml,&temp2,local_size*ml,&Q1,local_size*ml,&Q2,ml*ml,&B,ml*ml,&tempB,5*ml,&work);
624:   PetscArrayzero(B,ml*ml);
625: #if defined(PETSC_USE_COMPLEX)
626:   PetscMalloc1(5*ml,&rwork);
627: #endif
628:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);

630:   for (i=0;i<ml;i++) B[i*ml+i]=1;

632:   for (k=0;k<2;k++) {
633:     PetscBLASIntCast(local_size,&m);
634:     PetscBLASIntCast(ml,&l);
635:     n = l; lda = m; ldb = m; ldc = l;
636:     if (k == 0) {
637:       PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,s_data,&lda,s_data,&ldb,&beta,temp,&ldc));
638:     } else if ((k%2)==1) {
639:       PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,Q1,&lda,Q1,&ldb,&beta,temp,&ldc));
640:     } else {
641:       PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,Q2,&lda,Q2,&ldb,&beta,temp,&ldc));
642:     }
643:     PetscArrayzero(temp2,ml*ml);
644:     PetscMPIIntCast(ml*ml,&len);
645:     MPI_Allreduce(temp,temp2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)S));

647:     PetscBLASIntCast(ml,&m);
648:     n = m; lda = m; lwork = 5*m, ldu = 1; ldvt = 1;
649: #if defined(PETSC_USE_COMPLEX)
650:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
651: #else
652:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
653: #endif
654:     SlepcCheckLapackInfo("gesvd",info);

656:     PetscBLASIntCast(local_size,&l);
657:     PetscBLASIntCast(ml,&n);
658:     m = n; lda = l; ldb = m; ldc = l;
659:     if (k==0) {
660:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,s_data,&lda,temp2,&ldb,&beta,Q1,&ldc));
661:     } else if ((k%2)==1) {
662:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,Q1,&lda,temp2,&ldb,&beta,Q2,&ldc));
663:     } else {
664:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,Q2,&lda,temp2,&ldb,&beta,Q1,&ldc));
665:     }

667:     PetscBLASIntCast(ml,&l);
668:     m = l; n = l; lda = l; ldb = m; ldc = l;
669:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,B,&lda,temp2,&ldb,&beta,tempB,&ldc));
670:     for (i=0;i<ml;i++) {
671:       sigma[i] = sqrt(sigma[i]);
672:       for (j=0;j<local_size;j++) {
673:         if ((k%2)==1) Q2[j+i*local_size]/=sigma[i];
674:         else Q1[j+i*local_size]/=sigma[i];
675:       }
676:       for (j=0;j<ml;j++) {
677:         B[j+i*ml]=tempB[j+i*ml]*sigma[i];
678:       }
679:     }
680:   }

682:   PetscBLASIntCast(ml,&m);
683:   n = m; lda = m; ldu=1; ldvt=1;
684: #if defined(PETSC_USE_COMPLEX)
685:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
686: #else
687:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
688: #endif
689:   SlepcCheckLapackInfo("gesvd",info);

691:   PetscBLASIntCast(local_size,&l);
692:   PetscBLASIntCast(ml,&n);
693:   m = n; lda = l; ldb = m; ldc = l;
694:   if ((k%2)==1) {
695:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q1,&lda,B,&ldb,&beta,s_data,&ldc));
696:   } else {
697:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q2,&lda,B,&ldb,&beta,s_data,&ldc));
698:   }

700:   PetscFPTrapPop();
701:   BVRestoreArray(S,&s_data);

703:   (*K) = 0;
704:   for (i=0;i<ml;i++) {
705:     if (sigma[i]/PetscMax(sigma[0],1)>delta) (*K)++;
706:   }
707:   PetscFree7(temp,temp2,Q1,Q2,B,tempB,work);
708: #if defined(PETSC_USE_COMPLEX)
709:   PetscFree(rwork);
710: #endif
711:   return(0);
712: }

714: static PetscErrorCode isGhost(EPS eps,PetscInt ld,PetscInt nv,PetscBool *fl)
715: {
717:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
718:   PetscInt       i,j;
719:   PetscScalar    *pX;
720:   PetscReal      *tau,s1,s2,tau_max=0.0;

723:   PetscMalloc1(nv,&tau);
724:   DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
725:   DSGetArray(eps->ds,DS_MAT_X,&pX);

727:   for (i=0;i<nv;i++) {
728:     s1 = 0;
729:     s2 = 0;
730:     for (j=0;j<nv;j++) {
731:       s1 += PetscAbsScalar(PetscPowScalarInt(pX[i*ld+j],2));
732:       s2 += PetscPowRealInt(PetscAbsScalar(pX[i*ld+j]),2)/ctx->sigma[j];
733:     }
734:     tau[i] = s1/s2;
735:     tau_max = PetscMax(tau_max,tau[i]);
736:   }
737:   DSRestoreArray(eps->ds,DS_MAT_X,&pX);
738:   for (i=0;i<nv;i++) {
739:     tau[i] /= tau_max;
740:   }
741:   for (i=0;i<nv;i++) {
742:     if (tau[i]>=ctx->spurious_threshold) fl[i] = PETSC_TRUE;
743:     else fl[i] = PETSC_FALSE;
744:   }
745:   PetscFree(tau);
746:   return(0);
747: }

749: static PetscErrorCode rescale_eig(EPS eps,PetscInt nv)
750: {
752:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
753:   PetscInt       i;
754:   PetscScalar    center;
755:   PetscReal      radius,a,b,c,d,rgscale;
756: #if defined(PETSC_USE_COMPLEX)
757:   PetscReal      start_ang,end_ang,vscale,theta;
758: #endif
759:   PetscBool      isring,isellipse,isinterval;

762:   PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
763:   PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
764:   PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
765:   RGGetScale(eps->rg,&rgscale);
766:   if (isinterval) {
767:     RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
768:     if (c==d) {
769:       for (i=0;i<nv;i++) {
770: #if defined(PETSC_USE_COMPLEX)
771:         eps->eigr[i] = PetscRealPart(eps->eigr[i]);
772: #else
773:         eps->eigi[i] = 0;
774: #endif
775:       }
776:     }
777:   }
778:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
779:     if (isellipse) {
780:       RGEllipseGetParameters(eps->rg,&center,&radius,NULL);
781:       for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
782:     } else if (isinterval) {
783:       RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d);
784:       if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
785:         for (i=0;i<nv;i++) {
786:           if (c==d) eps->eigr[i] = ((b-a)*(eps->eigr[i]+1.0)/2.0+a)*rgscale;
787:           if (a==b) {
788: #if defined(PETSC_USE_COMPLEX)
789:             eps->eigr[i] = ((d-c)*(eps->eigr[i]+1.0)/2.0+c)*rgscale*PETSC_i;
790: #else
791:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
792: #endif
793:           }
794:         }
795:       } else {
796:         center = (b+a)/2.0+(d+c)/2.0*PETSC_PI;
797:         radius = PetscSqrtReal(PetscPowRealInt((b-a)/2.0,2)+PetscPowRealInt((d-c)/2.0,2));
798:         for (i=0;i<nv;i++) eps->eigr[i] = center + radius*eps->eigr[i];
799:       }
800:     } else if (isring) {  /* only supported in complex scalars */
801: #if defined(PETSC_USE_COMPLEX)
802:       RGRingGetParameters(eps->rg,&center,&radius,&vscale,&start_ang,&end_ang,NULL);
803:       if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
804:         for (i=0;i<nv;i++) {
805:           theta = (start_ang*2.0+(end_ang-start_ang)*(PetscRealPart(eps->eigr[i])+1.0))*PETSC_PI;
806:           eps->eigr[i] = rgscale*center + (rgscale*radius+PetscImaginaryPart(eps->eigr[i]))*PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
807:         }
808:       } else {
809:         for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
810:       }
811: #endif
812:     }
813:   }
814:   return(0);
815: }

817: PetscErrorCode EPSSetUp_CISS(EPS eps)
818: {
820:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
821:   PetscBool      issinvert,istrivial,isring,isellipse,isinterval,flg,useconj;
822:   PetscReal      c,d;
823:   Mat            A;

826:   if (eps->ncv==PETSC_DEFAULT) {
827:     eps->ncv = ctx->L_max*ctx->M;
828:     if (eps->ncv>eps->n) {
829:       eps->ncv = eps->n;
830:       ctx->L_max = eps->ncv/ctx->M;
831:       if (!ctx->L_max) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot adjust solver parameters, try setting a smaller value of M (moment size)");
832:     }
833:   } else {
834:     EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
835:     ctx->L_max = eps->ncv/ctx->M;
836:     if (!ctx->L_max) {
837:       ctx->L_max = 1;
838:       eps->ncv = ctx->L_max*ctx->M;
839:     }
840:   }
841:   ctx->L = PetscMin(ctx->L,ctx->L_max);
842:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = 1;
843:   if (eps->mpd==PETSC_DEFAULT) eps->mpd = eps->ncv;
844:   if (!eps->which) eps->which = EPS_ALL;
845:   if (!eps->extraction) { EPSSetExtraction(eps,EPS_RITZ); }
846:   else if (eps->extraction!=EPS_RITZ) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
847:   if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
848:   if (eps->stopping!=EPSStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support user-defined stopping test");

850:   /* check region */
851:   RGIsTrivial(eps->rg,&istrivial);
852:   if (istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
853:   RGGetComplement(eps->rg,&flg);
854:   if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
855:   PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
856:   PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
857:   PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
858:   if (!isellipse && !isring && !isinterval) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Currently only implemented for interval, elliptic or ring regions");
859:   /* if useconj has changed, then reset subcomm data */
860:   EPSCISSSetUseConj(eps,&useconj);
861:   if (useconj!=ctx->useconj) { EPSCISSResetSubcomm(eps); }

863: #if !defined(PETSC_USE_COMPLEX)
864:   if (isring) {
865:     SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Ring region only supported for complex scalars");
866:   }
867: #endif
868:   if (isinterval) {
869:     RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
870: #if !defined(PETSC_USE_COMPLEX)
871:     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");
872: #endif
873:     if (!ctx->quad && c==d) ctx->quad = EPS_CISS_QUADRULE_CHEBYSHEV;
874:   }
875:   if (!ctx->quad) ctx->quad = EPS_CISS_QUADRULE_TRAPEZOIDAL;

877:   /* create split comm */
878:   if (!ctx->subcomm) { EPSCISSSetUpSubComm(eps,&ctx->num_solve_point); }

880:   EPSAllocateSolution(eps,0);
881:   if (ctx->weight) { PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma); }
882:   PetscMalloc4(ctx->N,&ctx->weight,ctx->N+1,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma);
883:   PetscLogObjectMemory((PetscObject)eps,3*ctx->N*sizeof(PetscScalar)+ctx->L_max*ctx->N*sizeof(PetscReal));

885:   /* allocate basis vectors */
886:   BVDestroy(&ctx->S);
887:   BVDuplicateResize(eps->V,ctx->L_max*ctx->M,&ctx->S);
888:   PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->S);
889:   BVDestroy(&ctx->V);
890:   BVDuplicateResize(eps->V,ctx->L_max,&ctx->V);
891:   PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->V);

893:   STGetMatrix(eps->st,0,&A);
894:   PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);
895:   if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Matrix type shell is not supported in this solver");

897:   if (!ctx->usest_set) ctx->usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
898:   if (ctx->usest && ctx->npart>1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The usest flag is not supported when partitions > 1");

900:   CISSRedundantMat(eps);
901:   if (ctx->pA) {
902:     CISSScatterVec(eps);
903:     BVDestroy(&ctx->pV);
904:     BVCreate(PetscObjectComm((PetscObject)ctx->xsub),&ctx->pV);
905:     BVSetSizesFromVec(ctx->pV,ctx->xsub,eps->n);
906:     BVSetFromOptions(ctx->pV);
907:     BVResize(ctx->pV,ctx->L_max,PETSC_FALSE);
908:     PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->pV);
909:   }

911:   if (ctx->usest) {
912:     PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&issinvert);
913:     if (!issinvert) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"If the usest flag is set, you must select the STSINVERT spectral transformation");
914:   }

916:   BVDestroy(&ctx->Y);
917:   if (ctx->pA) {
918:     BVCreate(PetscObjectComm((PetscObject)ctx->xsub),&ctx->Y);
919:     BVSetSizesFromVec(ctx->Y,ctx->xsub,eps->n);
920:     BVSetFromOptions(ctx->Y);
921:     BVResize(ctx->Y,ctx->num_solve_point*ctx->L_max,PETSC_FALSE);
922:   } else {
923:     BVDuplicateResize(eps->V,ctx->num_solve_point*ctx->L_max,&ctx->Y);
924:   }
925:   PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->Y);

927:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
928:     DSSetType(eps->ds,DSGNHEP);
929:   } else if (eps->isgeneralized) {
930:     if (eps->ishermitian && eps->ispositive) {
931:       DSSetType(eps->ds,DSGHEP);
932:     } else {
933:       DSSetType(eps->ds,DSGNHEP);
934:     }
935:   } else {
936:     if (eps->ishermitian) {
937:       DSSetType(eps->ds,DSHEP);
938:     } else {
939:       DSSetType(eps->ds,DSNHEP);
940:     }
941:   }
942:   DSAllocate(eps->ds,eps->ncv);
943:   EPSSetWorkVecs(eps,2);

945: #if !defined(PETSC_USE_COMPLEX)
946:   if (!eps->ishermitian) { PetscInfo(eps,"Warning: complex eigenvalues are not calculated exactly without --with-scalar-type=complex in PETSc\n"); }
947: #endif
948:   return(0);
949: }

951: PetscErrorCode EPSSolve_CISS(EPS eps)
952: {
954:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
955:   Mat            A,B,X,M,pA,pB;
956:   PetscInt       i,j,ld,nmat,L_add=0,nv=0,L_base=ctx->L,inner,nlocal,*inside;
957:   PetscScalar    *Mu,*H0,*H1=NULL,*rr,*temp;
958:   PetscReal      error,max_error,norm;
959:   PetscBool      *fl1;
960:   Vec            si,w[3];
961:   SlepcSC        sc;
962:   PetscRandom    rand;
963: #if defined(PETSC_USE_COMPLEX)
964:   PetscBool      isellipse;
965: #endif

968:   w[0] = eps->work[0];
969:   w[1] = NULL;
970:   w[2] = eps->work[1];
971:   /* override SC settings */
972:   DSGetSlepcSC(eps->ds,&sc);
973:   sc->comparison    = SlepcCompareLargestMagnitude;
974:   sc->comparisonctx = NULL;
975:   sc->map           = NULL;
976:   sc->mapobj        = NULL;
977:   VecGetLocalSize(w[0],&nlocal);
978:   DSGetLeadingDimension(eps->ds,&ld);
979:   STGetNumMatrices(eps->st,&nmat);
980:   STGetMatrix(eps->st,0,&A);
981:   if (nmat>1) { STGetMatrix(eps->st,1,&B); }
982:   else B = NULL;
983:   SetPathParameter(eps);
984:   CISSVecSetRandom(ctx->V,0,ctx->L);
985:   BVGetRandomContext(ctx->V,&rand);

987:   if (ctx->pA) {
988:     VecScatterVecs(eps,ctx->V,ctx->L);
989:     SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,0,ctx->L,PETSC_TRUE);
990:   } else {
991:     SolveLinearSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_TRUE);
992:   }
993: #if defined(PETSC_USE_COMPLEX)
994:   PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
995:   if (isellipse) {
996:     EstimateNumberEigs(eps,&L_add);
997:   } else {
998:     L_add = 0;
999:   }
1000: #else
1001:   L_add = 0;
1002: #endif
1003:   if (L_add>0) {
1004:     PetscInfo2(eps,"Changing L %D -> %D by Estimate #Eig\n",ctx->L,ctx->L+L_add);
1005:     CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
1006:     if (ctx->pA) {
1007:       VecScatterVecs(eps,ctx->V,ctx->L+L_add);
1008:       SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,ctx->L,ctx->L+L_add,PETSC_FALSE);
1009:     } else {
1010:       SolveLinearSystem(eps,A,B,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
1011:     }
1012:     ctx->L += L_add;
1013:   }
1014:   PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
1015:   for (i=0;i<ctx->refine_blocksize;i++) {
1016:     CalcMu(eps,Mu);
1017:     BlockHankel(eps,Mu,0,H0);
1018:     SVD_H0(eps,H0,&nv);
1019:     if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
1020:     L_add = L_base;
1021:     if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
1022:     PetscInfo2(eps,"Changing L %D -> %D by SVD(H0)\n",ctx->L,ctx->L+L_add);
1023:     CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
1024:     if (ctx->pA) {
1025:       VecScatterVecs(eps,ctx->V,ctx->L+L_add);
1026:       SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,ctx->L,ctx->L+L_add,PETSC_FALSE);
1027:     } else {
1028:       SolveLinearSystem(eps,A,B,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
1029:     }
1030:     ctx->L += L_add;
1031:   }
1032:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1033:     PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1);
1034:   }

1036:   while (eps->reason == EPS_CONVERGED_ITERATING) {
1037:     eps->its++;
1038:     for (inner=0;inner<=ctx->refine_inner;inner++) {
1039:       if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1040:         CalcMu(eps,Mu);
1041:         BlockHankel(eps,Mu,0,H0);
1042:         SVD_H0(eps,H0,&nv);
1043:         break;
1044:       } else {
1045:         ConstructS(eps);
1046:         BVSetActiveColumns(ctx->S,0,ctx->L);
1047:         BVCopy(ctx->S,ctx->V);
1048:         SVD_S(ctx->S,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv);
1049:         if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
1050:           if (ctx->pA) {
1051:             VecScatterVecs(eps,ctx->V,ctx->L);
1052:             SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,0,ctx->L,PETSC_FALSE);
1053:           } else {
1054:             SolveLinearSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_FALSE);
1055:           }
1056:         } else break;
1057:       }
1058:     }
1059:     eps->nconv = 0;
1060:     if (nv == 0) eps->reason = EPS_CONVERGED_TOL;
1061:     else {
1062:       DSSetDimensions(eps->ds,nv,0,0,0);
1063:       DSSetState(eps->ds,DS_STATE_RAW);

1065:       if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1066:         BlockHankel(eps,Mu,0,H0);
1067:         BlockHankel(eps,Mu,1,H1);
1068:         DSGetArray(eps->ds,DS_MAT_A,&temp);
1069:         for (j=0;j<nv;j++) {
1070:           for (i=0;i<nv;i++) {
1071:             temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
1072:           }
1073:         }
1074:         DSRestoreArray(eps->ds,DS_MAT_A,&temp);
1075:         DSGetArray(eps->ds,DS_MAT_B,&temp);
1076:         for (j=0;j<nv;j++) {
1077:           for (i=0;i<nv;i++) {
1078:             temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
1079:           }
1080:         }
1081:         DSRestoreArray(eps->ds,DS_MAT_B,&temp);
1082:       } else {
1083:         BVSetActiveColumns(ctx->S,0,nv);
1084:         DSGetMat(eps->ds,DS_MAT_A,&pA);
1085:         MatZeroEntries(pA);
1086:         BVMatProject(ctx->S,A,ctx->S,pA);
1087:         DSRestoreMat(eps->ds,DS_MAT_A,&pA);
1088:         if (B) {
1089:           DSGetMat(eps->ds,DS_MAT_B,&pB);
1090:           MatZeroEntries(pB);
1091:           BVMatProject(ctx->S,B,ctx->S,pB);
1092:           DSRestoreMat(eps->ds,DS_MAT_B,&pB);
1093:         }
1094:       }

1096:       DSSolve(eps->ds,eps->eigr,eps->eigi);
1097:       DSSynchronize(eps->ds,eps->eigr,eps->eigi);

1099:       PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr);
1100:       rescale_eig(eps,nv);
1101:       isGhost(eps,ld,nv,fl1);
1102:       RGCheckInside(eps->rg,nv,eps->eigr,eps->eigi,inside);
1103:       for (i=0;i<nv;i++) {
1104:         if (fl1[i] && inside[i]>=0) {
1105:           rr[i] = 1.0;
1106:           eps->nconv++;
1107:         } else rr[i] = 0.0;
1108:       }
1109:       DSSort(eps->ds,eps->eigr,eps->eigi,rr,NULL,&eps->nconv);
1110:       DSSynchronize(eps->ds,eps->eigr,eps->eigi);
1111:       rescale_eig(eps,nv);
1112:       PetscFree3(fl1,inside,rr);
1113:       BVSetActiveColumns(eps->V,0,nv);
1114:       if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1115:         ConstructS(eps);
1116:         BVSetActiveColumns(ctx->S,0,ctx->L);
1117:         BVCopy(ctx->S,ctx->V);
1118:         BVSetActiveColumns(ctx->S,0,nv);
1119:       }
1120:       BVCopy(ctx->S,eps->V);

1122:       DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
1123:       DSGetMat(eps->ds,DS_MAT_X,&X);
1124:       BVMultInPlace(ctx->S,X,0,eps->nconv);
1125:       if (eps->ishermitian) {
1126:         BVMultInPlace(eps->V,X,0,eps->nconv);
1127:       }
1128:       MatDestroy(&X);
1129:       max_error = 0.0;
1130:       for (i=0;i<eps->nconv;i++) {
1131:         BVGetColumn(ctx->S,i,&si);
1132:         EPSComputeResidualNorm_Private(eps,PETSC_FALSE,eps->eigr[i],eps->eigi[i],si,NULL,w,&error);
1133:         if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {  /* vector is not normalized */
1134:           VecNorm(si,NORM_2,&norm);
1135:           error /= norm;
1136:         }
1137:         (*eps->converged)(eps,eps->eigr[i],eps->eigi[i],error,&error,eps->convergedctx);
1138:         BVRestoreColumn(ctx->S,i,&si);
1139:         max_error = PetscMax(max_error,error);
1140:       }

1142:       if (max_error <= eps->tol) eps->reason = EPS_CONVERGED_TOL;
1143:       else if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
1144:       else {
1145:         if (eps->nconv > ctx->L) {
1146:           MatCreateSeqDense(PETSC_COMM_SELF,eps->nconv,ctx->L,NULL,&M);
1147:           MatDenseGetArray(M,&temp);
1148:           for (i=0;i<ctx->L*eps->nconv;i++) {
1149:             PetscRandomGetValue(rand,&temp[i]);
1150:             temp[i] = PetscRealPart(temp[i]);
1151:           }
1152:           MatDenseRestoreArray(M,&temp);
1153:           BVSetActiveColumns(ctx->S,0,eps->nconv);
1154:           BVMultInPlace(ctx->S,M,0,ctx->L);
1155:           MatDestroy(&M);
1156:           BVSetActiveColumns(ctx->S,0,ctx->L);
1157:           BVCopy(ctx->S,ctx->V);
1158:         }
1159:         if (ctx->pA) {
1160:           VecScatterVecs(eps,ctx->V,ctx->L);
1161:           SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,0,ctx->L,PETSC_FALSE);
1162:         } else {
1163:           SolveLinearSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_FALSE);
1164:         }
1165:       }
1166:     }
1167:   }
1168:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1169:     PetscFree(H1);
1170:   }
1171:   PetscFree2(Mu,H0);
1172:   return(0);
1173: }

1175: PetscErrorCode EPSComputeVectors_CISS(EPS eps)
1176: {
1178:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1179:   PetscReal      norm;
1180:   PetscInt       i;
1181:   Mat            B=NULL;

1184:   EPSComputeVectors_Schur(eps);
1185:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) { /* normalize */
1186:     if (eps->isgeneralized && eps->ishermitian && eps->ispositive) {
1187:       STGetMatrix(eps->st,1,&B);
1188:       BVSetMatrix(eps->V,B,PETSC_FALSE);
1189:     }
1190:     for (i=0;i<eps->nconv;i++) {
1191:       BVNormColumn(eps->V,i,NORM_2,&norm);
1192:       BVScaleColumn(eps->V,i,1.0/norm);
1193:     }
1194:     if (B) { BVSetMatrix(eps->V,NULL,PETSC_FALSE); }
1195:   }
1196:   return(0);
1197: }

1199: static PetscErrorCode EPSCISSSetSizes_CISS(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
1200: {
1202:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1203:   PetscInt       oN,onpart;

1206:   oN = ctx->N;
1207:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
1208:     if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
1209:   } else {
1210:     if (ip<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
1211:     if (ip%2) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
1212:     if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
1213:   }
1214:   if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
1215:     ctx->L = 16;
1216:   } else {
1217:     if (bs<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
1218:     ctx->L = bs;
1219:   }
1220:   if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
1221:     ctx->M = ctx->N/4;
1222:   } else {
1223:     if (ms<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
1224:     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");
1225:     ctx->M = ms;
1226:   }
1227:   onpart = ctx->npart;
1228:   if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
1229:     ctx->npart = 1;
1230:   } else {
1231:     if (npart<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The npart argument must be > 0");
1232:     ctx->npart = npart;
1233:   }
1234:   if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
1235:     ctx->L_max = 64;
1236:   } else {
1237:     if (bsmax<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
1238:     ctx->L_max = PetscMax(bsmax,ctx->L);
1239:   }
1240:   if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) { EPSCISSResetSubcomm(eps); }
1241:   ctx->isreal = realmats;
1242:   eps->state = EPS_STATE_INITIAL;
1243:   return(0);
1244: }

1246: /*@
1247:    EPSCISSSetSizes - Sets the values of various size parameters in the CISS solver.

1249:    Logically Collective on eps

1251:    Input Parameters:
1252: +  eps   - the eigenproblem solver context
1253: .  ip    - number of integration points
1254: .  bs    - block size
1255: .  ms    - moment size
1256: .  npart - number of partitions when splitting the communicator
1257: .  bsmax - max block size
1258: -  realmats - A and B are real

1260:    Options Database Keys:
1261: +  -eps_ciss_integration_points - Sets the number of integration points
1262: .  -eps_ciss_blocksize - Sets the block size
1263: .  -eps_ciss_moments - Sets the moment size
1264: .  -eps_ciss_partitions - Sets the number of partitions
1265: .  -eps_ciss_maxblocksize - Sets the maximum block size
1266: -  -eps_ciss_realmats - A and B are real

1268:    Note:
1269:    The default number of partitions is 1. This means the internal KSP object is shared
1270:    among all processes of the EPS communicator. Otherwise, the communicator is split
1271:    into npart communicators, so that npart KSP solves proceed simultaneously.

1273:    Level: advanced

1275: .seealso: EPSCISSGetSizes()
1276: @*/
1277: PetscErrorCode EPSCISSSetSizes(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
1278: {

1289:   PetscTryMethod(eps,"EPSCISSSetSizes_C",(EPS,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(eps,ip,bs,ms,npart,bsmax,realmats));
1290:   return(0);
1291: }

1293: static PetscErrorCode EPSCISSGetSizes_CISS(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
1294: {
1295:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1298:   if (ip) *ip = ctx->N;
1299:   if (bs) *bs = ctx->L;
1300:   if (ms) *ms = ctx->M;
1301:   if (npart) *npart = ctx->npart;
1302:   if (bsmax) *bsmax = ctx->L_max;
1303:   if (realmats) *realmats = ctx->isreal;
1304:   return(0);
1305: }

1307: /*@
1308:    EPSCISSGetSizes - Gets the values of various size parameters in the CISS solver.

1310:    Not Collective

1312:    Input Parameter:
1313: .  eps - the eigenproblem solver context

1315:    Output Parameters:
1316: +  ip    - number of integration points
1317: .  bs    - block size
1318: .  ms    - moment size
1319: .  npart - number of partitions when splitting the communicator
1320: .  bsmax - max block size
1321: -  realmats - A and B are real

1323:    Level: advanced

1325: .seealso: EPSCISSSetSizes()
1326: @*/
1327: PetscErrorCode EPSCISSGetSizes(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
1328: {

1333:   PetscUseMethod(eps,"EPSCISSGetSizes_C",(EPS,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(eps,ip,bs,ms,npart,bsmax,realmats));
1334:   return(0);
1335: }

1337: static PetscErrorCode EPSCISSSetThreshold_CISS(EPS eps,PetscReal delta,PetscReal spur)
1338: {
1339:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1342:   if (delta == PETSC_DEFAULT) {
1343:     ctx->delta = 1e-12;
1344:   } else {
1345:     if (delta<=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
1346:     ctx->delta = delta;
1347:   }
1348:   if (spur == PETSC_DEFAULT) {
1349:     ctx->spurious_threshold = 1e-4;
1350:   } else {
1351:     if (spur<=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
1352:     ctx->spurious_threshold = spur;
1353:   }
1354:   return(0);
1355: }

1357: /*@
1358:    EPSCISSSetThreshold - Sets the values of various threshold parameters in
1359:    the CISS solver.

1361:    Logically Collective on eps

1363:    Input Parameters:
1364: +  eps   - the eigenproblem solver context
1365: .  delta - threshold for numerical rank
1366: -  spur  - spurious threshold (to discard spurious eigenpairs)

1368:    Options Database Keys:
1369: +  -eps_ciss_delta - Sets the delta
1370: -  -eps_ciss_spurious_threshold - Sets the spurious threshold

1372:    Level: advanced

1374: .seealso: EPSCISSGetThreshold()
1375: @*/
1376: PetscErrorCode EPSCISSSetThreshold(EPS eps,PetscReal delta,PetscReal spur)
1377: {

1384:   PetscTryMethod(eps,"EPSCISSSetThreshold_C",(EPS,PetscReal,PetscReal),(eps,delta,spur));
1385:   return(0);
1386: }

1388: static PetscErrorCode EPSCISSGetThreshold_CISS(EPS eps,PetscReal *delta,PetscReal *spur)
1389: {
1390:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1393:   if (delta) *delta = ctx->delta;
1394:   if (spur)  *spur = ctx->spurious_threshold;
1395:   return(0);
1396: }

1398: /*@
1399:    EPSCISSGetThreshold - Gets the values of various threshold parameters
1400:    in the CISS solver.

1402:    Not Collective

1404:    Input Parameter:
1405: .  eps - the eigenproblem solver context

1407:    Output Parameters:
1408: +  delta - threshold for numerical rank
1409: -  spur  - spurious threshold (to discard spurious eigenpairs)

1411:    Level: advanced

1413: .seealso: EPSCISSSetThreshold()
1414: @*/
1415: PetscErrorCode EPSCISSGetThreshold(EPS eps,PetscReal *delta,PetscReal *spur)
1416: {

1421:   PetscUseMethod(eps,"EPSCISSGetThreshold_C",(EPS,PetscReal*,PetscReal*),(eps,delta,spur));
1422:   return(0);
1423: }

1425: static PetscErrorCode EPSCISSSetRefinement_CISS(EPS eps,PetscInt inner,PetscInt blsize)
1426: {
1427:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1430:   if (inner == PETSC_DEFAULT) {
1431:     ctx->refine_inner = 0;
1432:   } else {
1433:     if (inner<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
1434:     ctx->refine_inner = inner;
1435:   }
1436:   if (blsize == PETSC_DEFAULT) {
1437:     ctx->refine_blocksize = 0;
1438:   } else {
1439:     if (blsize<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
1440:     ctx->refine_blocksize = blsize;
1441:   }
1442:   return(0);
1443: }

1445: /*@
1446:    EPSCISSSetRefinement - Sets the values of various refinement parameters
1447:    in the CISS solver.

1449:    Logically Collective on eps

1451:    Input Parameters:
1452: +  eps    - the eigenproblem solver context
1453: .  inner  - number of iterative refinement iterations (inner loop)
1454: -  blsize - number of iterative refinement iterations (blocksize loop)

1456:    Options Database Keys:
1457: +  -eps_ciss_refine_inner - Sets number of inner iterations
1458: -  -eps_ciss_refine_blocksize - Sets number of blocksize iterations

1460:    Level: advanced

1462: .seealso: EPSCISSGetRefinement()
1463: @*/
1464: PetscErrorCode EPSCISSSetRefinement(EPS eps,PetscInt inner,PetscInt blsize)
1465: {

1472:   PetscTryMethod(eps,"EPSCISSSetRefinement_C",(EPS,PetscInt,PetscInt),(eps,inner,blsize));
1473:   return(0);
1474: }

1476: static PetscErrorCode EPSCISSGetRefinement_CISS(EPS eps,PetscInt *inner,PetscInt *blsize)
1477: {
1478:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1481:   if (inner)  *inner = ctx->refine_inner;
1482:   if (blsize) *blsize = ctx->refine_blocksize;
1483:   return(0);
1484: }

1486: /*@
1487:    EPSCISSGetRefinement - Gets the values of various refinement parameters
1488:    in the CISS solver.

1490:    Not Collective

1492:    Input Parameter:
1493: .  eps - the eigenproblem solver context

1495:    Output Parameters:
1496: +  inner  - number of iterative refinement iterations (inner loop)
1497: -  blsize - number of iterative refinement iterations (blocksize loop)

1499:    Level: advanced

1501: .seealso: EPSCISSSetRefinement()
1502: @*/
1503: PetscErrorCode EPSCISSGetRefinement(EPS eps, PetscInt *inner, PetscInt *blsize)
1504: {

1509:   PetscUseMethod(eps,"EPSCISSGetRefinement_C",(EPS,PetscInt*,PetscInt*),(eps,inner,blsize));
1510:   return(0);
1511: }

1513: static PetscErrorCode EPSCISSSetUseST_CISS(EPS eps,PetscBool usest)
1514: {
1515:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1518:   ctx->usest     = usest;
1519:   ctx->usest_set = PETSC_TRUE;
1520:   eps->state     = EPS_STATE_INITIAL;
1521:   return(0);
1522: }

1524: /*@
1525:    EPSCISSSetUseST - Sets a flag indicating that the CISS solver will
1526:    use the ST object for the linear solves.

1528:    Logically Collective on eps

1530:    Input Parameters:
1531: +  eps    - the eigenproblem solver context
1532: -  usest  - boolean flag to use the ST object or not

1534:    Options Database Keys:
1535: .  -eps_ciss_usest <bool> - whether the ST object will be used or not

1537:    Level: advanced

1539: .seealso: EPSCISSGetUseST()
1540: @*/
1541: PetscErrorCode EPSCISSSetUseST(EPS eps,PetscBool usest)
1542: {

1548:   PetscTryMethod(eps,"EPSCISSSetUseST_C",(EPS,PetscBool),(eps,usest));
1549:   return(0);
1550: }

1552: static PetscErrorCode EPSCISSGetUseST_CISS(EPS eps,PetscBool *usest)
1553: {
1554:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1557:   *usest = ctx->usest;
1558:   return(0);
1559: }

1561: /*@
1562:    EPSCISSGetUseST - Gets the flag for using the ST object
1563:    in the CISS solver.

1565:    Not Collective

1567:    Input Parameter:
1568: .  eps - the eigenproblem solver context

1570:    Output Parameters:
1571: .  usest - boolean flag indicating if the ST object is being used

1573:    Level: advanced

1575: .seealso: EPSCISSSetUseST()
1576: @*/
1577: PetscErrorCode EPSCISSGetUseST(EPS eps,PetscBool *usest)
1578: {

1584:   PetscUseMethod(eps,"EPSCISSGetUseST_C",(EPS,PetscBool*),(eps,usest));
1585:   return(0);
1586: }

1588: static PetscErrorCode EPSCISSSetQuadRule_CISS(EPS eps,EPSCISSQuadRule quad)
1589: {
1590:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1593:   ctx->quad = quad;
1594:   return(0);
1595: }

1597: /*@
1598:    EPSCISSSetQuadRule - Sets the quadrature rule used in the CISS solver.

1600:    Logically Collective on eps

1602:    Input Parameters:
1603: +  eps  - the eigenproblem solver context
1604: -  quad - the quadrature rule

1606:    Options Database Key:
1607: .  -eps_ciss_quadrule - Sets the quadrature rule (either 'trapezoidal' or
1608:                            'chebyshev')

1610:    Notes:
1611:    By default, the trapezoidal rule is used (EPS_CISS_QUADRULE_TRAPEZOIDAL).

1613:    If the 'chebyshev' option is specified (EPS_CISS_QUADRULE_CHEBYSHEV), then
1614:    Chebyshev points are used as quadrature points.

1616:    Level: advanced

1618: .seealso: EPSCISSGetQuadRule(), EPSCISSQuadRule
1619: @*/
1620: PetscErrorCode EPSCISSSetQuadRule(EPS eps,EPSCISSQuadRule quad)
1621: {

1627:   PetscTryMethod(eps,"EPSCISSSetQuadRule_C",(EPS,EPSCISSQuadRule),(eps,quad));
1628:   return(0);
1629: }

1631: static PetscErrorCode EPSCISSGetQuadRule_CISS(EPS eps,EPSCISSQuadRule *quad)
1632: {
1633:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1636:   *quad = ctx->quad;
1637:   return(0);
1638: }

1640: /*@
1641:    EPSCISSGetQuadRule - Gets the quadrature rule used in the CISS solver.

1643:    Not Collective

1645:    Input Parameter:
1646: .  eps - the eigenproblem solver context

1648:    Output Parameters:
1649: .  quad - quadrature rule

1651:    Level: advanced

1653: .seealso: EPSCISSSetQuadRule() EPSCISSQuadRule
1654: @*/
1655: PetscErrorCode EPSCISSGetQuadRule(EPS eps,EPSCISSQuadRule *quad)
1656: {

1662:   PetscUseMethod(eps,"EPSCISSGetQuadRule_C",(EPS,EPSCISSQuadRule*),(eps,quad));
1663:   return(0);
1664: }

1666: static PetscErrorCode EPSCISSSetExtraction_CISS(EPS eps,EPSCISSExtraction extraction)
1667: {
1668:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1671:   ctx->extraction = extraction;
1672:   return(0);
1673: }

1675: /*@
1676:    EPSCISSSetExtraction - Sets the extraction technique used in the CISS solver.

1678:    Logically Collective on eps

1680:    Input Parameters:
1681: +  eps        - the eigenproblem solver context
1682: -  extraction - the extraction technique

1684:    Options Database Key:
1685: .  -eps_ciss_extraction - Sets the extraction technique (either 'ritz' or
1686:                            'hankel')

1688:    Notes:
1689:    By default, the Rayleigh-Ritz extraction is used (EPS_CISS_EXTRACTION_RITZ).

1691:    If the 'hankel' option is specified (EPS_CISS_EXTRACTION_HANKEL), then
1692:    the Block Hankel method is used for extracting eigenpairs.

1694:    Level: advanced

1696: .seealso: EPSCISSGetExtraction(), EPSCISSExtraction
1697: @*/
1698: PetscErrorCode EPSCISSSetExtraction(EPS eps,EPSCISSExtraction extraction)
1699: {

1705:   PetscTryMethod(eps,"EPSCISSSetExtraction_C",(EPS,EPSCISSExtraction),(eps,extraction));
1706:   return(0);
1707: }

1709: static PetscErrorCode EPSCISSGetExtraction_CISS(EPS eps,EPSCISSExtraction *extraction)
1710: {
1711:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1714:   *extraction = ctx->extraction;
1715:   return(0);
1716: }

1718: /*@
1719:    EPSCISSGetExtraction - Gets the extraction technique used in the CISS solver.

1721:    Not Collective

1723:    Input Parameter:
1724: .  eps - the eigenproblem solver context

1726:    Output Parameters:
1727: +  extraction - extraction technique

1729:    Level: advanced

1731: .seealso: EPSCISSSetExtraction() EPSCISSExtraction
1732: @*/
1733: PetscErrorCode EPSCISSGetExtraction(EPS eps,EPSCISSExtraction *extraction)
1734: {

1740:   PetscUseMethod(eps,"EPSCISSGetExtraction_C",(EPS,EPSCISSExtraction*),(eps,extraction));
1741:   return(0);
1742: }

1744: static PetscErrorCode EPSCISSGetKSPs_CISS(EPS eps,PetscInt *nsolve,KSP **ksp)
1745: {
1747:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1748:   PetscInt       i;
1749:   PC             pc;

1752:   if (!ctx->ksp) {
1753:     if (!ctx->subcomm) {  /* initialize subcomm first */
1754:       EPSCISSSetUseConj(eps,&ctx->useconj);
1755:       EPSCISSSetUpSubComm(eps,&ctx->num_solve_point);
1756:     }
1757:     PetscMalloc1(ctx->num_solve_point,&ctx->ksp);
1758:     for (i=0;i<ctx->num_solve_point;i++) {
1759:       KSPCreate(PetscSubcommChild(ctx->subcomm),&ctx->ksp[i]);
1760:       PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)eps,1);
1761:       KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)eps)->prefix);
1762:       KSPAppendOptionsPrefix(ctx->ksp[i],"eps_ciss_");
1763:       PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->ksp[i]);
1764:       PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)eps)->options);
1765:       KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1766:       KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1767:       KSPGetPC(ctx->ksp[i],&pc);
1768:       KSPSetType(ctx->ksp[i],KSPPREONLY);
1769:       PCSetType(pc,PCLU);
1770:     }
1771:   }
1772:   if (nsolve) *nsolve = ctx->num_solve_point;
1773:   if (ksp)    *ksp    = ctx->ksp;
1774:   return(0);
1775: }

1777: /*@C
1778:    EPSCISSGetKSPs - Retrieve the array of linear solver objects associated with
1779:    the CISS solver.

1781:    Not Collective

1783:    Input Parameter:
1784: .  eps - the eigenproblem solver solver

1786:    Output Parameters:
1787: +  nsolve - number of solver objects
1788: -  ksp - array of linear solver object

1790:    Notes:
1791:    The number of KSP solvers is equal to the number of integration points divided by
1792:    the number of partitions. This value is halved in the case of real matrices with
1793:    a region centered at the real axis.

1795:    Level: advanced

1797: .seealso: EPSCISSSetSizes()
1798: @*/
1799: PetscErrorCode EPSCISSGetKSPs(EPS eps,PetscInt *nsolve,KSP **ksp)
1800: {

1805:   PetscUseMethod(eps,"EPSCISSGetKSPs_C",(EPS,PetscInt*,KSP**),(eps,nsolve,ksp));
1806:   return(0);
1807: }

1809: PetscErrorCode EPSReset_CISS(EPS eps)
1810: {
1812:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1813:   PetscInt       i;

1816:   BVDestroy(&ctx->S);
1817:   BVDestroy(&ctx->V);
1818:   BVDestroy(&ctx->Y);
1819:   if (!ctx->usest) {
1820:     for (i=0;i<ctx->num_solve_point;i++) {
1821:       KSPReset(ctx->ksp[i]);
1822:     }
1823:   }
1824:   VecScatterDestroy(&ctx->scatterin);
1825:   VecDestroy(&ctx->xsub);
1826:   VecDestroy(&ctx->xdup);
1827:   if (ctx->pA) {
1828:     MatDestroy(&ctx->pA);
1829:     MatDestroy(&ctx->pB);
1830:     BVDestroy(&ctx->pV);
1831:   }
1832:   return(0);
1833: }

1835: PetscErrorCode EPSSetFromOptions_CISS(PetscOptionItems *PetscOptionsObject,EPS eps)
1836: {
1837:   PetscErrorCode    ierr;
1838:   PetscReal         r3,r4;
1839:   PetscInt          i,i1,i2,i3,i4,i5,i6,i7;
1840:   PetscBool         b1,b2,flg;
1841:   EPS_CISS          *ctx = (EPS_CISS*)eps->data;
1842:   EPSCISSQuadRule   quad;
1843:   EPSCISSExtraction extraction;

1846:   PetscOptionsHead(PetscOptionsObject,"EPS CISS Options");

1848:     EPSCISSGetSizes(eps,&i1,&i2,&i3,&i4,&i5,&b1);
1849:     PetscOptionsInt("-eps_ciss_integration_points","Number of integration points","EPSCISSSetSizes",i1,&i1,NULL);
1850:     PetscOptionsInt("-eps_ciss_blocksize","Block size","EPSCISSSetSizes",i2,&i2,NULL);
1851:     PetscOptionsInt("-eps_ciss_moments","Moment size","EPSCISSSetSizes",i3,&i3,NULL);
1852:     PetscOptionsInt("-eps_ciss_partitions","Number of partitions","EPSCISSSetSizes",i4,&i4,NULL);
1853:     PetscOptionsInt("-eps_ciss_maxblocksize","Maximum block size","EPSCISSSetSizes",i5,&i5,NULL);
1854:     PetscOptionsBool("-eps_ciss_realmats","True if A and B are real","EPSCISSSetSizes",b1,&b1,NULL);
1855:     EPSCISSSetSizes(eps,i1,i2,i3,i4,i5,b1);

1857:     EPSCISSGetThreshold(eps,&r3,&r4);
1858:     PetscOptionsReal("-eps_ciss_delta","Threshold for numerical rank","EPSCISSSetThreshold",r3,&r3,NULL);
1859:     PetscOptionsReal("-eps_ciss_spurious_threshold","Threshold for the spurious eigenpairs","EPSCISSSetThreshold",r4,&r4,NULL);
1860:     EPSCISSSetThreshold(eps,r3,r4);

1862:     EPSCISSGetRefinement(eps,&i6,&i7);
1863:     PetscOptionsInt("-eps_ciss_refine_inner","Number of inner iterative refinement iterations","EPSCISSSetRefinement",i6,&i6,NULL);
1864:     PetscOptionsInt("-eps_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","EPSCISSSetRefinement",i7,&i7,NULL);
1865:     EPSCISSSetRefinement(eps,i6,i7);

1867:     EPSCISSGetUseST(eps,&b2);
1868:     PetscOptionsBool("-eps_ciss_usest","Use ST for linear solves","EPSCISSSetUseST",b2,&b2,&flg);
1869:     if (flg) { EPSCISSSetUseST(eps,b2); }

1871:     PetscOptionsEnum("-eps_ciss_quadrule","Quadrature rule","EPSCISSSetQuadRule",EPSCISSQuadRules,(PetscEnum)ctx->quad,(PetscEnum*)&quad,&flg);
1872:     if (flg) { EPSCISSSetQuadRule(eps,quad); }

1874:     PetscOptionsEnum("-eps_ciss_extraction","Extraction technique","EPSCISSSetExtraction",EPSCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg);
1875:     if (flg) { EPSCISSSetExtraction(eps,extraction); }

1877:   PetscOptionsTail();

1879:   if (!eps->rg) { EPSGetRG(eps,&eps->rg); }
1880:   RGSetFromOptions(eps->rg); /* this is necessary here to set useconj */
1881:   if (!ctx->ksp) { EPSCISSGetKSPs(eps,&ctx->num_solve_point,&ctx->ksp); }
1882:   for (i=0;i<ctx->num_solve_point;i++) {
1883:     KSPSetFromOptions(ctx->ksp[i]);
1884:   }
1885:   PetscSubcommSetFromOptions(ctx->subcomm);
1886:   return(0);
1887: }

1889: PetscErrorCode EPSDestroy_CISS(EPS eps)
1890: {
1892:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1895:   EPSCISSResetSubcomm(eps);
1896:   PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1897:   PetscFree(eps->data);
1898:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",NULL);
1899:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",NULL);
1900:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",NULL);
1901:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",NULL);
1902:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",NULL);
1903:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",NULL);
1904:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",NULL);
1905:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",NULL);
1906:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",NULL);
1907:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",NULL);
1908:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",NULL);
1909:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",NULL);
1910:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",NULL);
1911:   return(0);
1912: }

1914: PetscErrorCode EPSView_CISS(EPS eps,PetscViewer viewer)
1915: {
1917:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1918:   PetscBool      isascii;

1921:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1922:   if (isascii) {
1923:     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);
1924:     if (ctx->isreal) {
1925:       PetscViewerASCIIPrintf(viewer,"  exploiting symmetry of integration points\n");
1926:     }
1927:     PetscViewerASCIIPrintf(viewer,"  threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1928:     PetscViewerASCIIPrintf(viewer,"  iterative refinement { inner: %D, blocksize: %D }\n",ctx->refine_inner, ctx->refine_blocksize);
1929:     PetscViewerASCIIPrintf(viewer,"  extraction: %s\n",EPSCISSExtractions[ctx->extraction]);
1930:     PetscViewerASCIIPrintf(viewer,"  quadrature rule: %s\n",EPSCISSQuadRules[ctx->quad]);
1931:     if (ctx->usest) {
1932:       PetscViewerASCIIPrintf(viewer,"  using ST for linear solves\n");
1933:     } else {
1934:       if (!ctx->ksp) { EPSCISSGetKSPs(eps,&ctx->num_solve_point,&ctx->ksp); }
1935:       PetscViewerASCIIPushTab(viewer);
1936:       KSPView(ctx->ksp[0],viewer);
1937:       PetscViewerASCIIPopTab(viewer);
1938:     }
1939:   }
1940:   return(0);
1941: }

1943: PetscErrorCode EPSSetDefaultST_CISS(EPS eps)
1944: {
1946:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1947:   PetscBool      usest = ctx->usest;

1950:   if (!((PetscObject)eps->st)->type_name) {
1951:     if (!ctx->usest_set) usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
1952:     if (usest) {
1953:       STSetType(eps->st,STSINVERT);
1954:     } else {
1955:       /* we are not going to use ST, so avoid factorizing the matrix */
1956:       STSetType(eps->st,STSHIFT);
1957:     }
1958:   }
1959:   return(0);
1960: }

1962: SLEPC_EXTERN PetscErrorCode EPSCreate_CISS(EPS eps)
1963: {
1965:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1968:   PetscNewLog(eps,&ctx);
1969:   eps->data = ctx;

1971:   eps->useds = PETSC_TRUE;
1972:   eps->categ = EPS_CATEGORY_CONTOUR;

1974:   eps->ops->solve          = EPSSolve_CISS;
1975:   eps->ops->setup          = EPSSetUp_CISS;
1976:   eps->ops->setfromoptions = EPSSetFromOptions_CISS;
1977:   eps->ops->destroy        = EPSDestroy_CISS;
1978:   eps->ops->reset          = EPSReset_CISS;
1979:   eps->ops->view           = EPSView_CISS;
1980:   eps->ops->computevectors = EPSComputeVectors_CISS;
1981:   eps->ops->setdefaultst   = EPSSetDefaultST_CISS;

1983:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",EPSCISSSetSizes_CISS);
1984:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",EPSCISSGetSizes_CISS);
1985:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",EPSCISSSetThreshold_CISS);
1986:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",EPSCISSGetThreshold_CISS);
1987:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",EPSCISSSetRefinement_CISS);
1988:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",EPSCISSGetRefinement_CISS);
1989:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",EPSCISSSetUseST_CISS);
1990:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",EPSCISSGetUseST_CISS);
1991:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",EPSCISSSetQuadRule_CISS);
1992:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",EPSCISSGetQuadRule_CISS);
1993:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",EPSCISSSetExtraction_CISS);
1994:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",EPSCISSGetExtraction_CISS);
1995:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",EPSCISSGetKSPs_CISS);
1996:   /* set default values of parameters */
1997:   ctx->N                  = 32;
1998:   ctx->L                  = 16;
1999:   ctx->M                  = ctx->N/4;
2000:   ctx->delta              = 1e-12;
2001:   ctx->L_max              = 64;
2002:   ctx->spurious_threshold = 1e-4;
2003:   ctx->usest              = PETSC_TRUE;
2004:   ctx->usest_set          = PETSC_FALSE;
2005:   ctx->isreal             = PETSC_FALSE;
2006:   ctx->refine_inner       = 0;
2007:   ctx->refine_blocksize   = 0;
2008:   ctx->npart              = 1;
2009:   ctx->quad               = (EPSCISSQuadRule)0;
2010:   ctx->extraction         = EPS_CISS_EXTRACTION_RITZ;
2011:   return(0);
2012: }