Actual source code: linear.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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_NA, MatCreateExplicit_Linear_NB },
260: { MatCreateExplicit_Linear_SA, MatCreateExplicit_Linear_SB },
261: { MatCreateExplicit_Linear_HA, MatCreateExplicit_Linear_HB },
262: };
265: if (pep->stopping!=PEPStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"User-defined stopping test not supported");
266: pep->lineariz = PETSC_TRUE;
267: STGetTransform(pep->st,&transf);
268: PetscObjectTypeCompare((PetscObject)pep->st,STSHIFT,&shift);
269: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
270: if (!shift && !sinv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Only STSHIFT and STSINVERT spectral transformations can be used");
271: if (!pep->which) {
272: if (sinv) pep->which = PEP_TARGET_MAGNITUDE;
273: else pep->which = PEP_LARGEST_MAGNITUDE;
274: }
275: if (pep->which==PEP_ALL) SETERRQ(PetscObjectComm((PetscObject)pep),1,"Wrong value of pep->which");
276: STSetUp(pep->st);
277: if (!ctx->eps) { PEPLinearGetEPS(pep,&ctx->eps); }
278: EPSGetST(ctx->eps,&st);
279: if (!transf && !ctx->usereps) { EPSSetTarget(ctx->eps,pep->target); }
280: if (sinv && !transf && !ctx->usereps) { STSetDefaultShift(st,pep->target); }
281: /* compute scale factor if not set by user */
282: PEPComputeScaleFactor(pep);
284: if (ctx->explicitmatrix) {
285: if (transf) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option is not implemented with st-transform flag active");
286: if (pep->nmat!=3) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option only available for quadratic problems");
287: if (pep->basis!=PEP_BASIS_MONOMIAL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option not implemented for non-monomial bases");
288: 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");
289: if (sinv && !transf) { STSetType(st,STSINVERT); }
290: RGPushScale(pep->rg,1.0/pep->sfactor);
291: STGetMatrixTransformed(pep->st,0,&ctx->K);
292: STGetMatrixTransformed(pep->st,1,&ctx->C);
293: STGetMatrixTransformed(pep->st,2,&ctx->M);
294: ctx->sfactor = pep->sfactor;
295: ctx->dsfactor = pep->dsfactor;
297: MatDestroy(&ctx->A);
298: MatDestroy(&ctx->B);
299: VecDestroy(&ctx->w[0]);
300: VecDestroy(&ctx->w[1]);
301: VecDestroy(&ctx->w[2]);
302: VecDestroy(&ctx->w[3]);
304: switch (pep->problem_type) {
305: case PEP_GENERAL: i = 0; break;
306: case PEP_HERMITIAN:
307: case PEP_HYPERBOLIC: i = 1; break;
308: case PEP_GYROSCOPIC: i = 2; break;
309: }
311: (*fcreate[i][0])(PetscObjectComm((PetscObject)pep),ctx,&ctx->A);
312: (*fcreate[i][1])(PetscObjectComm((PetscObject)pep),ctx,&ctx->B);
313: PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->A);
314: PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->B);
316: } else { /* implicit matrix */
317: 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");
318: if (!((PetscObject)(ctx->eps))->type_name) {
319: EPSSetType(ctx->eps,EPSKRYLOVSCHUR);
320: } else {
321: PetscObjectTypeCompare((PetscObject)ctx->eps,EPSKRYLOVSCHUR,&ks);
322: if (!ks) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option only implemented for Krylov-Schur");
323: }
324: 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");
325: STSetType(st,STSHELL);
326: STShellSetContext(st,(PetscObject)ctx);
327: if (!transf) { STShellSetBackTransform(st,BackTransform_Linear); }
328: else { STShellSetBackTransform(st,BackTransform_Skip); }
329: MatCreateVecsEmpty(pep->A[0],&ctx->w[0],&ctx->w[1]);
330: MatCreateVecsEmpty(pep->A[0],&ctx->w[2],&ctx->w[3]);
331: MatCreateVecs(pep->A[0],&ctx->w[4],&ctx->w[5]);
332: PetscLogObjectParents(pep,6,ctx->w);
333: MatCreateShell(PetscObjectComm((PetscObject)pep),deg*pep->nloc,deg*pep->nloc,deg*pep->n,deg*pep->n,ctx,&ctx->A);
334: if (sinv && !transf) {
335: MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Sinvert);
336: } else {
337: MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Shift);
338: }
339: STShellSetApply(st,Apply_Linear);
340: PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->A);
341: ctx->pep = pep;
343: PEPBasisCoefficients(pep,pep->pbc);
344: if (!transf) {
345: PetscMalloc1(pep->nmat,&pep->solvematcoeffs);
346: PetscLogObjectMemory((PetscObject)pep,pep->nmat*sizeof(PetscScalar));
347: if (sinv) {
348: PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL);
349: } else {
350: for (i=0;i<deg;i++) pep->solvematcoeffs[i] = 0.0;
351: pep->solvematcoeffs[deg] = 1.0;
352: }
353: STScaleShift(pep->st,1.0/pep->sfactor);
354: RGPushScale(pep->rg,1.0/pep->sfactor);
355: }
356: if (pep->sfactor!=1.0) {
357: for (i=0;i<pep->nmat;i++) {
358: pep->pbc[pep->nmat+i] /= pep->sfactor;
359: pep->pbc[2*pep->nmat+i] /= pep->sfactor*pep->sfactor;
360: }
361: }
362: }
364: EPSSetOperators(ctx->eps,ctx->A,ctx->B);
365: EPSGetProblemType(ctx->eps,&ptype);
366: if (!ptype) {
367: if (ctx->explicitmatrix) {
368: EPSSetProblemType(ctx->eps,EPS_GNHEP);
369: } else {
370: EPSSetProblemType(ctx->eps,EPS_NHEP);
371: }
372: }
373: if (!ctx->usereps) {
374: if (transf) which = EPS_LARGEST_MAGNITUDE;
375: else {
376: switch (pep->which) {
377: case PEP_LARGEST_MAGNITUDE: which = EPS_LARGEST_MAGNITUDE; break;
378: case PEP_SMALLEST_MAGNITUDE: which = EPS_SMALLEST_MAGNITUDE; break;
379: case PEP_LARGEST_REAL: which = EPS_LARGEST_REAL; break;
380: case PEP_SMALLEST_REAL: which = EPS_SMALLEST_REAL; break;
381: case PEP_LARGEST_IMAGINARY: which = EPS_LARGEST_IMAGINARY; break;
382: case PEP_SMALLEST_IMAGINARY: which = EPS_SMALLEST_IMAGINARY; break;
383: case PEP_TARGET_MAGNITUDE: which = EPS_TARGET_MAGNITUDE; break;
384: case PEP_TARGET_REAL: which = EPS_TARGET_REAL; break;
385: case PEP_TARGET_IMAGINARY: which = EPS_TARGET_IMAGINARY; break;
386: case PEP_ALL: which = EPS_ALL; break;
387: case PEP_WHICH_USER: which = EPS_WHICH_USER;
388: EPSSetEigenvalueComparison(ctx->eps,pep->sc->comparison,pep->sc->comparisonctx);
389: break;
390: }
391: }
392: EPSSetWhichEigenpairs(ctx->eps,which);
394: EPSSetDimensions(ctx->eps,pep->nev,pep->ncv?pep->ncv:PETSC_DEFAULT,pep->mpd?pep->mpd:PETSC_DEFAULT);
395: EPSSetTolerances(ctx->eps,pep->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL:pep->tol,pep->max_it?pep->max_it:PETSC_DEFAULT);
396: }
397: RGIsTrivial(pep->rg,&istrivial);
398: if (!istrivial) {
399: if (transf) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"PEPLINEAR does not support a nontrivial region with st-transform");
400: EPSSetRG(ctx->eps,pep->rg);
401: }
402: /* Transfer the trackall option from pep to eps */
403: PEPGetTrackAll(pep,&trackall);
404: EPSSetTrackAll(ctx->eps,trackall);
406: /* temporary change of target */
407: if (pep->sfactor!=1.0) {
408: EPSGetTarget(ctx->eps,&sigma);
409: EPSSetTarget(ctx->eps,sigma/pep->sfactor);
410: }
412: /* process initial vector */
413: if (pep->nini<0) {
414: VecCreateMPI(PetscObjectComm((PetscObject)ctx->eps),deg*pep->nloc,deg*pep->n,&veps);
415: VecGetArray(veps,&epsarray);
416: for (i=0;i<deg;i++) {
417: if (i<-pep->nini) {
418: VecGetArray(pep->IS[i],&peparray);
419: PetscMemcpy(epsarray+i*pep->nloc,peparray,pep->nloc*sizeof(PetscScalar));
420: VecRestoreArray(pep->IS[i],&peparray);
421: } else {
422: if (!w) { VecDuplicate(pep->IS[0],&w); }
423: VecSetRandom(w,NULL);
424: VecGetArray(w,&peparray);
425: PetscMemcpy(epsarray+i*pep->nloc,peparray,pep->nloc*sizeof(PetscScalar));
426: VecRestoreArray(w,&peparray);
427: }
428: }
429: VecRestoreArray(veps,&epsarray);
430: EPSSetInitialSpace(ctx->eps,1,&veps);
431: VecDestroy(&veps);
432: VecDestroy(&w);
433: SlepcBasisDestroy_Private(&pep->nini,&pep->IS);
434: }
436: EPSSetUp(ctx->eps);
437: EPSGetDimensions(ctx->eps,NULL,&pep->ncv,&pep->mpd);
438: EPSGetTolerances(ctx->eps,NULL,&pep->max_it);
439: PEPAllocateSolution(pep,0);
440: return(0);
441: }
443: /*
444: PEPLinearExtract_Residual - Auxiliary routine that copies the solution of the
445: linear eigenproblem to the PEP object. The eigenvector of the generalized
446: problem is supposed to be
447: z = [ x ]
448: [ l*x ]
449: The eigenvector is taken from z(1:n) or z(n+1:2*n) depending on the explicitly
450: computed residual norm.
451: Finally, x is normalized so that ||x||_2 = 1.
452: */
453: static PetscErrorCode PEPLinearExtract_Residual(PEP pep,EPS eps)
454: {
455: PetscErrorCode ierr;
456: PetscInt i,k;
457: const PetscScalar *px;
458: PetscScalar *er=pep->eigr,*ei=pep->eigi;
459: PetscReal rn1,rn2;
460: Vec xr,xi=NULL,wr;
461: Mat A;
462: #if !defined(PETSC_USE_COMPLEX)
463: Vec wi;
464: const PetscScalar *py;
465: #endif
468: #if defined(PETSC_USE_COMPLEX)
469: PEPSetWorkVecs(pep,2);
470: #else
471: PEPSetWorkVecs(pep,4);
472: #endif
473: EPSGetOperators(eps,&A,NULL);
474: MatCreateVecs(A,&xr,NULL);
475: MatCreateVecsEmpty(pep->A[0],&wr,NULL);
476: #if !defined(PETSC_USE_COMPLEX)
477: VecDuplicate(xr,&xi);
478: VecDuplicateEmpty(wr,&wi);
479: #endif
480: for (i=0;i<pep->nconv;i++) {
481: EPSGetEigenpair(eps,i,NULL,NULL,xr,xi);
482: #if !defined(PETSC_USE_COMPLEX)
483: if (ei[i]!=0.0) { /* complex conjugate pair */
484: VecGetArrayRead(xr,&px);
485: VecGetArrayRead(xi,&py);
486: VecPlaceArray(wr,px);
487: VecPlaceArray(wi,py);
488: VecNormalizeComplex(wr,wi,PETSC_TRUE,NULL);
489: PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,wi,pep->work,&rn1);
490: BVInsertVec(pep->V,i,wr);
491: BVInsertVec(pep->V,i+1,wi);
492: for (k=1;k<pep->nmat-1;k++) {
493: VecResetArray(wr);
494: VecResetArray(wi);
495: VecPlaceArray(wr,px+k*pep->nloc);
496: VecPlaceArray(wi,py+k*pep->nloc);
497: VecNormalizeComplex(wr,wi,PETSC_TRUE,NULL);
498: PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,wi,pep->work,&rn2);
499: if (rn1>rn2) {
500: BVInsertVec(pep->V,i,wr);
501: BVInsertVec(pep->V,i+1,wi);
502: rn1 = rn2;
503: }
504: }
505: VecResetArray(wr);
506: VecResetArray(wi);
507: VecRestoreArrayRead(xr,&px);
508: VecRestoreArrayRead(xi,&py);
509: i++;
510: } else /* real eigenvalue */
511: #endif
512: {
513: VecGetArrayRead(xr,&px);
514: VecPlaceArray(wr,px);
515: VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL);
516: PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn1);
517: BVInsertVec(pep->V,i,wr);
518: for (k=1;k<pep->nmat-1;k++) {
519: VecResetArray(wr);
520: VecPlaceArray(wr,px+k*pep->nloc);
521: VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL);
522: PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn2);
523: if (rn1>rn2) {
524: BVInsertVec(pep->V,i,wr);
525: rn1 = rn2;
526: }
527: }
528: VecResetArray(wr);
529: VecRestoreArrayRead(xr,&px);
530: }
531: }
532: VecDestroy(&wr);
533: VecDestroy(&xr);
534: #if !defined(PETSC_USE_COMPLEX)
535: VecDestroy(&wi);
536: VecDestroy(&xi);
537: #endif
538: return(0);
539: }
541: /*
542: PEPLinearExtract_None - Same as PEPLinearExtract_Norm but always takes
543: the first block.
544: */
545: static PetscErrorCode PEPLinearExtract_None(PEP pep,EPS eps)
546: {
547: PetscErrorCode ierr;
548: PetscInt i;
549: const PetscScalar *px;
550: Mat A;
551: Vec xr,xi=NULL,w;
552: #if !defined(PETSC_USE_COMPLEX)
553: PetscScalar *ei=pep->eigi;
554: #endif
557: EPSGetOperators(eps,&A,NULL);
558: MatCreateVecs(A,&xr,NULL);
559: #if !defined(PETSC_USE_COMPLEX)
560: VecDuplicate(xr,&xi);
561: #endif
562: MatCreateVecsEmpty(pep->A[0],&w,NULL);
563: for (i=0;i<pep->nconv;i++) {
564: EPSGetEigenvector(eps,i,xr,xi);
565: #if !defined(PETSC_USE_COMPLEX)
566: if (ei[i]!=0.0) { /* complex conjugate pair */
567: VecGetArrayRead(xr,&px);
568: VecPlaceArray(w,px);
569: BVInsertVec(pep->V,i,w);
570: VecResetArray(w);
571: VecRestoreArrayRead(xr,&px);
572: VecGetArrayRead(xi,&px);
573: VecPlaceArray(w,px);
574: BVInsertVec(pep->V,i+1,w);
575: VecResetArray(w);
576: VecRestoreArrayRead(xi,&px);
577: i++;
578: } else /* real eigenvalue */
579: #endif
580: {
581: VecGetArrayRead(xr,&px);
582: VecPlaceArray(w,px);
583: BVInsertVec(pep->V,i,w);
584: VecResetArray(w);
585: VecRestoreArrayRead(xr,&px);
586: }
587: }
588: VecDestroy(&w);
589: VecDestroy(&xr);
590: #if !defined(PETSC_USE_COMPLEX)
591: VecDestroy(&xi);
592: #endif
593: return(0);
594: }
596: /*
597: PEPLinearExtract_Norm - Auxiliary routine that copies the solution of the
598: linear eigenproblem to the PEP object. The eigenvector of the generalized
599: problem is supposed to be
600: z = [ x ]
601: [ l*x ]
602: If |l|<1.0, the eigenvector is taken from z(1:n), otherwise from z(n+1:2*n).
603: Finally, x is normalized so that ||x||_2 = 1.
604: */
605: static PetscErrorCode PEPLinearExtract_Norm(PEP pep,EPS eps)
606: {
607: PetscErrorCode ierr;
608: PetscInt i,offset;
609: const PetscScalar *px;
610: PetscScalar *er=pep->eigr;
611: Mat A;
612: Vec xr,xi=NULL,w;
613: #if !defined(PETSC_USE_COMPLEX)
614: PetscScalar *ei=pep->eigi;
615: #endif
618: EPSGetOperators(eps,&A,NULL);
619: MatCreateVecs(A,&xr,NULL);
620: #if !defined(PETSC_USE_COMPLEX)
621: VecDuplicate(xr,&xi);
622: #endif
623: MatCreateVecsEmpty(pep->A[0],&w,NULL);
624: for (i=0;i<pep->nconv;i++) {
625: EPSGetEigenpair(eps,i,NULL,NULL,xr,xi);
626: if (SlepcAbsEigenvalue(er[i],ei[i])>1.0) offset = (pep->nmat-2)*pep->nloc;
627: else offset = 0;
628: #if !defined(PETSC_USE_COMPLEX)
629: if (ei[i]!=0.0) { /* complex conjugate pair */
630: VecGetArrayRead(xr,&px);
631: VecPlaceArray(w,px+offset);
632: BVInsertVec(pep->V,i,w);
633: VecResetArray(w);
634: VecRestoreArrayRead(xr,&px);
635: VecGetArrayRead(xi,&px);
636: VecPlaceArray(w,px+offset);
637: BVInsertVec(pep->V,i+1,w);
638: VecResetArray(w);
639: VecRestoreArrayRead(xi,&px);
640: i++;
641: } else /* real eigenvalue */
642: #endif
643: {
644: VecGetArrayRead(xr,&px);
645: VecPlaceArray(w,px+offset);
646: BVInsertVec(pep->V,i,w);
647: VecResetArray(w);
648: VecRestoreArrayRead(xr,&px);
649: }
650: }
651: VecDestroy(&w);
652: VecDestroy(&xr);
653: #if !defined(PETSC_USE_COMPLEX)
654: VecDestroy(&xi);
655: #endif
656: return(0);
657: }
659: PetscErrorCode PEPExtractVectors_Linear(PEP pep)
660: {
662: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
665: switch (pep->extract) {
666: case PEP_EXTRACT_NONE:
667: PEPLinearExtract_None(pep,ctx->eps);
668: break;
669: case PEP_EXTRACT_NORM:
670: PEPLinearExtract_Norm(pep,ctx->eps);
671: break;
672: case PEP_EXTRACT_RESIDUAL:
673: PEPLinearExtract_Residual(pep,ctx->eps);
674: break;
675: case PEP_EXTRACT_STRUCTURED:
676: SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Extraction not implemented in this solver");
677: }
678: return(0);
679: }
681: PetscErrorCode PEPSolve_Linear(PEP pep)
682: {
684: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
685: PetscScalar sigma;
686: PetscBool flg;
687: PetscInt i;
690: EPSSolve(ctx->eps);
691: EPSGetConverged(ctx->eps,&pep->nconv);
692: EPSGetIterationNumber(ctx->eps,&pep->its);
693: EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&pep->reason);
695: /* recover eigenvalues */
696: for (i=0;i<pep->nconv;i++) {
697: EPSGetEigenpair(ctx->eps,i,&pep->eigr[i],&pep->eigi[i],NULL,NULL);
698: pep->eigr[i] *= pep->sfactor;
699: pep->eigi[i] *= pep->sfactor;
700: }
702: /* restore target */
703: EPSGetTarget(ctx->eps,&sigma);
704: EPSSetTarget(ctx->eps,sigma*pep->sfactor);
706: STGetTransform(pep->st,&flg);
707: if (flg && pep->ops->backtransform) {
708: (*pep->ops->backtransform)(pep);
709: }
710: if (pep->sfactor!=1.0) {
711: /* Restore original values */
712: for (i=0;i<pep->nmat;i++){
713: pep->pbc[pep->nmat+i] *= pep->sfactor;
714: pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
715: }
716: if (!flg && !ctx->explicitmatrix) {
717: STScaleShift(pep->st,pep->sfactor);
718: }
719: }
720: if (ctx->explicitmatrix || !flg) {
721: RGPopScale(pep->rg);
722: }
723: return(0);
724: }
726: static PetscErrorCode EPSMonitor_Linear(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
727: {
728: PEP pep = (PEP)ctx;
732: PEPMonitor(pep,its,nconv,eigr,eigi,errest,nest);
733: return(0);
734: }
736: PetscErrorCode PEPSetFromOptions_Linear(PetscOptionItems *PetscOptionsObject,PEP pep)
737: {
739: PetscBool set,val;
740: PetscInt k;
741: PetscReal array[2]={0,0};
742: PetscBool flg;
743: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
746: PetscOptionsHead(PetscOptionsObject,"PEP Linear Options");
748: k = 2;
749: PetscOptionsRealArray("-pep_linear_linearization","Parameters of the linearization","PEPLinearSetLinearization",array,&k,&flg);
750: if (flg) {
751: PEPLinearSetLinearization(pep,array[0],array[1]);
752: }
754: PetscOptionsBool("-pep_linear_explicitmatrix","Use explicit matrix in linearization","PEPLinearSetExplicitMatrix",ctx->explicitmatrix,&val,&set);
755: if (set) { PEPLinearSetExplicitMatrix(pep,val); }
757: PetscOptionsTail();
759: if (!ctx->eps) { PEPLinearGetEPS(pep,&ctx->eps); }
760: EPSSetFromOptions(ctx->eps);
761: return(0);
762: }
764: static PetscErrorCode PEPLinearSetLinearization_Linear(PEP pep,PetscReal alpha,PetscReal beta)
765: {
766: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
769: if (beta==0.0 && alpha==0.0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
770: ctx->alpha = alpha;
771: ctx->beta = beta;
772: return(0);
773: }
775: /*@
776: PEPLinearSetLinearization - Set the coefficients that define
777: the linearization of a quadratic eigenproblem.
779: Logically Collective on PEP
781: Input Parameters:
782: + pep - polynomial eigenvalue solver
783: . alpha - first parameter of the linearization
784: - beta - second parameter of the linearization
786: Options Database Key:
787: . -pep_linear_linearization <alpha,beta> - Sets the coefficients
789: Notes:
790: Cannot pass zero for both alpha and beta. The default values are
791: alpha=1 and beta=0.
793: Level: advanced
795: .seealso: PEPLinearGetLinearization()
796: @*/
797: PetscErrorCode PEPLinearSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
798: {
805: PetscTryMethod(pep,"PEPLinearSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
806: return(0);
807: }
809: static PetscErrorCode PEPLinearGetLinearization_Linear(PEP pep,PetscReal *alpha,PetscReal *beta)
810: {
811: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
814: if (alpha) *alpha = ctx->alpha;
815: if (beta) *beta = ctx->beta;
816: return(0);
817: }
819: /*@
820: PEPLinearGetLinearization - Returns the coefficients that define
821: the linearization of a quadratic eigenproblem.
823: Not Collective
825: Input Parameter:
826: . pep - polynomial eigenvalue solver
828: Output Parameters:
829: + alpha - the first parameter of the linearization
830: - beta - the second parameter of the linearization
832: Level: advanced
834: .seealso: PEPLinearSetLinearization()
835: @*/
836: PetscErrorCode PEPLinearGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
837: {
842: PetscUseMethod(pep,"PEPLinearGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
843: return(0);
844: }
846: static PetscErrorCode PEPLinearSetExplicitMatrix_Linear(PEP pep,PetscBool explicitmatrix)
847: {
848: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
851: if (ctx->explicitmatrix != explicitmatrix) {
852: ctx->explicitmatrix = explicitmatrix;
853: pep->state = PEP_STATE_INITIAL;
854: }
855: return(0);
856: }
858: /*@
859: PEPLinearSetExplicitMatrix - Indicate if the matrices A and B for the
860: linearization of the problem must be built explicitly.
862: Logically Collective on PEP
864: Input Parameters:
865: + pep - polynomial eigenvalue solver
866: - explicit - boolean flag indicating if the matrices are built explicitly
868: Options Database Key:
869: . -pep_linear_explicitmatrix <boolean> - Indicates the boolean flag
871: Level: advanced
873: .seealso: PEPLinearGetExplicitMatrix()
874: @*/
875: PetscErrorCode PEPLinearSetExplicitMatrix(PEP pep,PetscBool explicitmatrix)
876: {
882: PetscTryMethod(pep,"PEPLinearSetExplicitMatrix_C",(PEP,PetscBool),(pep,explicitmatrix));
883: return(0);
884: }
886: static PetscErrorCode PEPLinearGetExplicitMatrix_Linear(PEP pep,PetscBool *explicitmatrix)
887: {
888: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
891: *explicitmatrix = ctx->explicitmatrix;
892: return(0);
893: }
895: /*@
896: PEPLinearGetExplicitMatrix - Returns the flag indicating if the matrices
897: A and B for the linearization are built explicitly.
899: Not Collective
901: Input Parameter:
902: . pep - polynomial eigenvalue solver
904: Output Parameter:
905: . explicitmatrix - the mode flag
907: Level: advanced
909: .seealso: PEPLinearSetExplicitMatrix()
910: @*/
911: PetscErrorCode PEPLinearGetExplicitMatrix(PEP pep,PetscBool *explicitmatrix)
912: {
918: PetscUseMethod(pep,"PEPLinearGetExplicitMatrix_C",(PEP,PetscBool*),(pep,explicitmatrix));
919: return(0);
920: }
922: static PetscErrorCode PEPLinearSetEPS_Linear(PEP pep,EPS eps)
923: {
925: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
928: PetscObjectReference((PetscObject)eps);
929: EPSDestroy(&ctx->eps);
930: ctx->eps = eps;
931: ctx->usereps = PETSC_TRUE;
932: PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->eps);
933: pep->state = PEP_STATE_INITIAL;
934: return(0);
935: }
937: /*@
938: PEPLinearSetEPS - Associate an eigensolver object (EPS) to the
939: polynomial eigenvalue solver.
941: Collective on PEP
943: Input Parameters:
944: + pep - polynomial eigenvalue solver
945: - eps - the eigensolver object
947: Level: advanced
949: .seealso: PEPLinearGetEPS()
950: @*/
951: PetscErrorCode PEPLinearSetEPS(PEP pep,EPS eps)
952: {
959: PetscTryMethod(pep,"PEPLinearSetEPS_C",(PEP,EPS),(pep,eps));
960: return(0);
961: }
963: static PetscErrorCode PEPLinearGetEPS_Linear(PEP pep,EPS *eps)
964: {
966: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
969: if (!ctx->eps) {
970: EPSCreate(PetscObjectComm((PetscObject)pep),&ctx->eps);
971: PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)pep,1);
972: EPSSetOptionsPrefix(ctx->eps,((PetscObject)pep)->prefix);
973: EPSAppendOptionsPrefix(ctx->eps,"pep_linear_");
974: PetscLogObjectParent((PetscObject)pep,(PetscObject)ctx->eps);
975: PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)pep)->options);
976: EPSMonitorSet(ctx->eps,EPSMonitor_Linear,pep,NULL);
977: }
978: *eps = ctx->eps;
979: return(0);
980: }
982: /*@
983: PEPLinearGetEPS - Retrieve the eigensolver object (EPS) associated
984: to the polynomial eigenvalue solver.
986: Not Collective
988: Input Parameter:
989: . pep - polynomial eigenvalue solver
991: Output Parameter:
992: . eps - the eigensolver object
994: Level: advanced
996: .seealso: PEPLinearSetEPS()
997: @*/
998: PetscErrorCode PEPLinearGetEPS(PEP pep,EPS *eps)
999: {
1005: PetscUseMethod(pep,"PEPLinearGetEPS_C",(PEP,EPS*),(pep,eps));
1006: return(0);
1007: }
1009: PetscErrorCode PEPView_Linear(PEP pep,PetscViewer viewer)
1010: {
1012: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
1013: PetscBool isascii;
1016: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1017: if (isascii) {
1018: if (!ctx->eps) { PEPLinearGetEPS(pep,&ctx->eps); }
1019: PetscViewerASCIIPrintf(viewer," %s matrices\n",ctx->explicitmatrix? "explicit": "implicit");
1020: PetscViewerASCIIPrintf(viewer," linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta);
1021: PetscViewerASCIIPushTab(viewer);
1022: EPSView(ctx->eps,viewer);
1023: PetscViewerASCIIPopTab(viewer);
1024: }
1025: return(0);
1026: }
1028: PetscErrorCode PEPReset_Linear(PEP pep)
1029: {
1031: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
1034: if (!ctx->eps) { EPSReset(ctx->eps); }
1035: MatDestroy(&ctx->A);
1036: MatDestroy(&ctx->B);
1037: VecDestroy(&ctx->w[0]);
1038: VecDestroy(&ctx->w[1]);
1039: VecDestroy(&ctx->w[2]);
1040: VecDestroy(&ctx->w[3]);
1041: VecDestroy(&ctx->w[4]);
1042: VecDestroy(&ctx->w[5]);
1043: return(0);
1044: }
1046: PetscErrorCode PEPDestroy_Linear(PEP pep)
1047: {
1049: PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
1052: EPSDestroy(&ctx->eps);
1053: PetscFree(pep->data);
1054: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",NULL);
1055: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",NULL);
1056: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",NULL);
1057: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",NULL);
1058: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",NULL);
1059: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",NULL);
1060: return(0);
1061: }
1063: SLEPC_EXTERN PetscErrorCode PEPCreate_Linear(PEP pep)
1064: {
1066: PEP_LINEAR *ctx;
1069: PetscNewLog(pep,&ctx);
1070: ctx->explicitmatrix = PETSC_FALSE;
1071: pep->data = (void*)ctx;
1072: ctx->alpha = 1.0;
1073: ctx->beta = 0.0;
1075: pep->ops->solve = PEPSolve_Linear;
1076: pep->ops->setup = PEPSetUp_Linear;
1077: pep->ops->setfromoptions = PEPSetFromOptions_Linear;
1078: pep->ops->destroy = PEPDestroy_Linear;
1079: pep->ops->reset = PEPReset_Linear;
1080: pep->ops->view = PEPView_Linear;
1081: pep->ops->backtransform = PEPBackTransform_Default;
1082: pep->ops->computevectors = PEPComputeVectors_Default;
1083: pep->ops->extractvectors = PEPExtractVectors_Linear;
1085: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",PEPLinearSetLinearization_Linear);
1086: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",PEPLinearGetLinearization_Linear);
1087: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",PEPLinearSetEPS_Linear);
1088: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",PEPLinearGetEPS_Linear);
1089: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",PEPLinearSetExplicitMatrix_Linear);
1090: PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",PEPLinearGetExplicitMatrix_Linear);
1091: return(0);
1092: }