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