Actual source code: nleigs-fullb.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: Full basis for the linearization of the rational approximation of non-linear eigenproblems
12: */
14: #include <slepc/private/nepimpl.h> /*I "slepcnep.h" I*/
15: #include "nleigs.h"
17: static PetscErrorCode MatMult_FullBasis_Sinvert(Mat M,Vec x,Vec y)
18: {
19: PetscErrorCode ierr;
20: NEP_NLEIGS *ctx;
21: NEP nep;
22: const PetscScalar *px;
23: PetscScalar *beta,*s,*xi,*t,*py,sigma;
24: PetscInt nmat,d,i,k,m;
25: Vec xx,xxx,yy,yyy,w,ww,www;
28: MatShellGetContext(M,(void**)&nep);
29: ctx = (NEP_NLEIGS*)nep->data;
30: beta = ctx->beta; s = ctx->s; xi = ctx->xi;
31: sigma = ctx->shifts[0];
32: nmat = ctx->nmat;
33: d = nmat-1;
34: m = nep->nloc;
35: PetscMalloc1(ctx->nmat,&t);
36: xx = ctx->w[0]; xxx = ctx->w[1]; yy = ctx->w[2]; yyy=ctx->w[3];
37: w = nep->work[0]; ww = nep->work[1]; www = nep->work[2];
38: VecGetArrayRead(x,&px);
39: VecGetArray(y,&py);
40: VecPlaceArray(xx,px+(d-1)*m);
41: VecPlaceArray(xxx,px+(d-2)*m);
42: VecPlaceArray(yy,py+(d-2)*m);
43: VecCopy(xxx,yy);
44: VecAXPY(yy,beta[d-1]/xi[d-2],xx);
45: VecScale(yy,1.0/(s[d-2]-sigma));
46: VecResetArray(xx);
47: VecResetArray(xxx);
48: VecResetArray(yy);
49: for (i=d-3;i>=0;i--) {
50: VecPlaceArray(xx,px+(i+1)*m);
51: VecPlaceArray(xxx,px+i*m);
52: VecPlaceArray(yy,py+i*m);
53: VecPlaceArray(yyy,py+(i+1)*m);
54: VecCopy(xxx,yy);
55: VecAXPY(yy,beta[i+1]/xi[i],xx);
56: VecAXPY(yy,-beta[i+1]*(1.0-sigma/xi[i]),yyy);
57: VecScale(yy,1.0/(s[i]-sigma));
58: VecResetArray(xx);
59: VecResetArray(xxx);
60: VecResetArray(yy);
61: VecResetArray(yyy);
62: }
63: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
64: VecZeroEntries(w);
65: for (k=0;k<nep->nt;k++) {
66: VecZeroEntries(ww);
67: VecPlaceArray(xx,px+(d-1)*m);
68: VecAXPY(ww,-ctx->coeffD[k+nep->nt*d]/beta[d],xx);
69: VecResetArray(xx);
70: for (i=0;i<d-1;i++) {
71: VecPlaceArray(yy,py+i*m);
72: VecAXPY(ww,-ctx->coeffD[nep->nt*i+k],yy);
73: VecResetArray(yy);
74: }
75: MatMult(nep->A[k],ww,www);
76: VecAXPY(w,1.0,www);
77: }
78: } else {
79: VecPlaceArray(xx,px+(d-1)*m);
80: MatMult(ctx->D[d],xx,w);
81: VecScale(w,-1.0/beta[d]);
82: VecResetArray(xx);
83: for (i=0;i<d-1;i++) {
84: VecPlaceArray(yy,py+i*m);
85: MatMult(ctx->D[i],yy,ww);
86: VecResetArray(yy);
87: VecAXPY(w,-1.0,ww);
88: }
89: }
90: VecPlaceArray(yy,py+(d-1)*m);
91: KSPSolve(ctx->ksp[0],w,yy);
92: NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t);
93: for (i=0;i<d-1;i++) {
94: VecPlaceArray(yyy,py+i*m);
95: VecAXPY(yyy,t[i],yy);
96: VecResetArray(yyy);
97: }
98: VecScale(yy,t[d-1]);
99: VecResetArray(yy);
100: VecRestoreArrayRead(x,&px);
101: VecRestoreArray(y,&py);
102: PetscFree(t);
103: return(0);
104: }
106: static PetscErrorCode MatMultTranspose_FullBasis_Sinvert(Mat M,Vec x,Vec y)
107: {
108: PetscErrorCode ierr;
109: NEP_NLEIGS *ctx;
110: NEP nep;
111: const PetscScalar *px;
112: PetscScalar *beta,*s,*xi,*t,*py,sigma;
113: PetscInt nmat,d,i,k,m;
114: Vec xx,yy,yyy,w,z0;
117: MatShellGetContext(M,(void**)&nep);
118: ctx = (NEP_NLEIGS*)nep->data;
119: beta = ctx->beta; s = ctx->s; xi = ctx->xi;
120: sigma = ctx->shifts[0];
121: nmat = ctx->nmat;
122: d = nmat-1;
123: m = nep->nloc;
124: PetscMalloc1(ctx->nmat,&t);
125: xx = ctx->w[0]; yy = ctx->w[1]; yyy=ctx->w[2];
126: w = nep->work[0]; z0 = nep->work[1];
127: VecGetArrayRead(x,&px);
128: VecGetArray(y,&py);
129: NEPNLEIGSEvalNRTFunct(nep,d,sigma,t);
130: VecPlaceArray(xx,px+(d-1)*m);
131: VecCopy(xx,w);
132: VecScale(w,t[d-1]);
133: VecResetArray(xx);
134: for (i=0;i<d-1;i++) {
135: VecPlaceArray(xx,px+i*m);
136: VecAXPY(w,t[i],xx);
137: VecResetArray(xx);
138: }
139: KSPSolveTranspose(ctx->ksp[0],w,z0);
141: VecPlaceArray(yy,py);
142: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
143: VecZeroEntries(yy);
144: for (k=0;k<nep->nt;k++) {
145: MatMult(nep->A[k],z0,w);
146: VecAXPY(yy,ctx->coeffD[k],w);
147: }
148: } else {
149: MatMultTranspose(ctx->D[0],z0,yy);
150: }
151: VecPlaceArray(xx,px);
152: VecAXPY(yy,-1.0,xx);;
153: VecResetArray(xx);
154: VecScale(yy,-1.0/(s[0]-sigma));
155: VecResetArray(yy);
156: for (i=2;i<d;i++) {
157: VecPlaceArray(yy,py+(i-1)*m);
158: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
159: VecZeroEntries(yy);
160: for (k=0;k<nep->nt;k++) {
161: MatMult(nep->A[k],z0,w);
162: VecAXPY(yy,ctx->coeffD[k+(i-1)*nep->nt],w);
163: }
164: } else {
165: MatMultTranspose(ctx->D[i-1],z0,yy);
166: }
167: VecPlaceArray(yyy,py+(i-2)*m);
168: VecAXPY(yy,beta[i-1]*(1.0-sigma/xi[i-2]),yyy);
169: VecResetArray(yyy);
170: VecPlaceArray(xx,px+(i-1)*m);
171: VecAXPY(yy,-1.0,xx);
172: VecResetArray(xx);
173: VecScale(yy,-1.0/(s[i-1]-sigma));
174: VecResetArray(yy);
175: }
176: VecPlaceArray(yy,py+(d-1)*m);
177: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
178: VecZeroEntries(yy);
179: for (k=0;k<nep->nt;k++) {
180: MatMult(nep->A[k],z0,w);
181: VecAXPY(yy,ctx->coeffD[k+d*nep->nt],w);
182: }
183: } else {
184: MatMultTranspose(ctx->D[d],z0,yy);
185: }
186: VecScale(yy,-1.0/beta[d]);
187: VecPlaceArray(yyy,py+(d-2)*m);
188: VecAXPY(yy,beta[d-1]/xi[d-2],yyy);
189: VecResetArray(yyy);
190: VecResetArray(yy);
192: for (i=d-2;i>0;i--) {
193: VecPlaceArray(yyy,py+(i-1)*m);
194: VecPlaceArray(yy,py+i*m);
195: VecAXPY(yy,beta[i]/xi[i-1],yyy);
196: VecResetArray(yyy);
197: VecResetArray(yy);
198: }
200: VecRestoreArrayRead(x,&px);
201: VecRestoreArray(y,&py);
202: PetscFree(t);
203: return(0);
204: }
206: static PetscErrorCode BackTransform_FullBasis(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
207: {
209: NEP nep;
212: STShellGetContext(st,(void**)&nep);
213: NEPNLEIGSBackTransform((PetscObject)nep,n,eigr,eigi);
214: return(0);
215: }
217: static PetscErrorCode Apply_FullBasis(ST st,Vec x,Vec y)
218: {
220: NEP nep;
221: NEP_NLEIGS *ctx;
224: STShellGetContext(st,(void**)&nep);
225: ctx = (NEP_NLEIGS*)nep->data;
226: MatMult(ctx->A,x,y);
227: return(0);
228: }
230: static PetscErrorCode ApplyTranspose_FullBasis(ST st,Vec x,Vec y)
231: {
233: NEP nep;
234: NEP_NLEIGS *ctx;
237: STShellGetContext(st,(void**)&nep);
238: ctx = (NEP_NLEIGS*)nep->data;
239: MatMultTranspose(ctx->A,x,y);
240: return(0);
241: }
243: PetscErrorCode NEPSetUp_NLEIGS_FullBasis(NEP nep)
244: {
246: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
247: ST st;
248: Mat Q;
249: PetscInt i=0,deg=ctx->nmat-1;
250: PetscBool trackall,istrivial,ks;
251: PetscScalar *epsarray,*neparray;
252: Vec veps,w=NULL;
253: EPSWhich which;
256: if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
257: EPSGetST(ctx->eps,&st);
258: EPSSetTarget(ctx->eps,nep->target);
259: STSetDefaultShift(st,nep->target);
260: if (!((PetscObject)(ctx->eps))->type_name) {
261: EPSSetType(ctx->eps,EPSKRYLOVSCHUR);
262: } else {
263: PetscObjectTypeCompare((PetscObject)ctx->eps,EPSKRYLOVSCHUR,&ks);
264: if (!ks) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Full-basis option only implemented for Krylov-Schur");
265: }
266: STSetType(st,STSHELL);
267: STShellSetContext(st,(PetscObject)nep);
268: STShellSetBackTransform(st,BackTransform_FullBasis);
269: KSPGetOperators(ctx->ksp[0],&Q,NULL);
270: MatCreateVecsEmpty(Q,&ctx->w[0],&ctx->w[1]);
271: MatCreateVecsEmpty(Q,&ctx->w[2],&ctx->w[3]);
272: PetscLogObjectParents(nep,6,ctx->w);
273: MatCreateShell(PetscObjectComm((PetscObject)nep),deg*nep->nloc,deg*nep->nloc,deg*nep->n,deg*nep->n,nep,&ctx->A);
274: MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_FullBasis_Sinvert);
275: MatShellSetOperation(ctx->A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_FullBasis_Sinvert);
276: STShellSetApply(st,Apply_FullBasis);
277: STShellSetApplyTranspose(st,ApplyTranspose_FullBasis);
278: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->A);
279: EPSSetOperators(ctx->eps,ctx->A,NULL);
280: EPSSetProblemType(ctx->eps,EPS_NHEP);
281: switch (nep->which) {
282: case NEP_TARGET_MAGNITUDE: which = EPS_TARGET_MAGNITUDE; break;
283: case NEP_TARGET_REAL: which = EPS_TARGET_REAL; break;
284: case NEP_TARGET_IMAGINARY: which = EPS_TARGET_IMAGINARY; break;
285: case NEP_WHICH_USER: which = EPS_WHICH_USER;
286: EPSSetEigenvalueComparison(ctx->eps,nep->sc->comparison,nep->sc->comparisonctx);
287: break;
288: default: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Should set a target selection in NEPSetWhichEigenvalues()");
289: }
290: EPSSetWhichEigenpairs(ctx->eps,which);
291: RGIsTrivial(nep->rg,&istrivial);
292: if (!istrivial) { EPSSetRG(ctx->eps,nep->rg);}
293: EPSSetDimensions(ctx->eps,nep->nev,nep->ncv?nep->ncv:PETSC_DEFAULT,nep->mpd?nep->mpd:PETSC_DEFAULT);
294: EPSSetTolerances(ctx->eps,nep->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL:nep->tol,nep->max_it?nep->max_it:PETSC_DEFAULT);
295: EPSSetTwoSided(ctx->eps,nep->twosided);
296: /* Transfer the trackall option from pep to eps */
297: NEPGetTrackAll(nep,&trackall);
298: EPSSetTrackAll(ctx->eps,trackall);
300: /* process initial vector */
301: if (nep->nini<0) {
302: VecCreateMPI(PetscObjectComm((PetscObject)ctx->eps),deg*nep->nloc,deg*nep->n,&veps);
303: VecGetArray(veps,&epsarray);
304: for (i=0;i<deg;i++) {
305: if (i<-nep->nini) {
306: VecGetArray(nep->IS[i],&neparray);
307: PetscMemcpy(epsarray+i*nep->nloc,neparray,nep->nloc*sizeof(PetscScalar));
308: VecRestoreArray(nep->IS[i],&neparray);
309: } else {
310: if (!w) { VecDuplicate(nep->IS[0],&w); }
311: VecSetRandom(w,NULL);
312: VecGetArray(w,&neparray);
313: PetscMemcpy(epsarray+i*nep->nloc,neparray,nep->nloc*sizeof(PetscScalar));
314: VecRestoreArray(w,&neparray);
315: }
316: }
317: VecRestoreArray(veps,&epsarray);
318: EPSSetInitialSpace(ctx->eps,1,&veps);
319: VecDestroy(&veps);
320: VecDestroy(&w);
321: SlepcBasisDestroy_Private(&nep->nini,&nep->IS);
322: }
324: EPSSetUp(ctx->eps);
325: EPSGetDimensions(ctx->eps,NULL,&nep->ncv,&nep->mpd);
326: EPSGetTolerances(ctx->eps,NULL,&nep->max_it);
327: NEPAllocateSolution(nep,0);
328: return(0);
329: }
331: /*
332: NEPNLEIGSExtract_None - Extracts the first block of the basis
333: and normalizes the columns.
334: */
335: static PetscErrorCode NEPNLEIGSExtract_None(NEP nep,EPS eps)
336: {
337: PetscErrorCode ierr;
338: PetscInt i,k,m,d;
339: const PetscScalar *px;
340: PetscScalar sigma=nep->target,*b;
341: Mat A;
342: Vec xxr,xxi=NULL,w,t,xx;
343: PetscReal norm;
344: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
347: d = ctx->nmat-1;
348: EPSGetOperators(eps,&A,NULL);
349: MatCreateVecs(A,&xxr,NULL);
350: #if !defined(PETSC_USE_COMPLEX)
351: VecDuplicate(xxr,&xxi);
352: #endif
353: w = nep->work[0];
354: for (i=0;i<nep->nconv;i++) {
355: EPSGetEigenvector(eps,i,xxr,xxi);
356: VecGetArrayRead(xxr,&px);
357: VecPlaceArray(w,px);
358: BVInsertVec(nep->V,i,w);
359: BVNormColumn(nep->V,i,NORM_2,&norm);
360: BVScaleColumn(nep->V,i,1.0/norm);
361: VecResetArray(w);
362: VecRestoreArrayRead(xxr,&px);
363: }
364: if (nep->twosided) {
365: PetscMalloc1(ctx->nmat,&b);
366: NEPNLEIGSEvalNRTFunct(nep,d,sigma,b);
367: m = nep->nloc;
368: xx = ctx->w[0];
369: w = nep->work[0]; t = nep->work[1];
370: for (k=0;k<nep->nconv;k++) {
371: EPSGetLeftEigenvector(eps,k,xxr,xxi);
372: VecGetArrayRead(xxr,&px);
373: VecPlaceArray(xx,px+(d-1)*m);
374: VecCopy(xx,w);
375: VecScale(w,PetscConj(b[d-1]));
376: VecResetArray(xx);
377: for (i=0;i<d-1;i++) {
378: VecPlaceArray(xx,px+i*m);
379: VecAXPY(w,PetscConj(b[i]),xx);
380: VecResetArray(xx);
381: }
382: VecConjugate(w);
383: KSPSolveTranspose(ctx->ksp[0],w,t);
384: VecConjugate(t);
385: BVInsertVec(nep->W,k,t);
386: BVNormColumn(nep->W,k,NORM_2,&norm);
387: BVScaleColumn(nep->W,k,1.0/norm);
388: VecRestoreArrayRead(xxr,&px);
389: }
390: PetscFree(b);
391: }
392: VecDestroy(&xxr);
393: #if !defined(PETSC_USE_COMPLEX)
394: VecDestroy(&xxi);
395: #endif
396: return(0);
397: }
399: PetscErrorCode NEPSolve_NLEIGS_FullBasis(NEP nep)
400: {
402: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
403: PetscInt i;
404: PetscScalar eigi=0.0;
407: EPSSolve(ctx->eps);
408: EPSGetConverged(ctx->eps,&nep->nconv);
409: EPSGetIterationNumber(ctx->eps,&nep->its);
410: EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&nep->reason);
412: /* recover eigenvalues */
413: for (i=0;i<nep->nconv;i++) {
414: EPSGetEigenpair(ctx->eps,i,&nep->eigr[i],&eigi,NULL,NULL);
415: #if !defined(PETSC_USE_COMPLEX)
416: if (eigi!=0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Complex value requires complex arithmetic");
417: #endif
418: }
419: NEPNLEIGSExtract_None(nep,ctx->eps);
420: return(0);
421: }
423: PetscErrorCode NEPNLEIGSSetEPS_NLEIGS(NEP nep,EPS eps)
424: {
426: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
429: PetscObjectReference((PetscObject)eps);
430: EPSDestroy(&ctx->eps);
431: ctx->eps = eps;
432: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->eps);
433: nep->state = NEP_STATE_INITIAL;
434: return(0);
435: }
437: /*@
438: NEPNLEIGSSetEPS - Associate an eigensolver object (EPS) to the NLEIGS solver.
440: Collective on NEP
442: Input Parameters:
443: + nep - nonlinear eigenvalue solver
444: - eps - the eigensolver object
446: Level: advanced
448: .seealso: NEPNLEIGSGetEPS()
449: @*/
450: PetscErrorCode NEPNLEIGSSetEPS(NEP nep,EPS eps)
451: {
458: PetscTryMethod(nep,"NEPNLEIGSSetEPS_C",(NEP,EPS),(nep,eps));
459: return(0);
460: }
462: static PetscErrorCode EPSMonitor_NLEIGS(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
463: {
464: NEP nep = (NEP)ctx;
465: PetscInt i,nv = PetscMin(nest,nep->ncv);
469: for (i=0;i<nv;i++) {
470: nep->eigr[i] = eigr[i];
471: nep->eigi[i] = eigi[i];
472: nep->errest[i] = errest[i];
473: }
474: NEPNLEIGSBackTransform((PetscObject)nep,nv,nep->eigr,nep->eigi);
475: NEPMonitor(nep,its,nconv,nep->eigr,nep->eigi,nep->errest,nest);
476: return(0);
477: }
479: PetscErrorCode NEPNLEIGSGetEPS_NLEIGS(NEP nep,EPS *eps)
480: {
482: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
485: if (!ctx->eps) {
486: EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->eps);
487: PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)nep,1);
488: EPSSetOptionsPrefix(ctx->eps,((PetscObject)nep)->prefix);
489: EPSAppendOptionsPrefix(ctx->eps,"nep_nleigs_");
490: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->eps);
491: PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)nep)->options);
492: EPSMonitorSet(ctx->eps,EPSMonitor_NLEIGS,nep,NULL);
493: }
494: *eps = ctx->eps;
495: return(0);
496: }
498: /*@
499: NEPNLEIGSGetEPS - Retrieve the eigensolver object (EPS) associated
500: to the nonlinear eigenvalue solver.
502: Not Collective
504: Input Parameter:
505: . nep - nonlinear eigenvalue solver
507: Output Parameter:
508: . eps - the eigensolver object
510: Level: advanced
512: .seealso: NEPNLEIGSSetEPS()
513: @*/
514: PetscErrorCode NEPNLEIGSGetEPS(NEP nep,EPS *eps)
515: {
521: PetscUseMethod(nep,"NEPNLEIGSGetEPS_C",(NEP,EPS*),(nep,eps));
522: return(0);
523: }