Actual source code: linear.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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:    Explicit linearization for polynomial eigenproblems
 12: */

 14: #include <slepc/private/pepimpl.h>
 15: #include "linear.h"

 17: static PetscErrorCode MatMult_Linear_Shift(Mat M,Vec x,Vec y)
 18: {
 19:   PetscErrorCode    ierr;
 20:   PEP_LINEAR        *ctx;
 21:   PEP               pep;
 22:   const PetscScalar *px;
 23:   PetscScalar       *py,a,sigma=0.0;
 24:   PetscInt          nmat,deg,i,m;
 25:   Vec               x1,x2,x3,y1,aux;
 26:   PetscReal         *ca,*cb,*cg;
 27:   PetscBool         flg;

 30:   MatShellGetContext(M,&ctx);
 31:   pep = ctx->pep;
 32:   STGetTransform(pep->st,&flg);
 33:   if (!flg) {
 34:     STGetShift(pep->st,&sigma);
 35:   }
 36:   nmat = pep->nmat;
 37:   deg = nmat-1;
 38:   m = pep->nloc;
 39:   ca = pep->pbc;
 40:   cb = pep->pbc+nmat;
 41:   cg = pep->pbc+2*nmat;
 42:   x1=ctx->w[0];x2=ctx->w[1];x3=ctx->w[2];y1=ctx->w[3];aux=ctx->w[4];

 44:   VecSet(y,0.0);
 45:   VecGetArrayRead(x,&px);
 46:   VecGetArray(y,&py);
 47:   a = 1.0;

 49:   /* first block */
 50:   VecPlaceArray(x2,px);
 51:   VecPlaceArray(x3,px+m);
 52:   VecPlaceArray(y1,py);
 53:   VecAXPY(y1,cb[0]-sigma,x2);
 54:   VecAXPY(y1,ca[0],x3);
 55:   VecResetArray(x2);
 56:   VecResetArray(x3);
 57:   VecResetArray(y1);

 59:   /* inner blocks */
 60:   for (i=1;i<deg-1;i++) {
 61:     VecPlaceArray(x1,px+(i-1)*m);
 62:     VecPlaceArray(x2,px+i*m);
 63:     VecPlaceArray(x3,px+(i+1)*m);
 64:     VecPlaceArray(y1,py+i*m);
 65:     VecAXPY(y1,cg[i],x1);
 66:     VecAXPY(y1,cb[i]-sigma,x2);
 67:     VecAXPY(y1,ca[i],x3);
 68:     VecResetArray(x1);
 69:     VecResetArray(x2);
 70:     VecResetArray(x3);
 71:     VecResetArray(y1);
 72:   }

 74:   /* last block */
 75:   VecPlaceArray(y1,py+(deg-1)*m);
 76:   for (i=0;i<deg;i++) {
 77:     VecPlaceArray(x1,px+i*m);
 78:     STMatMult(pep->st,i,x1,aux);
 79:     VecAXPY(y1,a,aux);
 80:     VecResetArray(x1);
 81:     a *= pep->sfactor;
 82:   }
 83:   VecCopy(y1,aux);
 84:   STMatSolve(pep->st,aux,y1);
 85:   VecScale(y1,-ca[deg-1]/a);
 86:   VecPlaceArray(x1,px+(deg-2)*m);
 87:   VecPlaceArray(x2,px+(deg-1)*m);
 88:   VecAXPY(y1,cg[deg-1],x1);
 89:   VecAXPY(y1,cb[deg-1]-sigma,x2);
 90:   VecResetArray(x1);
 91:   VecResetArray(x2);
 92:   VecResetArray(y1);

 94:   VecRestoreArrayRead(x,&px);
 95:   VecRestoreArray(y,&py);
 96:   return(0);
 97: }

 99: static PetscErrorCode MatMult_Linear_Sinvert(Mat M,Vec x,Vec y)
100: {
101:   PetscErrorCode    ierr;
102:   PEP_LINEAR        *ctx;
103:   PEP               pep;
104:   const PetscScalar *px;
105:   PetscScalar       *py,a,sigma,t=1.0,tp=0.0,tt;
106:   PetscInt          nmat,deg,i,m;
107:   Vec               x1,y1,y2,y3,aux,aux2;
108:   PetscReal         *ca,*cb,*cg;

111:   MatShellGetContext(M,&ctx);
112:   pep = ctx->pep;
113:   nmat = pep->nmat;
114:   deg = nmat-1;
115:   m = pep->nloc;
116:   ca = pep->pbc;
117:   cb = pep->pbc+nmat;
118:   cg = pep->pbc+2*nmat;
119:   x1=ctx->w[0];y1=ctx->w[1];y2=ctx->w[2];y3=ctx->w[3];aux=ctx->w[4];aux2=ctx->w[5];
120:   EPSGetTarget(ctx->eps,&sigma);
121:   VecSet(y,0.0);
122:   VecGetArrayRead(x,&px);
123:   VecGetArray(y,&py);
124:   a = pep->sfactor;

126:   /* first block */
127:   VecPlaceArray(x1,px);
128:   VecPlaceArray(y1,py+m);
129:   VecCopy(x1,y1);
130:   VecScale(y1,1.0/ca[0]);
131:   VecResetArray(x1);
132:   VecResetArray(y1);

134:   /* second block */
135:   if (deg>2) {
136:     VecPlaceArray(x1,px+m);
137:     VecPlaceArray(y1,py+m);
138:     VecPlaceArray(y2,py+2*m);
139:     VecCopy(x1,y2);
140:     VecAXPY(y2,sigma-cb[1],y1);
141:     VecScale(y2,1.0/ca[1]);
142:     VecResetArray(x1);
143:     VecResetArray(y1);
144:     VecResetArray(y2);
145:   }

147:   /* inner blocks */
148:   for (i=2;i<deg-1;i++) {
149:     VecPlaceArray(x1,px+i*m);
150:     VecPlaceArray(y1,py+(i-1)*m);
151:     VecPlaceArray(y2,py+i*m);
152:     VecPlaceArray(y3,py+(i+1)*m);
153:     VecCopy(x1,y3);
154:     VecAXPY(y3,sigma-cb[i],y2);
155:     VecAXPY(y3,-cg[i],y1);
156:     VecScale(y3,1.0/ca[i]);
157:     VecResetArray(x1);
158:     VecResetArray(y1);
159:     VecResetArray(y2);
160:     VecResetArray(y3);
161:   }

163:   /* last block */
164:   VecPlaceArray(y1,py);
165:   for (i=0;i<deg-2;i++) {
166:     VecPlaceArray(y2,py+(i+1)*m);
167:     STMatMult(pep->st,i+1,y2,aux);
168:     VecAXPY(y1,a,aux);
169:     VecResetArray(y2);
170:     a *= pep->sfactor;
171:   }
172:   i = deg-2;
173:   VecPlaceArray(y2,py+(i+1)*m);
174:   VecPlaceArray(y3,py+i*m);
175:   VecCopy(y2,aux2);
176:   VecAXPY(aux2,cg[i+1]/ca[i+1],y3);
177:   STMatMult(pep->st,i+1,aux2,aux);
178:   VecAXPY(y1,a,aux);
179:   VecResetArray(y2);
180:   VecResetArray(y3);
181:   a *= pep->sfactor;
182:   i = deg-1;
183:   VecPlaceArray(x1,px+i*m);
184:   VecPlaceArray(y3,py+i*m);
185:   VecCopy(x1,aux2);
186:   VecAXPY(aux2,sigma-cb[i],y3);
187:   VecScale(aux2,1.0/ca[i]);
188:   STMatMult(pep->st,i+1,aux2,aux);
189:   VecAXPY(y1,a,aux);
190:   VecResetArray(x1);
191:   VecResetArray(y3);

193:   VecCopy(y1,aux);
194:   STMatSolve(pep->st,aux,y1);
195:   VecScale(y1,-1.0);

197:   /* final update */
198:   for (i=1;i<deg;i++) {
199:     VecPlaceArray(y2,py+i*m);
200:     tt = t;
201:     t = ((sigma-cb[i-1])*t-cg[i-1]*tp)/ca[i-1]; /* i-th basis polynomial */
202:     tp = tt;
203:     VecAXPY(y2,t,y1);
204:     VecResetArray(y2);
205:   }
206:   VecResetArray(y1);

208:   VecRestoreArrayRead(x,&px);
209:   VecRestoreArray(y,&py);
210:   return(0);
211: }

213: static PetscErrorCode BackTransform_Linear(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
214: {
216:   PEP_LINEAR     *ctx;
217:   ST             stctx;

220:   STShellGetContext(st,&ctx);
221:   PEPGetST(ctx->pep,&stctx);
222:   STBackTransform(stctx,n,eigr,eigi);
223:   return(0);
224: }

226: /*
227:    Dummy backtransform operation
228:  */
229: static PetscErrorCode BackTransform_Skip(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
230: {
232:   return(0);
233: }

235: static PetscErrorCode Apply_Linear(ST st,Vec x,Vec y)
236: {
238:   PEP_LINEAR     *ctx;

241:   STShellGetContext(st,&ctx);
242:   MatMult(ctx->A,x,y);
243:   return(0);
244: }

246: PetscErrorCode PEPSetUp_Linear(PEP pep)
247: {
249:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
250:   ST             st;
251:   PetscInt       i=0,deg=pep->nmat-1;
252:   EPSWhich       which = EPS_LARGEST_MAGNITUDE;
253:   EPSProblemType ptype;
254:   PetscBool      trackall,istrivial,transf,sinv,ks;
255:   PetscScalar    sigma,*epsarray,*peparray;
256:   Vec            veps,w=NULL;
257:   /* function tables */
258:   PetscErrorCode (*fcreate[][2])(MPI_Comm,PEP_LINEAR*,Mat*) = {
259:     { MatCreateExplicit_Linear_NA, MatCreateExplicit_Linear_NB },
260:     { MatCreateExplicit_Linear_SA, MatCreateExplicit_Linear_SB },
261:     { MatCreateExplicit_Linear_HA, MatCreateExplicit_Linear_HB },
262:   };

265:   PEPCheckShiftSinvert(pep);
266:   PEPCheckUnsupported(pep,PEP_FEATURE_STOPPING);
267:   PEPCheckIgnored(pep,PEP_FEATURE_CONVERGENCE);
268:   STGetTransform(pep->st,&transf);
269:   PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
270:   if (!pep->which) { PEPSetWhichEigenpairs_Default(pep); }
271:   if (pep->which==PEP_ALL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
272:   STSetUp(pep->st);
273:   if (!ctx->eps) { PEPLinearGetEPS(pep,&ctx->eps); }
274:   EPSGetST(ctx->eps,&st);
275:   if (!transf && !ctx->usereps) { EPSSetTarget(ctx->eps,pep->target); }
276:   if (sinv && !transf && !ctx->usereps) { STSetDefaultShift(st,pep->target); }
277:   /* compute scale factor if not set by user */
278:   PEPComputeScaleFactor(pep);

280:   if (ctx->explicitmatrix) {
281:     PEPCheckQuadraticCondition(pep,PETSC_TRUE," (with explicit matrix)");
282:     PEPCheckUnsupportedCondition(pep,PEP_FEATURE_NONMONOMIAL,PETSC_TRUE," (with explicit matrix)");
283:     if (transf) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option is not implemented with st-transform flag active");
284:     if (pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Diagonal scaling not allowed in PEPLINEAR with explicit matrices");
285:     if (sinv && !transf) { STSetType(st,STSINVERT); }
286:     RGPushScale(pep->rg,1.0/pep->sfactor);
287:     STGetMatrixTransformed(pep->st,0,&ctx->K);
288:     STGetMatrixTransformed(pep->st,1,&ctx->C);
289:     STGetMatrixTransformed(pep->st,2,&ctx->M);
290:     ctx->sfactor = pep->sfactor;
291:     ctx->dsfactor = pep->dsfactor;

293:     MatDestroy(&ctx->A);
294:     MatDestroy(&ctx->B);
295:     VecDestroy(&ctx->w[0]);
296:     VecDestroy(&ctx->w[1]);
297:     VecDestroy(&ctx->w[2]);
298:     VecDestroy(&ctx->w[3]);

300:     switch (pep->problem_type) {
301:       case PEP_GENERAL:    i = 0; break;
302:       case PEP_HERMITIAN:
303:       case PEP_HYPERBOLIC: i = 1; break;
304:       case PEP_GYROSCOPIC: i = 2; break;
305:     }

307:     (*fcreate[i][0])(PetscObjectComm((PetscObject)pep),ctx,&ctx->A);
308:     (*fcreate[i][1])(PetscObjectComm((PetscObject)pep),ctx,&ctx->B);
309:     PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->A);
310:     PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->B);

312:   } else {   /* implicit matrix */
313:     if (pep->problem_type!=PEP_GENERAL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Must use the explicit matrix option if problem type is not general");
314:     if (!((PetscObject)(ctx->eps))->type_name) {
315:       EPSSetType(ctx->eps,EPSKRYLOVSCHUR);
316:     } else {
317:       PetscObjectTypeCompare((PetscObject)ctx->eps,EPSKRYLOVSCHUR,&ks);
318:       if (!ks) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option only implemented for Krylov-Schur");
319:     }
320:     if (ctx->alpha!=1.0 || ctx->beta!=0.0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option does not support changing alpha,beta parameters of the linearization");
321:     STSetType(st,STSHELL);
322:     STShellSetContext(st,ctx);
323:     if (!transf) { STShellSetBackTransform(st,BackTransform_Linear); }
324:     else { STShellSetBackTransform(st,BackTransform_Skip); }
325:     MatCreateVecsEmpty(pep->A[0],&ctx->w[0],&ctx->w[1]);
326:     MatCreateVecsEmpty(pep->A[0],&ctx->w[2],&ctx->w[3]);
327:     MatCreateVecs(pep->A[0],&ctx->w[4],&ctx->w[5]);
328:     PetscLogObjectParents(pep,6,ctx->w);
329:     MatCreateShell(PetscObjectComm((PetscObject)pep),deg*pep->nloc,deg*pep->nloc,deg*pep->n,deg*pep->n,ctx,&ctx->A);
330:     if (sinv && !transf) {
331:       MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Sinvert);
332:     } else {
333:       MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Shift);
334:     }
335:     STShellSetApply(st,Apply_Linear);
336:     PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->A);
337:     ctx->pep = pep;

339:     PEPBasisCoefficients(pep,pep->pbc);
340:     if (!transf) {
341:       PetscMalloc1(pep->nmat,&pep->solvematcoeffs);
342:       PetscLogObjectMemory((PetscObject)pep,pep->nmat*sizeof(PetscScalar));
343:       if (sinv) {
344:         PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL);
345:       } else {
346:         for (i=0;i<deg;i++) pep->solvematcoeffs[i] = 0.0;
347:         pep->solvematcoeffs[deg] = 1.0;
348:       }
349:       STScaleShift(pep->st,1.0/pep->sfactor);
350:       RGPushScale(pep->rg,1.0/pep->sfactor);
351:     }
352:     if (pep->sfactor!=1.0) {
353:       for (i=0;i<pep->nmat;i++) {
354:         pep->pbc[pep->nmat+i] /= pep->sfactor;
355:         pep->pbc[2*pep->nmat+i] /= pep->sfactor*pep->sfactor;
356:       }
357:     }
358:   }

360:   EPSSetOperators(ctx->eps,ctx->A,ctx->B);
361:   EPSGetProblemType(ctx->eps,&ptype);
362:   if (!ptype) {
363:     if (ctx->explicitmatrix) {
364:       EPSSetProblemType(ctx->eps,EPS_GNHEP);
365:     } else {
366:       EPSSetProblemType(ctx->eps,EPS_NHEP);
367:     }
368:   }
369:   if (!ctx->usereps) {
370:     if (transf) which = EPS_LARGEST_MAGNITUDE;
371:     else {
372:       switch (pep->which) {
373:         case PEP_LARGEST_MAGNITUDE:  which = EPS_LARGEST_MAGNITUDE; break;
374:         case PEP_SMALLEST_MAGNITUDE: which = EPS_SMALLEST_MAGNITUDE; break;
375:         case PEP_LARGEST_REAL:       which = EPS_LARGEST_REAL; break;
376:         case PEP_SMALLEST_REAL:      which = EPS_SMALLEST_REAL; break;
377:         case PEP_LARGEST_IMAGINARY:  which = EPS_LARGEST_IMAGINARY; break;
378:         case PEP_SMALLEST_IMAGINARY: which = EPS_SMALLEST_IMAGINARY; break;
379:         case PEP_TARGET_MAGNITUDE:   which = EPS_TARGET_MAGNITUDE; break;
380:         case PEP_TARGET_REAL:        which = EPS_TARGET_REAL; break;
381:         case PEP_TARGET_IMAGINARY:   which = EPS_TARGET_IMAGINARY; break;
382:         case PEP_ALL:                which = EPS_ALL; break;
383:         case PEP_WHICH_USER:         which = EPS_WHICH_USER;
384:           EPSSetEigenvalueComparison(ctx->eps,pep->sc->comparison,pep->sc->comparisonctx);
385:           break;
386:       }
387:     }
388:     EPSSetWhichEigenpairs(ctx->eps,which);

390:     EPSSetDimensions(ctx->eps,pep->nev,pep->ncv,pep->mpd);
391:     EPSSetTolerances(ctx->eps,SlepcDefaultTol(pep->tol),pep->max_it);
392:   }
393:   RGIsTrivial(pep->rg,&istrivial);
394:   if (!istrivial) {
395:     if (transf) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"PEPLINEAR does not support a nontrivial region with st-transform");
396:     EPSSetRG(ctx->eps,pep->rg);
397:   }
398:   /* Transfer the trackall option from pep to eps */
399:   PEPGetTrackAll(pep,&trackall);
400:   EPSSetTrackAll(ctx->eps,trackall);

402:   /* temporary change of target */
403:   if (pep->sfactor!=1.0) {
404:     EPSGetTarget(ctx->eps,&sigma);
405:     EPSSetTarget(ctx->eps,sigma/pep->sfactor);
406:   }

408:   /* process initial vector */
409:   if (pep->nini<0) {
410:     VecCreateMPI(PetscObjectComm((PetscObject)ctx->eps),deg*pep->nloc,deg*pep->n,&veps);
411:     VecGetArray(veps,&epsarray);
412:     for (i=0;i<deg;i++) {
413:       if (i<-pep->nini) {
414:         VecGetArray(pep->IS[i],&peparray);
415:         PetscArraycpy(epsarray+i*pep->nloc,peparray,pep->nloc);
416:         VecRestoreArray(pep->IS[i],&peparray);
417:       } else {
418:         if (!w) { VecDuplicate(pep->IS[0],&w); }
419:         VecSetRandom(w,NULL);
420:         VecGetArray(w,&peparray);
421:         PetscArraycpy(epsarray+i*pep->nloc,peparray,pep->nloc);
422:         VecRestoreArray(w,&peparray);
423:       }
424:     }
425:     VecRestoreArray(veps,&epsarray);
426:     EPSSetInitialSpace(ctx->eps,1,&veps);
427:     VecDestroy(&veps);
428:     VecDestroy(&w);
429:     SlepcBasisDestroy_Private(&pep->nini,&pep->IS);
430:   }

432:   EPSSetUp(ctx->eps);
433:   EPSGetDimensions(ctx->eps,NULL,&pep->ncv,&pep->mpd);
434:   EPSGetTolerances(ctx->eps,NULL,&pep->max_it);
435:   PEPAllocateSolution(pep,0);
436:   return(0);
437: }

439: /*
440:    PEPLinearExtract_Residual - Auxiliary routine that copies the solution of the
441:    linear eigenproblem to the PEP object. The eigenvector of the generalized
442:    problem is supposed to be
443:                                z = [  x  ]
444:                                    [ l*x ]
445:    The eigenvector is taken from z(1:n) or z(n+1:2*n) depending on the explicitly
446:    computed residual norm.
447:    Finally, x is normalized so that ||x||_2 = 1.
448: */
449: static PetscErrorCode PEPLinearExtract_Residual(PEP pep,EPS eps)
450: {
451:   PetscErrorCode    ierr;
452:   PetscInt          i,k;
453:   const PetscScalar *px;
454:   PetscScalar       *er=pep->eigr,*ei=pep->eigi;
455:   PetscReal         rn1,rn2;
456:   Vec               xr,xi=NULL,wr;
457:   Mat               A;
458: #if !defined(PETSC_USE_COMPLEX)
459:   Vec               wi;
460:   const PetscScalar *py;
461: #endif

464: #if defined(PETSC_USE_COMPLEX)
465:   PEPSetWorkVecs(pep,2);
466: #else
467:   PEPSetWorkVecs(pep,4);
468: #endif
469:   EPSGetOperators(eps,&A,NULL);
470:   MatCreateVecs(A,&xr,NULL);
471:   MatCreateVecsEmpty(pep->A[0],&wr,NULL);
472: #if !defined(PETSC_USE_COMPLEX)
473:   VecDuplicate(xr,&xi);
474:   VecDuplicateEmpty(wr,&wi);
475: #endif
476:   for (i=0;i<pep->nconv;i++) {
477:     EPSGetEigenpair(eps,i,NULL,NULL,xr,xi);
478: #if !defined(PETSC_USE_COMPLEX)
479:     if (ei[i]!=0.0) {   /* complex conjugate pair */
480:       VecGetArrayRead(xr,&px);
481:       VecGetArrayRead(xi,&py);
482:       VecPlaceArray(wr,px);
483:       VecPlaceArray(wi,py);
484:       VecNormalizeComplex(wr,wi,PETSC_TRUE,NULL);
485:       PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,wi,pep->work,&rn1);
486:       BVInsertVec(pep->V,i,wr);
487:       BVInsertVec(pep->V,i+1,wi);
488:       for (k=1;k<pep->nmat-1;k++) {
489:         VecResetArray(wr);
490:         VecResetArray(wi);
491:         VecPlaceArray(wr,px+k*pep->nloc);
492:         VecPlaceArray(wi,py+k*pep->nloc);
493:         VecNormalizeComplex(wr,wi,PETSC_TRUE,NULL);
494:         PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,wi,pep->work,&rn2);
495:         if (rn1>rn2) {
496:           BVInsertVec(pep->V,i,wr);
497:           BVInsertVec(pep->V,i+1,wi);
498:           rn1 = rn2;
499:         }
500:       }
501:       VecResetArray(wr);
502:       VecResetArray(wi);
503:       VecRestoreArrayRead(xr,&px);
504:       VecRestoreArrayRead(xi,&py);
505:       i++;
506:     } else   /* real eigenvalue */
507: #endif
508:     {
509:       VecGetArrayRead(xr,&px);
510:       VecPlaceArray(wr,px);
511:       VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL);
512:       PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn1);
513:       BVInsertVec(pep->V,i,wr);
514:       for (k=1;k<pep->nmat-1;k++) {
515:         VecResetArray(wr);
516:         VecPlaceArray(wr,px+k*pep->nloc);
517:         VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL);
518:         PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn2);
519:         if (rn1>rn2) {
520:           BVInsertVec(pep->V,i,wr);
521:           rn1 = rn2;
522:         }
523:       }
524:       VecResetArray(wr);
525:       VecRestoreArrayRead(xr,&px);
526:     }
527:   }
528:   VecDestroy(&wr);
529:   VecDestroy(&xr);
530: #if !defined(PETSC_USE_COMPLEX)
531:   VecDestroy(&wi);
532:   VecDestroy(&xi);
533: #endif
534:   return(0);
535: }

537: /*
538:    PEPLinearExtract_None - Same as PEPLinearExtract_Norm but always takes
539:    the first block.
540: */
541: static PetscErrorCode PEPLinearExtract_None(PEP pep,EPS eps)
542: {
543:   PetscErrorCode    ierr;
544:   PetscInt          i;
545:   const PetscScalar *px;
546:   Mat               A;
547:   Vec               xr,xi=NULL,w;
548: #if !defined(PETSC_USE_COMPLEX)
549:   PetscScalar       *ei=pep->eigi;
550: #endif

553:   EPSGetOperators(eps,&A,NULL);
554:   MatCreateVecs(A,&xr,NULL);
555: #if !defined(PETSC_USE_COMPLEX)
556:   VecDuplicate(xr,&xi);
557: #endif
558:   MatCreateVecsEmpty(pep->A[0],&w,NULL);
559:   for (i=0;i<pep->nconv;i++) {
560:     EPSGetEigenvector(eps,i,xr,xi);
561: #if !defined(PETSC_USE_COMPLEX)
562:     if (ei[i]!=0.0) {   /* complex conjugate pair */
563:       VecGetArrayRead(xr,&px);
564:       VecPlaceArray(w,px);
565:       BVInsertVec(pep->V,i,w);
566:       VecResetArray(w);
567:       VecRestoreArrayRead(xr,&px);
568:       VecGetArrayRead(xi,&px);
569:       VecPlaceArray(w,px);
570:       BVInsertVec(pep->V,i+1,w);
571:       VecResetArray(w);
572:       VecRestoreArrayRead(xi,&px);
573:       i++;
574:     } else   /* real eigenvalue */
575: #endif
576:     {
577:       VecGetArrayRead(xr,&px);
578:       VecPlaceArray(w,px);
579:       BVInsertVec(pep->V,i,w);
580:       VecResetArray(w);
581:       VecRestoreArrayRead(xr,&px);
582:     }
583:   }
584:   VecDestroy(&w);
585:   VecDestroy(&xr);
586: #if !defined(PETSC_USE_COMPLEX)
587:   VecDestroy(&xi);
588: #endif
589:   return(0);
590: }

592: /*
593:    PEPLinearExtract_Norm - Auxiliary routine that copies the solution of the
594:    linear eigenproblem to the PEP object. The eigenvector of the generalized
595:    problem is supposed to be
596:                                z = [  x  ]
597:                                    [ l*x ]
598:    If |l|<1.0, the eigenvector is taken from z(1:n), otherwise from z(n+1:2*n).
599:    Finally, x is normalized so that ||x||_2 = 1.
600: */
601: static PetscErrorCode PEPLinearExtract_Norm(PEP pep,EPS eps)
602: {
603:   PetscErrorCode    ierr;
604:   PetscInt          i,offset;
605:   const PetscScalar *px;
606:   PetscScalar       *er=pep->eigr;
607:   Mat               A;
608:   Vec               xr,xi=NULL,w;
609: #if !defined(PETSC_USE_COMPLEX)
610:   PetscScalar       *ei=pep->eigi;
611: #endif

614:   EPSGetOperators(eps,&A,NULL);
615:   MatCreateVecs(A,&xr,NULL);
616: #if !defined(PETSC_USE_COMPLEX)
617:   VecDuplicate(xr,&xi);
618: #endif
619:   MatCreateVecsEmpty(pep->A[0],&w,NULL);
620:   for (i=0;i<pep->nconv;i++) {
621:     EPSGetEigenpair(eps,i,NULL,NULL,xr,xi);
622:     if (SlepcAbsEigenvalue(er[i],ei[i])>1.0) offset = (pep->nmat-2)*pep->nloc;
623:     else offset = 0;
624: #if !defined(PETSC_USE_COMPLEX)
625:     if (ei[i]!=0.0) {   /* complex conjugate pair */
626:       VecGetArrayRead(xr,&px);
627:       VecPlaceArray(w,px+offset);
628:       BVInsertVec(pep->V,i,w);
629:       VecResetArray(w);
630:       VecRestoreArrayRead(xr,&px);
631:       VecGetArrayRead(xi,&px);
632:       VecPlaceArray(w,px+offset);
633:       BVInsertVec(pep->V,i+1,w);
634:       VecResetArray(w);
635:       VecRestoreArrayRead(xi,&px);
636:       i++;
637:     } else /* real eigenvalue */
638: #endif
639:     {
640:       VecGetArrayRead(xr,&px);
641:       VecPlaceArray(w,px+offset);
642:       BVInsertVec(pep->V,i,w);
643:       VecResetArray(w);
644:       VecRestoreArrayRead(xr,&px);
645:     }
646:   }
647:   VecDestroy(&w);
648:   VecDestroy(&xr);
649: #if !defined(PETSC_USE_COMPLEX)
650:   VecDestroy(&xi);
651: #endif
652:   return(0);
653: }

655: PetscErrorCode PEPExtractVectors_Linear(PEP pep)
656: {
658:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

661:   switch (pep->extract) {
662:   case PEP_EXTRACT_NONE:
663:     PEPLinearExtract_None(pep,ctx->eps);
664:     break;
665:   case PEP_EXTRACT_NORM:
666:     PEPLinearExtract_Norm(pep,ctx->eps);
667:     break;
668:   case PEP_EXTRACT_RESIDUAL:
669:     PEPLinearExtract_Residual(pep,ctx->eps);
670:     break;
671:   case PEP_EXTRACT_STRUCTURED:
672:     SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Extraction not implemented in this solver");
673:   }
674:   return(0);
675: }

677: PetscErrorCode PEPSolve_Linear(PEP pep)
678: {
680:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
681:   PetscScalar    sigma;
682:   PetscBool      flg;
683:   PetscInt       i;

686:   EPSSolve(ctx->eps);
687:   EPSGetConverged(ctx->eps,&pep->nconv);
688:   EPSGetIterationNumber(ctx->eps,&pep->its);
689:   EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&pep->reason);

691:   /* recover eigenvalues */
692:   for (i=0;i<pep->nconv;i++) {
693:     EPSGetEigenpair(ctx->eps,i,&pep->eigr[i],&pep->eigi[i],NULL,NULL);
694:     pep->eigr[i] *= pep->sfactor;
695:     pep->eigi[i] *= pep->sfactor;
696:   }

698:   /* restore target */
699:   EPSGetTarget(ctx->eps,&sigma);
700:   EPSSetTarget(ctx->eps,sigma*pep->sfactor);

702:   STGetTransform(pep->st,&flg);
703:   if (flg && pep->ops->backtransform) {
704:     (*pep->ops->backtransform)(pep);
705:   }
706:   if (pep->sfactor!=1.0) {
707:     /* Restore original values */
708:     for (i=0;i<pep->nmat;i++) {
709:       pep->pbc[pep->nmat+i] *= pep->sfactor;
710:       pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
711:     }
712:     if (!flg && !ctx->explicitmatrix) {
713:       STScaleShift(pep->st,pep->sfactor);
714:     }
715:   }
716:   if (ctx->explicitmatrix || !flg) {
717:     RGPopScale(pep->rg);
718:   }
719:   return(0);
720: }

722: static PetscErrorCode EPSMonitor_Linear(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
723: {
724:   PEP            pep = (PEP)ctx;

728:   PEPMonitor(pep,its,nconv,eigr,eigi,errest,nest);
729:   return(0);
730: }

732: PetscErrorCode PEPSetFromOptions_Linear(PetscOptionItems *PetscOptionsObject,PEP pep)
733: {
735:   PetscBool      set,val;
736:   PetscInt       k;
737:   PetscReal      array[2]={0,0};
738:   PetscBool      flg;
739:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

742:   PetscOptionsHead(PetscOptionsObject,"PEP Linear Options");

744:     k = 2;
745:     PetscOptionsRealArray("-pep_linear_linearization","Parameters of the linearization","PEPLinearSetLinearization",array,&k,&flg);
746:     if (flg) {
747:       PEPLinearSetLinearization(pep,array[0],array[1]);
748:     }

750:     PetscOptionsBool("-pep_linear_explicitmatrix","Use explicit matrix in linearization","PEPLinearSetExplicitMatrix",ctx->explicitmatrix,&val,&set);
751:     if (set) { PEPLinearSetExplicitMatrix(pep,val); }

753:   PetscOptionsTail();

755:   if (!ctx->eps) { PEPLinearGetEPS(pep,&ctx->eps); }
756:   EPSSetFromOptions(ctx->eps);
757:   return(0);
758: }

760: static PetscErrorCode PEPLinearSetLinearization_Linear(PEP pep,PetscReal alpha,PetscReal beta)
761: {
762:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

765:   if (beta==0.0 && alpha==0.0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
766:   ctx->alpha = alpha;
767:   ctx->beta  = beta;
768:   return(0);
769: }

771: /*@
772:    PEPLinearSetLinearization - Set the coefficients that define
773:    the linearization of a quadratic eigenproblem.

775:    Logically Collective on pep

777:    Input Parameters:
778: +  pep   - polynomial eigenvalue solver
779: .  alpha - first parameter of the linearization
780: -  beta  - second parameter of the linearization

782:    Options Database Key:
783: .  -pep_linear_linearization <alpha,beta> - Sets the coefficients

785:    Notes:
786:    Cannot pass zero for both alpha and beta. The default values are
787:    alpha=1 and beta=0.

789:    Level: advanced

791: .seealso: PEPLinearGetLinearization()
792: @*/
793: PetscErrorCode PEPLinearSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
794: {

801:   PetscTryMethod(pep,"PEPLinearSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
802:   return(0);
803: }

805: static PetscErrorCode PEPLinearGetLinearization_Linear(PEP pep,PetscReal *alpha,PetscReal *beta)
806: {
807:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

810:   if (alpha) *alpha = ctx->alpha;
811:   if (beta)  *beta  = ctx->beta;
812:   return(0);
813: }

815: /*@
816:    PEPLinearGetLinearization - Returns the coefficients that define
817:    the linearization of a quadratic eigenproblem.

819:    Not Collective

821:    Input Parameter:
822: .  pep  - polynomial eigenvalue solver

824:    Output Parameters:
825: +  alpha - the first parameter of the linearization
826: -  beta  - the second parameter of the linearization

828:    Level: advanced

830: .seealso: PEPLinearSetLinearization()
831: @*/
832: PetscErrorCode PEPLinearGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
833: {

838:   PetscUseMethod(pep,"PEPLinearGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
839:   return(0);
840: }

842: static PetscErrorCode PEPLinearSetExplicitMatrix_Linear(PEP pep,PetscBool explicitmatrix)
843: {
844:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

847:   if (ctx->explicitmatrix != explicitmatrix) {
848:     ctx->explicitmatrix = explicitmatrix;
849:     pep->state = PEP_STATE_INITIAL;
850:   }
851:   return(0);
852: }

854: /*@
855:    PEPLinearSetExplicitMatrix - Indicate if the matrices A and B for the
856:    linearization of the problem must be built explicitly.

858:    Logically Collective on pep

860:    Input Parameters:
861: +  pep      - polynomial eigenvalue solver
862: -  explicit - boolean flag indicating if the matrices are built explicitly

864:    Options Database Key:
865: .  -pep_linear_explicitmatrix <boolean> - Indicates the boolean flag

867:    Level: advanced

869: .seealso: PEPLinearGetExplicitMatrix()
870: @*/
871: PetscErrorCode PEPLinearSetExplicitMatrix(PEP pep,PetscBool explicitmatrix)
872: {

878:   PetscTryMethod(pep,"PEPLinearSetExplicitMatrix_C",(PEP,PetscBool),(pep,explicitmatrix));
879:   return(0);
880: }

882: static PetscErrorCode PEPLinearGetExplicitMatrix_Linear(PEP pep,PetscBool *explicitmatrix)
883: {
884:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

887:   *explicitmatrix = ctx->explicitmatrix;
888:   return(0);
889: }

891: /*@
892:    PEPLinearGetExplicitMatrix - Returns the flag indicating if the matrices
893:    A and B for the linearization are built explicitly.

895:    Not Collective

897:    Input Parameter:
898: .  pep  - polynomial eigenvalue solver

900:    Output Parameter:
901: .  explicitmatrix - the mode flag

903:    Level: advanced

905: .seealso: PEPLinearSetExplicitMatrix()
906: @*/
907: PetscErrorCode PEPLinearGetExplicitMatrix(PEP pep,PetscBool *explicitmatrix)
908: {

914:   PetscUseMethod(pep,"PEPLinearGetExplicitMatrix_C",(PEP,PetscBool*),(pep,explicitmatrix));
915:   return(0);
916: }

918: static PetscErrorCode PEPLinearSetEPS_Linear(PEP pep,EPS eps)
919: {
921:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

924:   PetscObjectReference((PetscObject)eps);
925:   EPSDestroy(&ctx->eps);
926:   ctx->eps = eps;
927:   ctx->usereps = PETSC_TRUE;
928:   PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->eps);
929:   pep->state = PEP_STATE_INITIAL;
930:   return(0);
931: }

933: /*@
934:    PEPLinearSetEPS - Associate an eigensolver object (EPS) to the
935:    polynomial eigenvalue solver.

937:    Collective on pep

939:    Input Parameters:
940: +  pep - polynomial eigenvalue solver
941: -  eps - the eigensolver object

943:    Level: advanced

945: .seealso: PEPLinearGetEPS()
946: @*/
947: PetscErrorCode PEPLinearSetEPS(PEP pep,EPS eps)
948: {

955:   PetscTryMethod(pep,"PEPLinearSetEPS_C",(PEP,EPS),(pep,eps));
956:   return(0);
957: }

959: static PetscErrorCode PEPLinearGetEPS_Linear(PEP pep,EPS *eps)
960: {
962:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

965:   if (!ctx->eps) {
966:     EPSCreate(PetscObjectComm((PetscObject)pep),&ctx->eps);
967:     PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)pep,1);
968:     EPSSetOptionsPrefix(ctx->eps,((PetscObject)pep)->prefix);
969:     EPSAppendOptionsPrefix(ctx->eps,"pep_linear_");
970:     PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->eps);
971:     PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)pep)->options);
972:     EPSMonitorSet(ctx->eps,EPSMonitor_Linear,pep,NULL);
973:   }
974:   *eps = ctx->eps;
975:   return(0);
976: }

978: /*@
979:    PEPLinearGetEPS - Retrieve the eigensolver object (EPS) associated
980:    to the polynomial eigenvalue solver.

982:    Not Collective

984:    Input Parameter:
985: .  pep - polynomial eigenvalue solver

987:    Output Parameter:
988: .  eps - the eigensolver object

990:    Level: advanced

992: .seealso: PEPLinearSetEPS()
993: @*/
994: PetscErrorCode PEPLinearGetEPS(PEP pep,EPS *eps)
995: {

1001:   PetscUseMethod(pep,"PEPLinearGetEPS_C",(PEP,EPS*),(pep,eps));
1002:   return(0);
1003: }

1005: PetscErrorCode PEPView_Linear(PEP pep,PetscViewer viewer)
1006: {
1008:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
1009:   PetscBool      isascii;

1012:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1013:   if (isascii) {
1014:     if (!ctx->eps) { PEPLinearGetEPS(pep,&ctx->eps); }
1015:     PetscViewerASCIIPrintf(viewer,"  %s matrices\n",ctx->explicitmatrix? "explicit": "implicit");
1016:     PetscViewerASCIIPrintf(viewer,"  linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta);
1017:     PetscViewerASCIIPushTab(viewer);
1018:     EPSView(ctx->eps,viewer);
1019:     PetscViewerASCIIPopTab(viewer);
1020:   }
1021:   return(0);
1022: }

1024: PetscErrorCode PEPReset_Linear(PEP pep)
1025: {
1027:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

1030:   if (!ctx->eps) { EPSReset(ctx->eps); }
1031:   MatDestroy(&ctx->A);
1032:   MatDestroy(&ctx->B);
1033:   VecDestroy(&ctx->w[0]);
1034:   VecDestroy(&ctx->w[1]);
1035:   VecDestroy(&ctx->w[2]);
1036:   VecDestroy(&ctx->w[3]);
1037:   VecDestroy(&ctx->w[4]);
1038:   VecDestroy(&ctx->w[5]);
1039:   return(0);
1040: }

1042: PetscErrorCode PEPDestroy_Linear(PEP pep)
1043: {
1045:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

1048:   EPSDestroy(&ctx->eps);
1049:   PetscFree(pep->data);
1050:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",NULL);
1051:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",NULL);
1052:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",NULL);
1053:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",NULL);
1054:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",NULL);
1055:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",NULL);
1056:   return(0);
1057: }

1059: SLEPC_EXTERN PetscErrorCode PEPCreate_Linear(PEP pep)
1060: {
1062:   PEP_LINEAR     *ctx;

1065:   PetscNewLog(pep,&ctx);
1066:   pep->data = (void*)ctx;

1068:   pep->lineariz       = PETSC_TRUE;
1069:   ctx->explicitmatrix = PETSC_FALSE;
1070:   ctx->alpha          = 1.0;
1071:   ctx->beta           = 0.0;

1073:   pep->ops->solve          = PEPSolve_Linear;
1074:   pep->ops->setup          = PEPSetUp_Linear;
1075:   pep->ops->setfromoptions = PEPSetFromOptions_Linear;
1076:   pep->ops->destroy        = PEPDestroy_Linear;
1077:   pep->ops->reset          = PEPReset_Linear;
1078:   pep->ops->view           = PEPView_Linear;
1079:   pep->ops->backtransform  = PEPBackTransform_Default;
1080:   pep->ops->computevectors = PEPComputeVectors_Default;
1081:   pep->ops->extractvectors = PEPExtractVectors_Linear;

1083:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",PEPLinearSetLinearization_Linear);
1084:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",PEPLinearGetLinearization_Linear);
1085:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",PEPLinearSetEPS_Linear);
1086:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",PEPLinearGetEPS_Linear);
1087:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",PEPLinearSetExplicitMatrix_Linear);
1088:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",PEPLinearGetExplicitMatrix_Linear);
1089:   return(0);
1090: }