Actual source code: linear.c

slepc-3.9.2 2018-07-02
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2018, 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>         /*I "slepcpep.h" I*/
 15: #include "linearp.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,(void**)&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,(void**)&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,(void**)&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,(void**)&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,shift,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_N1A, MatCreateExplicit_Linear_N1B },   /* N1 */
260:     { MatCreateExplicit_Linear_N2A, MatCreateExplicit_Linear_N2B },   /* N2 */
261:     { MatCreateExplicit_Linear_S1A, MatCreateExplicit_Linear_S1B },   /* S1 */
262:     { MatCreateExplicit_Linear_S2A, MatCreateExplicit_Linear_S2B },   /* S2 */
263:     { MatCreateExplicit_Linear_H1A, MatCreateExplicit_Linear_H1B },   /* H1 */
264:     { MatCreateExplicit_Linear_H2A, MatCreateExplicit_Linear_H2B }    /* H2 */
265:   };

268:   if (pep->stopping!=PEPStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"User-defined stopping test not supported");
269:   pep->lineariz = PETSC_TRUE;
270:   if (!ctx->cform) ctx->cform = 1;
271:   STGetTransform(pep->st,&transf);
272:   PetscObjectTypeCompare((PetscObject)pep->st,STSHIFT,&shift);
273:   PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
274:   if (!shift && !sinv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Only STSHIFT and STSINVERT spectral transformations can be used");
275:   if (!pep->which) {
276:     if (sinv) pep->which = PEP_TARGET_MAGNITUDE;
277:     else pep->which = PEP_LARGEST_MAGNITUDE;
278:   }
279:   if (pep->which==PEP_ALL) SETERRQ(PetscObjectComm((PetscObject)pep),1,"Wrong value of pep->which");
280:   STSetUp(pep->st);
281:   if (!ctx->eps) { PEPLinearGetEPS(pep,&ctx->eps); }
282:   EPSGetST(ctx->eps,&st);
283:   if (!transf && !ctx->usereps) { EPSSetTarget(ctx->eps,pep->target); }
284:   if (sinv && !transf && !ctx->usereps) { STSetDefaultShift(st,pep->target); }
285:   /* compute scale factor if not set by user */
286:   PEPComputeScaleFactor(pep);

288:   if (ctx->explicitmatrix) {
289:     if (transf) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option is not implemented with st-transform flag active");
290:     if (pep->nmat!=3) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option only available for quadratic problems");
291:     if (pep->basis!=PEP_BASIS_MONOMIAL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option not implemented for non-monomial bases");
292:     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");
293:     if (sinv && !transf) { STSetType(st,STSINVERT); }
294:     RGPushScale(pep->rg,1.0/pep->sfactor);
295:     STGetMatrixTransformed(pep->st,0,&ctx->K);
296:     STGetMatrixTransformed(pep->st,1,&ctx->C);
297:     STGetMatrixTransformed(pep->st,2,&ctx->M);
298:     ctx->sfactor = pep->sfactor;
299:     ctx->dsfactor = pep->dsfactor;

301:     MatDestroy(&ctx->A);
302:     MatDestroy(&ctx->B);
303:     VecDestroy(&ctx->w[0]);
304:     VecDestroy(&ctx->w[1]);
305:     VecDestroy(&ctx->w[2]);
306:     VecDestroy(&ctx->w[3]);

308:     switch (pep->problem_type) {
309:       case PEP_GENERAL:    i = 0; break;
310:       case PEP_HERMITIAN:
311:       case PEP_HYPERBOLIC: i = 2; break;
312:       case PEP_GYROSCOPIC: i = 4; break;
313:     }
314:     i += ctx->cform-1;

316:     (*fcreate[i][0])(PetscObjectComm((PetscObject)pep),ctx,&ctx->A);
317:     (*fcreate[i][1])(PetscObjectComm((PetscObject)pep),ctx,&ctx->B);
318:     PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->A);
319:     PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->B);

321:   } else {   /* implicit matrix */
322:     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");
323:     if (!((PetscObject)(ctx->eps))->type_name) {
324:       EPSSetType(ctx->eps,EPSKRYLOVSCHUR);
325:     } else {
326:       PetscObjectTypeCompare((PetscObject)ctx->eps,EPSKRYLOVSCHUR,&ks);
327:       if (!ks) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option only implemented for Krylov-Schur");
328:     }
329:     if (ctx->cform!=1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option not available for 2nd companion form");
330:     STSetType(st,STSHELL);
331:     STShellSetContext(st,(PetscObject)ctx);
332:     if (!transf) { STShellSetBackTransform(st,BackTransform_Linear); }
333:     else { STShellSetBackTransform(st,BackTransform_Skip); }
334:     MatCreateVecsEmpty(pep->A[0],&ctx->w[0],&ctx->w[1]);
335:     MatCreateVecsEmpty(pep->A[0],&ctx->w[2],&ctx->w[3]);
336:     MatCreateVecs(pep->A[0],&ctx->w[4],&ctx->w[5]);
337:     PetscLogObjectParents(pep,6,ctx->w);
338:     MatCreateShell(PetscObjectComm((PetscObject)pep),deg*pep->nloc,deg*pep->nloc,deg*pep->n,deg*pep->n,ctx,&ctx->A);
339:     if (sinv && !transf) {
340:       MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Sinvert);
341:     } else {
342:       MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Shift);
343:     }
344:     STShellSetApply(st,Apply_Linear);
345:     PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->A);
346:     ctx->pep = pep;

348:     PEPBasisCoefficients(pep,pep->pbc);
349:     if (!transf) {
350:       PetscMalloc1(pep->nmat,&pep->solvematcoeffs);
351:       PetscLogObjectMemory((PetscObject)pep,pep->nmat*sizeof(PetscScalar));
352:       if (sinv) {
353:         PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL);
354:       } else {
355:         for (i=0;i<deg;i++) pep->solvematcoeffs[i] = 0.0;
356:         pep->solvematcoeffs[deg] = 1.0;
357:       }
358:       STScaleShift(pep->st,1.0/pep->sfactor);
359:       RGPushScale(pep->rg,1.0/pep->sfactor);
360:     }
361:     if (pep->sfactor!=1.0) {
362:       for (i=0;i<pep->nmat;i++) {
363:         pep->pbc[pep->nmat+i] /= pep->sfactor;
364:         pep->pbc[2*pep->nmat+i] /= pep->sfactor*pep->sfactor;
365:       }
366:     }
367:   }

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

399:     EPSSetDimensions(ctx->eps,pep->nev,pep->ncv?pep->ncv:PETSC_DEFAULT,pep->mpd?pep->mpd:PETSC_DEFAULT);
400:     EPSSetTolerances(ctx->eps,pep->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL:pep->tol,pep->max_it?pep->max_it:PETSC_DEFAULT);
401:   }
402:   RGIsTrivial(pep->rg,&istrivial);
403:   if (!istrivial) {
404:     if (transf) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"PEPLINEAR does not support a nontrivial region with st-transform");
405:     EPSSetRG(ctx->eps,pep->rg);
406:   }
407:   /* Transfer the trackall option from pep to eps */
408:   PEPGetTrackAll(pep,&trackall);
409:   EPSSetTrackAll(ctx->eps,trackall);

411:   /* temporary change of target */
412:   if (pep->sfactor!=1.0) {
413:     EPSGetTarget(ctx->eps,&sigma);
414:     EPSSetTarget(ctx->eps,sigma/pep->sfactor);
415:   }

417:   /* process initial vector */
418:   if (pep->nini<0) {
419:     VecCreateMPI(PetscObjectComm((PetscObject)ctx->eps),deg*pep->nloc,deg*pep->n,&veps);
420:     VecGetArray(veps,&epsarray);
421:     for (i=0;i<deg;i++) {
422:       if (i<-pep->nini) {
423:         VecGetArray(pep->IS[i],&peparray);
424:         PetscMemcpy(epsarray+i*pep->nloc,peparray,pep->nloc*sizeof(PetscScalar));
425:         VecRestoreArray(pep->IS[i],&peparray);
426:       } else {
427:         if (!w) { VecDuplicate(pep->IS[0],&w); }
428:         VecSetRandom(w,NULL);
429:         VecGetArray(w,&peparray);
430:         PetscMemcpy(epsarray+i*pep->nloc,peparray,pep->nloc*sizeof(PetscScalar));
431:         VecRestoreArray(w,&peparray);
432:       }
433:     }
434:     VecRestoreArray(veps,&epsarray);
435:     EPSSetInitialSpace(ctx->eps,1,&veps);
436:     VecDestroy(&veps);
437:     VecDestroy(&w);
438:     SlepcBasisDestroy_Private(&pep->nini,&pep->IS);
439:   }

441:   EPSSetUp(ctx->eps);
442:   EPSGetDimensions(ctx->eps,NULL,&pep->ncv,&pep->mpd);
443:   EPSGetTolerances(ctx->eps,NULL,&pep->max_it);
444:   PEPAllocateSolution(pep,0);
445:   return(0);
446: }

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

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

546: /*
547:    PEPLinearExtract_None - Same as PEPLinearExtract_Norm but always takes
548:    the first block.
549: */
550: static PetscErrorCode PEPLinearExtract_None(PEP pep,EPS eps)
551: {
552:   PetscErrorCode    ierr;
553:   PetscInt          i;
554:   const PetscScalar *px;
555:   Mat               A;
556:   Vec               xr,xi,w;
557: #if !defined(PETSC_USE_COMPLEX)
558:   PetscScalar       *ei=pep->eigi;
559: #endif

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

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

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

660: PetscErrorCode PEPExtractVectors_Linear(PEP pep)
661: {
663:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

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

682: PetscErrorCode PEPSolve_Linear(PEP pep)
683: {
685:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
686:   PetscScalar    sigma;
687:   PetscBool      flg;
688:   PetscInt       i;

691:   EPSSolve(ctx->eps);
692:   EPSGetConverged(ctx->eps,&pep->nconv);
693:   EPSGetIterationNumber(ctx->eps,&pep->its);
694:   EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&pep->reason);

696:   /* recover eigenvalues */
697:   for (i=0;i<pep->nconv;i++) {
698:     EPSGetEigenpair(ctx->eps,i,&pep->eigr[i],&pep->eigi[i],NULL,NULL);
699:     pep->eigr[i] *= pep->sfactor;
700:     pep->eigi[i] *= pep->sfactor;
701:   }

703:   /* restore target */
704:   EPSGetTarget(ctx->eps,&sigma);
705:   EPSSetTarget(ctx->eps,sigma*pep->sfactor);

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

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

733:   PEPMonitor(pep,its,nconv,eigr,eigi,errest,nest);
734:   return(0);
735: }

737: PetscErrorCode PEPSetFromOptions_Linear(PetscOptionItems *PetscOptionsObject,PEP pep)
738: {
740:   PetscBool      set,val;
741:   PetscInt       i;
742:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

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

747:     PetscOptionsInt("-pep_linear_cform","Number of the companion form (1 or 2)","PEPLinearSetCompanionForm",ctx->cform,&i,&set);
748:     if (set) { PEPLinearSetCompanionForm(pep,i); }

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 PEPLinearSetCompanionForm_Linear(PEP pep,PetscInt cform)
761: {
762:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

765:   if (!cform) return(0);
766:   if (cform==PETSC_DECIDE || cform==PETSC_DEFAULT) ctx->cform = 1;
767:   else {
768:     if (cform!=1 && cform!=2) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid value of argument 'cform'");
769:     ctx->cform = cform;
770:   }
771:   return(0);
772: }

774: /*@
775:    PEPLinearSetCompanionForm - Choose between the two companion forms available
776:    for the linearization of a quadratic eigenproblem.

778:    Logically Collective on PEP

780:    Input Parameters:
781: +  pep   - polynomial eigenvalue solver
782: -  cform - 1 or 2 (first or second companion form)

784:    Options Database Key:
785: .  -pep_linear_cform <int> - Choose the companion form

787:    Level: advanced

789: .seealso: PEPLinearGetCompanionForm()
790: @*/
791: PetscErrorCode PEPLinearSetCompanionForm(PEP pep,PetscInt cform)
792: {

798:   PetscTryMethod(pep,"PEPLinearSetCompanionForm_C",(PEP,PetscInt),(pep,cform));
799:   return(0);
800: }

802: static PetscErrorCode PEPLinearGetCompanionForm_Linear(PEP pep,PetscInt *cform)
803: {
804:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

807:   *cform = ctx->cform;
808:   return(0);
809: }

811: /*@
812:    PEPLinearGetCompanionForm - Returns the number of the companion form that
813:    will be used for the linearization of a quadratic eigenproblem.

815:    Not Collective

817:    Input Parameter:
818: .  pep  - polynomial eigenvalue solver

820:    Output Parameter:
821: .  cform - the companion form number (1 or 2)

823:    Level: advanced

825: .seealso: PEPLinearSetCompanionForm()
826: @*/
827: PetscErrorCode PEPLinearGetCompanionForm(PEP pep,PetscInt *cform)
828: {

834:   PetscUseMethod(pep,"PEPLinearGetCompanionForm_C",(PEP,PetscInt*),(pep,cform));
835:   return(0);
836: }

838: static PetscErrorCode PEPLinearSetExplicitMatrix_Linear(PEP pep,PetscBool explicitmatrix)
839: {
840:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

843:   if (ctx->explicitmatrix != explicitmatrix) {
844:     ctx->explicitmatrix = explicitmatrix;
845:     pep->state = PEP_STATE_INITIAL;
846:   }
847:   return(0);
848: }

850: /*@
851:    PEPLinearSetExplicitMatrix - Indicate if the matrices A and B for the
852:    linearization of the problem must be built explicitly.

854:    Logically Collective on PEP

856:    Input Parameters:
857: +  pep      - polynomial eigenvalue solver
858: -  explicit - boolean flag indicating if the matrices are built explicitly

860:    Options Database Key:
861: .  -pep_linear_explicitmatrix <boolean> - Indicates the boolean flag

863:    Level: advanced

865: .seealso: PEPLinearGetExplicitMatrix()
866: @*/
867: PetscErrorCode PEPLinearSetExplicitMatrix(PEP pep,PetscBool explicitmatrix)
868: {

874:   PetscTryMethod(pep,"PEPLinearSetExplicitMatrix_C",(PEP,PetscBool),(pep,explicitmatrix));
875:   return(0);
876: }

878: static PetscErrorCode PEPLinearGetExplicitMatrix_Linear(PEP pep,PetscBool *explicitmatrix)
879: {
880:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

883:   *explicitmatrix = ctx->explicitmatrix;
884:   return(0);
885: }

887: /*@
888:    PEPLinearGetExplicitMatrix - Returns the flag indicating if the matrices
889:    A and B for the linearization are built explicitly.

891:    Not Collective

893:    Input Parameter:
894: .  pep  - polynomial eigenvalue solver

896:    Output Parameter:
897: .  explicitmatrix - the mode flag

899:    Level: advanced

901: .seealso: PEPLinearSetExplicitMatrix()
902: @*/
903: PetscErrorCode PEPLinearGetExplicitMatrix(PEP pep,PetscBool *explicitmatrix)
904: {

910:   PetscUseMethod(pep,"PEPLinearGetExplicitMatrix_C",(PEP,PetscBool*),(pep,explicitmatrix));
911:   return(0);
912: }

914: static PetscErrorCode PEPLinearSetEPS_Linear(PEP pep,EPS eps)
915: {
917:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

920:   PetscObjectReference((PetscObject)eps);
921:   EPSDestroy(&ctx->eps);
922:   ctx->eps = eps;
923:   ctx->usereps = PETSC_TRUE;
924:   PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->eps);
925:   pep->state = PEP_STATE_INITIAL;
926:   return(0);
927: }

929: /*@
930:    PEPLinearSetEPS - Associate an eigensolver object (EPS) to the
931:    polynomial eigenvalue solver.

933:    Collective on PEP

935:    Input Parameters:
936: +  pep - polynomial eigenvalue solver
937: -  eps - the eigensolver object

939:    Level: advanced

941: .seealso: PEPLinearGetEPS()
942: @*/
943: PetscErrorCode PEPLinearSetEPS(PEP pep,EPS eps)
944: {

951:   PetscTryMethod(pep,"PEPLinearSetEPS_C",(PEP,EPS),(pep,eps));
952:   return(0);
953: }

955: static PetscErrorCode PEPLinearGetEPS_Linear(PEP pep,EPS *eps)
956: {
958:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

961:   if (!ctx->eps) {
962:     EPSCreate(PetscObjectComm((PetscObject)pep),&ctx->eps);
963:     PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)pep,1);
964:     EPSSetOptionsPrefix(ctx->eps,((PetscObject)pep)->prefix);
965:     EPSAppendOptionsPrefix(ctx->eps,"pep_linear_");
966:     PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->eps);
967:     EPSMonitorSet(ctx->eps,EPSMonitor_Linear,pep,NULL);
968:   }
969:   *eps = ctx->eps;
970:   return(0);
971: }

973: /*@
974:    PEPLinearGetEPS - Retrieve the eigensolver object (EPS) associated
975:    to the polynomial eigenvalue solver.

977:    Not Collective

979:    Input Parameter:
980: .  pep - polynomial eigenvalue solver

982:    Output Parameter:
983: .  eps - the eigensolver object

985:    Level: advanced

987: .seealso: PEPLinearSetEPS()
988: @*/
989: PetscErrorCode PEPLinearGetEPS(PEP pep,EPS *eps)
990: {

996:   PetscUseMethod(pep,"PEPLinearGetEPS_C",(PEP,EPS*),(pep,eps));
997:   return(0);
998: }

1000: PetscErrorCode PEPView_Linear(PEP pep,PetscViewer viewer)
1001: {
1003:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
1004:   PetscBool      isascii;

1007:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1008:   if (isascii) {
1009:     if (!ctx->eps) { PEPLinearGetEPS(pep,&ctx->eps); }
1010:     PetscViewerASCIIPrintf(viewer,"  %s matrices\n",ctx->explicitmatrix? "explicit": "implicit");
1011:     PetscViewerASCIIPrintf(viewer,"  %s companion form\n",ctx->cform==1? "1st": "2nd");
1012:     EPSView(ctx->eps,viewer);
1013:   }
1014:   return(0);
1015: }

1017: PetscErrorCode PEPReset_Linear(PEP pep)
1018: {
1020:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

1023:   if (!ctx->eps) { EPSReset(ctx->eps); }
1024:   MatDestroy(&ctx->A);
1025:   MatDestroy(&ctx->B);
1026:   VecDestroy(&ctx->w[0]);
1027:   VecDestroy(&ctx->w[1]);
1028:   VecDestroy(&ctx->w[2]);
1029:   VecDestroy(&ctx->w[3]);
1030:   VecDestroy(&ctx->w[4]);
1031:   VecDestroy(&ctx->w[5]);
1032:   return(0);
1033: }

1035: PetscErrorCode PEPDestroy_Linear(PEP pep)
1036: {
1038:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

1041:   EPSDestroy(&ctx->eps);
1042:   PetscFree(pep->data);
1043:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetCompanionForm_C",NULL);
1044:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetCompanionForm_C",NULL);
1045:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",NULL);
1046:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",NULL);
1047:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",NULL);
1048:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",NULL);
1049:   return(0);
1050: }

1052: PETSC_EXTERN PetscErrorCode PEPCreate_Linear(PEP pep)
1053: {
1055:   PEP_LINEAR     *ctx;

1058:   PetscNewLog(pep,&ctx);
1059:   ctx->explicitmatrix = PETSC_FALSE;
1060:   pep->data = (void*)ctx;

1062:   pep->ops->solve          = PEPSolve_Linear;
1063:   pep->ops->setup          = PEPSetUp_Linear;
1064:   pep->ops->setfromoptions = PEPSetFromOptions_Linear;
1065:   pep->ops->destroy        = PEPDestroy_Linear;
1066:   pep->ops->reset          = PEPReset_Linear;
1067:   pep->ops->view           = PEPView_Linear;
1068:   pep->ops->backtransform  = PEPBackTransform_Default;
1069:   pep->ops->computevectors = PEPComputeVectors_Default;
1070:   pep->ops->extractvectors = PEPExtractVectors_Linear;

1072:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetCompanionForm_C",PEPLinearSetCompanionForm_Linear);
1073:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetCompanionForm_C",PEPLinearGetCompanionForm_Linear);
1074:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",PEPLinearSetEPS_Linear);
1075:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",PEPLinearGetEPS_Linear);
1076:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",PEPLinearSetExplicitMatrix_Linear);
1077:   PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",PEPLinearGetExplicitMatrix_Linear);
1078:   return(0);
1079: }