Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 : /*
11 : SLEPc eigensolver: "ciss"
12 :
13 : Method: Contour Integral Spectral Slicing
14 :
15 : Algorithm:
16 :
17 : Contour integral based on Sakurai-Sugiura method to construct a
18 : subspace, with various eigenpair extractions (Rayleigh-Ritz,
19 : explicit moment).
20 :
21 : Based on code contributed by Y. Maeda, T. Sakurai.
22 :
23 : References:
24 :
25 : [1] J. Asakura, T. Sakurai, H. Tadano, T. Ikegami, K. Kimura, "A
26 : numerical method for nonlinear eigenvalue problems using contour
27 : integrals", JSIAM Lett. 1:52-55, 2009.
28 :
29 : [2] S. Yokota and T. Sakurai, "A projection method for nonlinear
30 : eigenvalue problems using contour integrals", JSIAM Lett.
31 : 5:41-44, 2013.
32 : */
33 :
34 : #include <slepc/private/nepimpl.h> /*I "slepcnep.h" I*/
35 : #include <slepc/private/slepccontour.h>
36 :
37 : typedef struct _n_nep_ciss_project *NEP_CISS_PROJECT;
38 : typedef struct {
39 : /* parameters */
40 : PetscInt N; /* number of integration points (32) */
41 : PetscInt L; /* block size (16) */
42 : PetscInt M; /* moment degree (N/4 = 4) */
43 : PetscReal delta; /* threshold of singular value (1e-12) */
44 : PetscInt L_max; /* maximum number of columns of the source matrix V */
45 : PetscReal spurious_threshold; /* discard spurious eigenpairs */
46 : PetscBool isreal; /* T(z) is real for real z */
47 : PetscInt npart; /* number of partitions */
48 : PetscInt refine_inner;
49 : PetscInt refine_blocksize;
50 : NEPCISSExtraction extraction;
51 : /* private data */
52 : SlepcContourData contour;
53 : PetscReal *sigma; /* threshold for numerical rank */
54 : PetscScalar *weight;
55 : PetscScalar *omega;
56 : PetscScalar *pp;
57 : BV V;
58 : BV S;
59 : BV Y;
60 : PetscBool useconj;
61 : Mat J; /* auxiliary matrix when using subcomm */
62 : BV pV;
63 : NEP_CISS_PROJECT dsctxf;
64 : PetscObjectId rgid;
65 : PetscObjectState rgstate;
66 : } NEP_CISS;
67 :
68 : struct _n_nep_ciss_project {
69 : NEP nep;
70 : BV Q;
71 : };
72 :
73 408 : static PetscErrorCode NEPContourDSComputeMatrix(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat,void *ctx)
74 : {
75 408 : NEP_CISS_PROJECT proj = (NEP_CISS_PROJECT)ctx;
76 408 : Mat M,fun;
77 :
78 408 : PetscFunctionBegin;
79 408 : if (!deriv) {
80 408 : PetscCall(NEPComputeFunction(proj->nep,lambda,proj->nep->function,proj->nep->function));
81 408 : fun = proj->nep->function;
82 : } else {
83 0 : PetscCall(NEPComputeJacobian(proj->nep,lambda,proj->nep->jacobian));
84 0 : fun = proj->nep->jacobian;
85 : }
86 408 : PetscCall(DSGetMat(ds,mat,&M));
87 408 : PetscCall(BVMatProject(proj->Q,fun,proj->Q,M));
88 408 : PetscCall(DSRestoreMat(ds,mat,&M));
89 408 : PetscFunctionReturn(PETSC_SUCCESS);
90 : }
91 :
92 224 : static PetscErrorCode NEPComputeFunctionSubcomm(NEP nep,PetscScalar lambda,Mat T,Mat P,PetscBool deriv)
93 : {
94 224 : PetscInt i;
95 224 : PetscScalar alpha;
96 224 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
97 :
98 224 : PetscFunctionBegin;
99 224 : PetscAssert(nep->fui!=NEP_USER_INTERFACE_CALLBACK,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Should not arrive here with callbacks");
100 224 : PetscCall(MatZeroEntries(T));
101 224 : if (!deriv && T != P) PetscCall(MatZeroEntries(P));
102 832 : for (i=0;i<nep->nt;i++) {
103 608 : if (!deriv) PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
104 304 : else PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
105 608 : PetscCall(MatAXPY(T,alpha,ctx->contour->pA[i],nep->mstr));
106 608 : if (!deriv && T != P) PetscCall(MatAXPY(P,alpha,ctx->contour->pP[i],nep->mstrp));
107 : }
108 224 : PetscFunctionReturn(PETSC_SUCCESS);
109 : }
110 :
111 : /*
112 : Set up KSP solvers for every integration point
113 : */
114 40 : static PetscErrorCode NEPCISSSetUp(NEP nep,Mat T,Mat P)
115 : {
116 40 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
117 40 : SlepcContourData contour;
118 40 : PetscInt i,p_id;
119 40 : Mat Amat,Pmat;
120 :
121 40 : PetscFunctionBegin;
122 40 : if (!ctx->contour || !ctx->contour->ksp) PetscCall(NEPCISSGetKSPs(nep,NULL,NULL));
123 40 : contour = ctx->contour;
124 40 : PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Something went wrong with NEPCISSGetKSPs()");
125 840 : for (i=0;i<contour->npoints;i++) {
126 800 : p_id = i*contour->subcomm->n + contour->subcomm->color;
127 800 : PetscCall(MatDuplicate(T,MAT_DO_NOT_COPY_VALUES,&Amat));
128 800 : if (T != P) PetscCall(MatDuplicate(P,MAT_DO_NOT_COPY_VALUES,&Pmat)); else Pmat = Amat;
129 800 : if (contour->subcomm->n == 1 || nep->fui==NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPComputeFunction(nep,ctx->omega[p_id],Amat,Pmat));
130 112 : else PetscCall(NEPComputeFunctionSubcomm(nep,ctx->omega[p_id],Amat,Pmat,PETSC_FALSE));
131 800 : PetscCall(NEP_KSPSetOperators(contour->ksp[i],Amat,Pmat));
132 800 : PetscCall(MatDestroy(&Amat));
133 800 : if (T != P) PetscCall(MatDestroy(&Pmat));
134 : }
135 40 : PetscFunctionReturn(PETSC_SUCCESS);
136 : }
137 :
138 : /*
139 : Y_i = F(z_i)^{-1}Fp(z_i)V for every integration point, Y=[Y_i] is in the context
140 : */
141 41 : static PetscErrorCode NEPCISSSolve(NEP nep,Mat dT,BV V,PetscInt L_start,PetscInt L_end)
142 : {
143 41 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
144 41 : SlepcContourData contour = ctx->contour;
145 41 : PetscInt i,p_id;
146 41 : Mat MV,BMV=NULL,MC;
147 :
148 41 : PetscFunctionBegin;
149 41 : PetscCall(BVSetActiveColumns(V,L_start,L_end));
150 41 : PetscCall(BVGetMat(V,&MV));
151 873 : for (i=0;i<contour->npoints;i++) {
152 832 : p_id = i*contour->subcomm->n + contour->subcomm->color;
153 832 : if (contour->subcomm->n==1 || nep->fui==NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPComputeJacobian(nep,ctx->omega[p_id],dT));
154 112 : else PetscCall(NEPComputeFunctionSubcomm(nep,ctx->omega[p_id],dT,NULL,PETSC_TRUE));
155 832 : PetscCall(BVSetActiveColumns(ctx->Y,i*ctx->L+L_start,i*ctx->L+L_end));
156 832 : PetscCall(BVGetMat(ctx->Y,&MC));
157 832 : if (!i) {
158 41 : PetscCall(MatProductCreate(dT,MV,NULL,&BMV));
159 41 : PetscCall(MatProductSetType(BMV,MATPRODUCT_AB));
160 41 : PetscCall(MatProductSetFromOptions(BMV));
161 41 : PetscCall(MatProductSymbolic(BMV));
162 : }
163 832 : PetscCall(MatProductNumeric(BMV));
164 832 : PetscCall(KSPMatSolve(contour->ksp[i],BMV,MC));
165 832 : PetscCall(BVRestoreMat(ctx->Y,&MC));
166 : }
167 41 : PetscCall(MatDestroy(&BMV));
168 41 : PetscCall(BVRestoreMat(V,&MV));
169 41 : PetscFunctionReturn(PETSC_SUCCESS);
170 : }
171 :
172 40 : static PetscErrorCode NEPSetUp_CISS(NEP nep)
173 : {
174 40 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
175 40 : SlepcContourData contour;
176 40 : PetscInt nwork;
177 40 : PetscBool istrivial,isellipse,flg;
178 40 : NEP_CISS_PROJECT dsctxf;
179 40 : PetscObjectId id;
180 40 : PetscObjectState state;
181 40 : Vec v0;
182 :
183 40 : PetscFunctionBegin;
184 40 : if (nep->ncv==PETSC_DETERMINE) nep->ncv = ctx->L_max*ctx->M;
185 : else {
186 18 : ctx->L_max = nep->ncv/ctx->M;
187 18 : if (!ctx->L_max) {
188 0 : ctx->L_max = 1;
189 0 : nep->ncv = ctx->L_max*ctx->M;
190 : }
191 : }
192 40 : ctx->L = PetscMin(ctx->L,ctx->L_max);
193 40 : if (nep->max_it==PETSC_DETERMINE) nep->max_it = 5;
194 40 : if (nep->mpd==PETSC_DETERMINE) nep->mpd = nep->ncv;
195 40 : if (!nep->which) nep->which = NEP_ALL;
196 40 : PetscCheck(nep->which==NEP_ALL,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only computing all eigenvalues");
197 40 : NEPCheckUnsupported(nep,NEP_FEATURE_STOPPING | NEP_FEATURE_TWOSIDED);
198 :
199 : /* check region */
200 40 : PetscCall(RGIsTrivial(nep->rg,&istrivial));
201 40 : PetscCheck(!istrivial,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
202 40 : PetscCall(RGGetComplement(nep->rg,&flg));
203 40 : PetscCheck(!flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
204 40 : PetscCall(PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse));
205 40 : PetscCheck(isellipse,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Currently only implemented for elliptic regions");
206 :
207 : /* if the region has changed, then reset contour data */
208 40 : PetscCall(PetscObjectGetId((PetscObject)nep->rg,&id));
209 40 : PetscCall(PetscObjectStateGet((PetscObject)nep->rg,&state));
210 40 : if (ctx->rgid && (id != ctx->rgid || state != ctx->rgstate)) {
211 0 : PetscCall(SlepcContourDataDestroy(&ctx->contour));
212 0 : PetscCall(PetscInfo(nep,"Resetting the contour data structure due to a change of region\n"));
213 0 : ctx->rgid = id; ctx->rgstate = state;
214 : }
215 :
216 : /* create contour data structure */
217 40 : if (!ctx->contour) {
218 0 : PetscCall(RGCanUseConjugates(nep->rg,ctx->isreal,&ctx->useconj));
219 0 : PetscCall(SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)nep,&ctx->contour));
220 : }
221 :
222 40 : PetscCall(NEPAllocateSolution(nep,0));
223 40 : if (ctx->weight) PetscCall(PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma));
224 40 : PetscCall(PetscMalloc4(ctx->N,&ctx->weight,ctx->N,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma));
225 :
226 : /* allocate basis vectors */
227 40 : PetscCall(BVDestroy(&ctx->S));
228 40 : PetscCall(BVDuplicateResize(nep->V,ctx->L*ctx->M,&ctx->S));
229 40 : PetscCall(BVDestroy(&ctx->V));
230 40 : PetscCall(BVDuplicateResize(nep->V,ctx->L,&ctx->V));
231 :
232 40 : contour = ctx->contour;
233 40 : if (contour->subcomm && contour->subcomm->n != 1 && nep->fui==NEP_USER_INTERFACE_CALLBACK) {
234 6 : PetscCall(NEPComputeFunction(nep,0,nep->function,nep->function_pre));
235 6 : PetscCall(SlepcContourRedundantMat(contour,1,&nep->function,(nep->function!=nep->function_pre)?&nep->function_pre:NULL));
236 34 : } else PetscCall(SlepcContourRedundantMat(contour,nep->nt,nep->A,nep->P));
237 40 : if (contour->pA) {
238 18 : if (!ctx->J) PetscCall(MatDuplicate(contour->pA[0],MAT_DO_NOT_COPY_VALUES,&ctx->J));
239 18 : PetscCall(BVGetColumn(ctx->V,0,&v0));
240 18 : PetscCall(SlepcContourScatterCreate(contour,v0));
241 18 : PetscCall(BVRestoreColumn(ctx->V,0,&v0));
242 18 : PetscCall(BVDestroy(&ctx->pV));
243 18 : PetscCall(BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->pV));
244 18 : PetscCall(BVSetSizesFromVec(ctx->pV,contour->xsub,nep->n));
245 18 : PetscCall(BVSetFromOptions(ctx->pV));
246 18 : PetscCall(BVResize(ctx->pV,ctx->L,PETSC_FALSE));
247 : }
248 :
249 40 : PetscCall(BVDestroy(&ctx->Y));
250 40 : if (contour->pA) {
251 18 : PetscCall(BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->Y));
252 18 : PetscCall(BVSetSizesFromVec(ctx->Y,contour->xsub,nep->n));
253 18 : PetscCall(BVSetFromOptions(ctx->Y));
254 18 : PetscCall(BVResize(ctx->Y,contour->npoints*ctx->L,PETSC_FALSE));
255 22 : } else PetscCall(BVDuplicateResize(nep->V,contour->npoints*ctx->L,&ctx->Y));
256 :
257 40 : if (ctx->extraction == NEP_CISS_EXTRACTION_RITZ) {
258 24 : PetscCall(DSSetMethod(nep->ds,1));
259 24 : PetscCall(DSNEPSetRG(nep->ds,nep->rg));
260 24 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscCall(DSNEPSetFN(nep->ds,nep->nt,nep->f));
261 : else {
262 6 : PetscCall(PetscNew(&dsctxf));
263 6 : PetscCall(DSNEPSetComputeMatrixFunction(nep->ds,NEPContourDSComputeMatrix,dsctxf));
264 6 : dsctxf->nep = nep;
265 6 : ctx->dsctxf = dsctxf;
266 : }
267 : }
268 40 : PetscCall(DSAllocate(nep->ds,nep->ncv));
269 40 : nwork = (nep->fui==NEP_USER_INTERFACE_SPLIT)? 2: 1;
270 40 : PetscCall(NEPSetWorkVecs(nep,nwork));
271 40 : PetscFunctionReturn(PETSC_SUCCESS);
272 : }
273 :
274 40 : static PetscErrorCode NEPSolve_CISS(NEP nep)
275 : {
276 40 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
277 40 : SlepcContourData contour = ctx->contour;
278 40 : Mat X,M,E,T,P,J;
279 40 : BV V;
280 40 : PetscInt i,j,ld,L_add=0,nv=0,L_base=ctx->L,inner,*inside;
281 40 : PetscScalar *Mu,*H0,*H1,*rr,*temp,center;
282 40 : PetscReal error,max_error,radius,rgscale,est_eig,eta;
283 40 : PetscBool isellipse,*fl1;
284 40 : Vec si;
285 40 : SlepcSC sc;
286 40 : PetscRandom rand;
287 :
288 40 : PetscFunctionBegin;
289 40 : PetscCall(DSGetSlepcSC(nep->ds,&sc));
290 40 : sc->comparison = SlepcCompareLargestMagnitude;
291 40 : sc->comparisonctx = NULL;
292 40 : sc->map = NULL;
293 40 : sc->mapobj = NULL;
294 40 : PetscCall(DSGetLeadingDimension(nep->ds,&ld));
295 40 : PetscCall(RGComputeQuadrature(nep->rg,RG_QUADRULE_TRAPEZOIDAL,ctx->N,ctx->omega,ctx->pp,ctx->weight));
296 40 : if (contour->pA) {
297 18 : T = contour->pA[0];
298 18 : P = ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && contour->pP))? contour->pP[0]: T;
299 : } else {
300 22 : T = nep->function;
301 22 : P = nep->function_pre? nep->function_pre: nep->function;
302 : }
303 40 : PetscCall(NEPCISSSetUp(nep,T,P));
304 40 : PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
305 40 : PetscCall(BVSetRandomSign(ctx->V));
306 40 : PetscCall(BVGetRandomContext(ctx->V,&rand));
307 40 : if (contour->pA) {
308 18 : J = ctx->J;
309 18 : V = ctx->pV;
310 18 : PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
311 : } else {
312 22 : J = nep->jacobian;
313 22 : V = ctx->V;
314 : }
315 40 : PetscCall(NEPCISSSolve(nep,J,V,0,ctx->L));
316 40 : PetscCall(PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse));
317 40 : if (isellipse) {
318 40 : PetscCall(BVTraceQuadrature(ctx->Y,ctx->V,ctx->L,ctx->L,ctx->weight,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj,&est_eig));
319 40 : PetscCall(PetscInfo(nep,"Estimated eigenvalue count: %f\n",(double)est_eig));
320 40 : eta = PetscPowReal(10.0,-PetscLog10Real(nep->tol)/ctx->N);
321 40 : L_add = PetscMax(0,(PetscInt)PetscCeilReal((est_eig*eta)/ctx->M)-ctx->L);
322 40 : if (L_add>ctx->L_max-ctx->L) {
323 0 : PetscCall(PetscInfo(nep,"Number of eigenvalues inside the contour path may be too large\n"));
324 0 : L_add = ctx->L_max-ctx->L;
325 : }
326 : }
327 : /* Updates L after estimate the number of eigenvalue */
328 40 : if (L_add>0) {
329 1 : PetscCall(PetscInfo(nep,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by Estimate #Eig\n",ctx->L,ctx->L+L_add));
330 1 : PetscCall(BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints));
331 1 : PetscCall(BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add));
332 1 : PetscCall(BVSetRandomSign(ctx->V));
333 1 : if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
334 1 : ctx->L += L_add;
335 1 : PetscCall(NEPCISSSolve(nep,J,V,ctx->L-L_add,ctx->L));
336 : }
337 :
338 40 : PetscCall(PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0));
339 40 : for (i=0;i<ctx->refine_blocksize;i++) {
340 3 : PetscCall(BVDotQuadrature(ctx->Y,V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj));
341 3 : PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
342 3 : PetscCall(PetscLogEventBegin(NEP_CISS_SVD,nep,0,0,0));
343 3 : PetscCall(SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv));
344 3 : PetscCall(PetscLogEventEnd(NEP_CISS_SVD,nep,0,0,0));
345 3 : if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
346 0 : L_add = L_base;
347 0 : if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
348 0 : PetscCall(PetscInfo(nep,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by SVD(H0)\n",ctx->L,ctx->L+L_add));
349 0 : PetscCall(BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints));
350 0 : PetscCall(BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add));
351 0 : PetscCall(BVSetRandomSign(ctx->V));
352 0 : if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
353 0 : ctx->L += L_add;
354 0 : PetscCall(NEPCISSSolve(nep,J,V,ctx->L-L_add,ctx->L));
355 0 : if (L_add) {
356 0 : PetscCall(PetscFree2(Mu,H0));
357 0 : PetscCall(PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0));
358 : }
359 : }
360 :
361 40 : PetscCall(RGGetScale(nep->rg,&rgscale));
362 40 : PetscCall(RGEllipseGetParameters(nep->rg,¢er,&radius,NULL));
363 :
364 40 : if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) PetscCall(PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1));
365 :
366 80 : while (nep->reason == NEP_CONVERGED_ITERATING) {
367 40 : nep->its++;
368 40 : for (inner=0;inner<=ctx->refine_inner;inner++) {
369 40 : if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
370 8 : PetscCall(BVDotQuadrature(ctx->Y,V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj));
371 8 : PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
372 8 : PetscCall(PetscLogEventBegin(NEP_CISS_SVD,nep,0,0,0));
373 8 : PetscCall(SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv));
374 8 : PetscCall(PetscLogEventEnd(NEP_CISS_SVD,nep,0,0,0));
375 : } else {
376 32 : PetscCall(BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj));
377 : /* compute SVD of S */
378 56 : PetscCall(BVSVDAndRank(ctx->S,ctx->M,ctx->L,ctx->delta,(ctx->extraction==NEP_CISS_EXTRACTION_CAA)?BV_SVD_METHOD_QR_CAA:BV_SVD_METHOD_QR,H0,ctx->sigma,&nv));
379 : }
380 40 : PetscCall(PetscInfo(nep,"Estimated rank: nv = %" PetscInt_FMT "\n",nv));
381 40 : if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
382 0 : PetscCall(BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj));
383 0 : PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
384 0 : PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
385 0 : PetscCall(BVCopy(ctx->S,ctx->V));
386 0 : if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
387 0 : PetscCall(NEPCISSSolve(nep,J,V,0,ctx->L));
388 : } else break;
389 : }
390 40 : nep->nconv = 0;
391 40 : if (nv == 0) { nep->reason = NEP_CONVERGED_TOL; break; }
392 : else {
393 : /* Extracting eigenpairs */
394 40 : PetscCall(DSSetDimensions(nep->ds,nv,0,0));
395 40 : PetscCall(DSSetState(nep->ds,DS_STATE_RAW));
396 40 : if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
397 8 : PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
398 8 : PetscCall(CISS_BlockHankel(Mu,1,ctx->L,ctx->M,H1));
399 8 : PetscCall(DSGetArray(nep->ds,DS_MAT_A,&temp));
400 51 : for (j=0;j<nv;j++)
401 290 : for (i=0;i<nv;i++)
402 247 : temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
403 8 : PetscCall(DSRestoreArray(nep->ds,DS_MAT_A,&temp));
404 8 : PetscCall(DSGetArray(nep->ds,DS_MAT_B,&temp));
405 51 : for (j=0;j<nv;j++)
406 290 : for (i=0;i<nv;i++)
407 247 : temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
408 8 : PetscCall(DSRestoreArray(nep->ds,DS_MAT_B,&temp));
409 32 : } else if (ctx->extraction == NEP_CISS_EXTRACTION_CAA) {
410 8 : PetscCall(BVSetActiveColumns(ctx->S,0,nv));
411 8 : PetscCall(DSGetArray(nep->ds,DS_MAT_A,&temp));
412 50 : for (i=0;i<nv;i++) PetscCall(PetscArraycpy(temp+i*ld,H0+i*nv,nv));
413 8 : PetscCall(DSRestoreArray(nep->ds,DS_MAT_A,&temp));
414 : } else {
415 24 : PetscCall(BVSetActiveColumns(ctx->S,0,nv));
416 24 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
417 66 : for (i=0;i<nep->nt;i++) {
418 48 : PetscCall(DSGetMat(nep->ds,DSMatExtra[i],&E));
419 48 : PetscCall(BVMatProject(ctx->S,nep->A[i],ctx->S,E));
420 48 : PetscCall(DSRestoreMat(nep->ds,DSMatExtra[i],&E));
421 : }
422 6 : } else { ctx->dsctxf->Q = ctx->S; }
423 : }
424 40 : PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
425 40 : PetscCall(DSSynchronize(nep->ds,nep->eigr,nep->eigi));
426 40 : PetscCall(DSGetDimensions(nep->ds,NULL,NULL,NULL,&nv));
427 40 : if (ctx->extraction == NEP_CISS_EXTRACTION_CAA || ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
428 101 : for (i=0;i<nv;i++) {
429 85 : nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
430 : }
431 : }
432 40 : PetscCall(PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr));
433 40 : PetscCall(DSVectors(nep->ds,DS_MAT_X,NULL,NULL));
434 40 : PetscCall(DSGetMat(nep->ds,DS_MAT_X,&X));
435 40 : PetscCall(SlepcCISS_isGhost(X,nv,ctx->sigma,ctx->spurious_threshold,fl1));
436 40 : PetscCall(DSRestoreMat(nep->ds,DS_MAT_X,&X));
437 40 : PetscCall(RGCheckInside(nep->rg,nv,nep->eigr,nep->eigi,inside));
438 228 : for (i=0;i<nv;i++) {
439 188 : if (fl1[i] && inside[i]>=0) {
440 167 : rr[i] = 1.0;
441 167 : nep->nconv++;
442 21 : } else rr[i] = 0.0;
443 : }
444 40 : PetscCall(DSSort(nep->ds,nep->eigr,nep->eigi,rr,NULL,&nep->nconv));
445 40 : PetscCall(DSSynchronize(nep->ds,nep->eigr,nep->eigi));
446 40 : if (ctx->extraction == NEP_CISS_EXTRACTION_CAA || ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
447 101 : for (i=0;i<nv;i++) nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
448 : }
449 40 : PetscCall(PetscFree3(fl1,inside,rr));
450 40 : PetscCall(BVSetActiveColumns(nep->V,0,nv));
451 40 : PetscCall(DSVectors(nep->ds,DS_MAT_X,NULL,NULL));
452 40 : if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
453 8 : PetscCall(BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj));
454 8 : PetscCall(BVSetActiveColumns(ctx->S,0,nv));
455 8 : PetscCall(BVCopy(ctx->S,nep->V));
456 8 : PetscCall(DSGetMat(nep->ds,DS_MAT_X,&X));
457 8 : PetscCall(BVMultInPlace(ctx->S,X,0,nep->nconv));
458 8 : PetscCall(BVMultInPlace(nep->V,X,0,nep->nconv));
459 8 : PetscCall(DSRestoreMat(nep->ds,DS_MAT_X,&X));
460 : } else {
461 32 : PetscCall(DSGetMat(nep->ds,DS_MAT_X,&X));
462 32 : PetscCall(BVMultInPlace(ctx->S,X,0,nep->nconv));
463 32 : PetscCall(DSRestoreMat(nep->ds,DS_MAT_X,&X));
464 32 : PetscCall(BVSetActiveColumns(ctx->S,0,nep->nconv));
465 32 : PetscCall(BVCopy(ctx->S,nep->V));
466 : }
467 : max_error = 0.0;
468 207 : for (i=0;i<nep->nconv;i++) {
469 167 : PetscCall(BVGetColumn(nep->V,i,&si));
470 167 : PetscCall(VecNormalize(si,NULL));
471 167 : PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,nep->eigr[i],si,nep->work,&error));
472 167 : PetscCall((*nep->converged)(nep,nep->eigr[i],0,error,&error,nep->convergedctx));
473 167 : PetscCall(BVRestoreColumn(nep->V,i,&si));
474 225 : max_error = PetscMax(max_error,error);
475 : }
476 40 : if (max_error <= nep->tol) nep->reason = NEP_CONVERGED_TOL;
477 0 : else if (nep->its > nep->max_it) nep->reason = NEP_DIVERGED_ITS;
478 : else {
479 0 : if (nep->nconv > ctx->L) nv = nep->nconv;
480 0 : else if (ctx->L > nv) nv = ctx->L;
481 0 : nv = PetscMin(nv,ctx->L*ctx->M);
482 0 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M));
483 0 : PetscCall(MatSetRandom(M,rand));
484 0 : PetscCall(BVSetActiveColumns(ctx->S,0,nv));
485 0 : PetscCall(BVMultInPlace(ctx->S,M,0,ctx->L));
486 0 : PetscCall(MatDestroy(&M));
487 0 : PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
488 0 : PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
489 0 : PetscCall(BVCopy(ctx->S,ctx->V));
490 0 : if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
491 80 : PetscCall(NEPCISSSolve(nep,J,V,0,ctx->L));
492 : }
493 : }
494 : }
495 40 : PetscCall(PetscFree2(Mu,H0));
496 40 : if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) PetscCall(PetscFree(H1));
497 40 : PetscFunctionReturn(PETSC_SUCCESS);
498 : }
499 :
500 38 : static PetscErrorCode NEPCISSSetSizes_CISS(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
501 : {
502 38 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
503 38 : PetscInt oN,oL,oM,oLmax,onpart;
504 38 : PetscMPIInt size;
505 :
506 38 : PetscFunctionBegin;
507 38 : oN = ctx->N;
508 38 : if (ip == PETSC_DETERMINE) {
509 13 : if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
510 25 : } else if (ip != PETSC_CURRENT) {
511 25 : PetscCheck(ip>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
512 25 : PetscCheck(ip%2==0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
513 25 : if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
514 : }
515 38 : oL = ctx->L;
516 38 : if (bs == PETSC_DETERMINE) {
517 13 : ctx->L = 16;
518 25 : } else if (bs != PETSC_CURRENT) {
519 25 : PetscCheck(bs>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
520 25 : ctx->L = bs;
521 : }
522 38 : oM = ctx->M;
523 38 : if (ms == PETSC_DETERMINE) {
524 13 : ctx->M = ctx->N/4;
525 25 : } else if (ms != PETSC_CURRENT) {
526 25 : PetscCheck(ms>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
527 25 : PetscCheck(ms<=ctx->N,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be less than or equal to the number of integration points");
528 25 : ctx->M = ms;
529 : }
530 38 : onpart = ctx->npart;
531 38 : if (npart == PETSC_DETERMINE) {
532 0 : ctx->npart = 1;
533 38 : } else if (npart != PETSC_CURRENT) {
534 38 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)nep),&size));
535 38 : PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
536 38 : ctx->npart = npart;
537 : }
538 38 : oLmax = ctx->L_max;
539 38 : if (bsmax == PETSC_DETERMINE) {
540 13 : ctx->L_max = 64;
541 25 : } else if (bsmax != PETSC_CURRENT) {
542 25 : PetscCheck(bsmax>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
543 25 : ctx->L_max = PetscMax(bsmax,ctx->L);
544 : }
545 38 : if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) {
546 31 : PetscCall(SlepcContourDataDestroy(&ctx->contour));
547 31 : PetscCall(PetscInfo(nep,"Resetting the contour data structure due to a change of parameters\n"));
548 31 : nep->state = NEP_STATE_INITIAL;
549 : }
550 38 : ctx->isreal = realmats;
551 38 : if (oL != ctx->L || oM != ctx->M || oLmax != ctx->L_max) nep->state = NEP_STATE_INITIAL;
552 38 : PetscFunctionReturn(PETSC_SUCCESS);
553 : }
554 :
555 : /*@
556 : NEPCISSSetSizes - Sets the values of various size parameters in the CISS solver.
557 :
558 : Logically Collective
559 :
560 : Input Parameters:
561 : + nep - the nonlinear eigensolver context
562 : . ip - number of integration points
563 : . bs - block size
564 : . ms - moment size
565 : . npart - number of partitions when splitting the communicator
566 : . bsmax - max block size
567 : - realmats - T(z) is real for real z
568 :
569 : Options Database Keys:
570 : + -nep_ciss_integration_points - Sets the number of integration points
571 : . -nep_ciss_blocksize - Sets the block size
572 : . -nep_ciss_moments - Sets the moment size
573 : . -nep_ciss_partitions - Sets the number of partitions
574 : . -nep_ciss_maxblocksize - Sets the maximum block size
575 : - -nep_ciss_realmats - T(z) is real for real z
576 :
577 : Notes:
578 : For all integer arguments, you can use PETSC_CURRENT to keep the current value, and
579 : PETSC_DETERMINE to set them to a default value.
580 :
581 : The default number of partitions is 1. This means the internal KSP object is shared
582 : among all processes of the NEP communicator. Otherwise, the communicator is split
583 : into npart communicators, so that npart KSP solves proceed simultaneously.
584 :
585 : The realmats flag can be set to true when T(.) is guaranteed to be real
586 : when the argument is a real value, for example, when all matrices in
587 : the split form are real. When set to true, the solver avoids some computations.
588 :
589 : Level: advanced
590 :
591 : .seealso: NEPCISSGetSizes()
592 : @*/
593 38 : PetscErrorCode NEPCISSSetSizes(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
594 : {
595 38 : PetscFunctionBegin;
596 38 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
597 114 : PetscValidLogicalCollectiveInt(nep,ip,2);
598 114 : PetscValidLogicalCollectiveInt(nep,bs,3);
599 114 : PetscValidLogicalCollectiveInt(nep,ms,4);
600 114 : PetscValidLogicalCollectiveInt(nep,npart,5);
601 114 : PetscValidLogicalCollectiveInt(nep,bsmax,6);
602 114 : PetscValidLogicalCollectiveBool(nep,realmats,7);
603 38 : PetscTryMethod(nep,"NEPCISSSetSizes_C",(NEP,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(nep,ip,bs,ms,npart,bsmax,realmats));
604 38 : PetscFunctionReturn(PETSC_SUCCESS);
605 : }
606 :
607 38 : static PetscErrorCode NEPCISSGetSizes_CISS(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
608 : {
609 38 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
610 :
611 38 : PetscFunctionBegin;
612 38 : if (ip) *ip = ctx->N;
613 38 : if (bs) *bs = ctx->L;
614 38 : if (ms) *ms = ctx->M;
615 38 : if (npart) *npart = ctx->npart;
616 38 : if (bsmax) *bsmax = ctx->L_max;
617 38 : if (realmats) *realmats = ctx->isreal;
618 38 : PetscFunctionReturn(PETSC_SUCCESS);
619 : }
620 :
621 : /*@
622 : NEPCISSGetSizes - Gets the values of various size parameters in the CISS solver.
623 :
624 : Not Collective
625 :
626 : Input Parameter:
627 : . nep - the nonlinear eigensolver context
628 :
629 : Output Parameters:
630 : + ip - number of integration points
631 : . bs - block size
632 : . ms - moment size
633 : . npart - number of partitions when splitting the communicator
634 : . bsmax - max block size
635 : - realmats - T(z) is real for real z
636 :
637 : Level: advanced
638 :
639 : .seealso: NEPCISSSetSizes()
640 : @*/
641 38 : PetscErrorCode NEPCISSGetSizes(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
642 : {
643 38 : PetscFunctionBegin;
644 38 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
645 38 : PetscUseMethod(nep,"NEPCISSGetSizes_C",(NEP,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(nep,ip,bs,ms,npart,bsmax,realmats));
646 38 : PetscFunctionReturn(PETSC_SUCCESS);
647 : }
648 :
649 3 : static PetscErrorCode NEPCISSSetThreshold_CISS(NEP nep,PetscReal delta,PetscReal spur)
650 : {
651 3 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
652 :
653 3 : PetscFunctionBegin;
654 3 : if (delta == (PetscReal)PETSC_DETERMINE) {
655 0 : ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
656 3 : } else if (delta != (PetscReal)PETSC_CURRENT) {
657 3 : PetscCheck(delta>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
658 3 : ctx->delta = delta;
659 : }
660 3 : if (spur == (PetscReal)PETSC_DETERMINE) {
661 0 : ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
662 3 : } else if (spur != (PetscReal)PETSC_CURRENT) {
663 3 : PetscCheck(spur>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
664 3 : ctx->spurious_threshold = spur;
665 : }
666 3 : PetscFunctionReturn(PETSC_SUCCESS);
667 : }
668 :
669 : /*@
670 : NEPCISSSetThreshold - Sets the values of various threshold parameters in
671 : the CISS solver.
672 :
673 : Logically Collective
674 :
675 : Input Parameters:
676 : + nep - the nonlinear eigensolver context
677 : . delta - threshold for numerical rank
678 : - spur - spurious threshold (to discard spurious eigenpairs)
679 :
680 : Options Database Keys:
681 : + -nep_ciss_delta - Sets the delta
682 : - -nep_ciss_spurious_threshold - Sets the spurious threshold
683 :
684 : Note:
685 : PETSC_CURRENT can be used to preserve the current value of any of the
686 : arguments, and PETSC_DETERMINE to set them to a default value.
687 :
688 : Level: advanced
689 :
690 : .seealso: NEPCISSGetThreshold()
691 : @*/
692 3 : PetscErrorCode NEPCISSSetThreshold(NEP nep,PetscReal delta,PetscReal spur)
693 : {
694 3 : PetscFunctionBegin;
695 3 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
696 9 : PetscValidLogicalCollectiveReal(nep,delta,2);
697 9 : PetscValidLogicalCollectiveReal(nep,spur,3);
698 3 : PetscTryMethod(nep,"NEPCISSSetThreshold_C",(NEP,PetscReal,PetscReal),(nep,delta,spur));
699 3 : PetscFunctionReturn(PETSC_SUCCESS);
700 : }
701 :
702 38 : static PetscErrorCode NEPCISSGetThreshold_CISS(NEP nep,PetscReal *delta,PetscReal *spur)
703 : {
704 38 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
705 :
706 38 : PetscFunctionBegin;
707 38 : if (delta) *delta = ctx->delta;
708 38 : if (spur) *spur = ctx->spurious_threshold;
709 38 : PetscFunctionReturn(PETSC_SUCCESS);
710 : }
711 :
712 : /*@
713 : NEPCISSGetThreshold - Gets the values of various threshold parameters in
714 : the CISS solver.
715 :
716 : Not Collective
717 :
718 : Input Parameter:
719 : . nep - the nonlinear eigensolver context
720 :
721 : Output Parameters:
722 : + delta - threshold for numerical rank
723 : - spur - spurious threshold (to discard spurious eigenpairs)
724 :
725 : Level: advanced
726 :
727 : .seealso: NEPCISSSetThreshold()
728 : @*/
729 38 : PetscErrorCode NEPCISSGetThreshold(NEP nep,PetscReal *delta,PetscReal *spur)
730 : {
731 38 : PetscFunctionBegin;
732 38 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
733 38 : PetscUseMethod(nep,"NEPCISSGetThreshold_C",(NEP,PetscReal*,PetscReal*),(nep,delta,spur));
734 38 : PetscFunctionReturn(PETSC_SUCCESS);
735 : }
736 :
737 3 : static PetscErrorCode NEPCISSSetRefinement_CISS(NEP nep,PetscInt inner,PetscInt blsize)
738 : {
739 3 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
740 :
741 3 : PetscFunctionBegin;
742 3 : if (inner == PETSC_DETERMINE) {
743 0 : ctx->refine_inner = 0;
744 3 : } else if (inner != PETSC_CURRENT) {
745 3 : PetscCheck(inner>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
746 3 : ctx->refine_inner = inner;
747 : }
748 3 : if (blsize == PETSC_DETERMINE) {
749 0 : ctx->refine_blocksize = 0;
750 3 : } else if (blsize != PETSC_CURRENT) {
751 3 : PetscCheck(blsize>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
752 3 : ctx->refine_blocksize = blsize;
753 : }
754 3 : PetscFunctionReturn(PETSC_SUCCESS);
755 : }
756 :
757 : /*@
758 : NEPCISSSetRefinement - Sets the values of various refinement parameters
759 : in the CISS solver.
760 :
761 : Logically Collective
762 :
763 : Input Parameters:
764 : + nep - the nonlinear eigensolver context
765 : . inner - number of iterative refinement iterations (inner loop)
766 : - blsize - number of iterative refinement iterations (blocksize loop)
767 :
768 : Options Database Keys:
769 : + -nep_ciss_refine_inner - Sets number of inner iterations
770 : - -nep_ciss_refine_blocksize - Sets number of blocksize iterations
771 :
772 : Note:
773 : PETSC_CURRENT can be used to preserve the current value of any of the
774 : arguments, and PETSC_DETERMINE to set them to a default value.
775 :
776 : Level: advanced
777 :
778 : .seealso: NEPCISSGetRefinement()
779 : @*/
780 3 : PetscErrorCode NEPCISSSetRefinement(NEP nep,PetscInt inner,PetscInt blsize)
781 : {
782 3 : PetscFunctionBegin;
783 3 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
784 9 : PetscValidLogicalCollectiveInt(nep,inner,2);
785 9 : PetscValidLogicalCollectiveInt(nep,blsize,3);
786 3 : PetscTryMethod(nep,"NEPCISSSetRefinement_C",(NEP,PetscInt,PetscInt),(nep,inner,blsize));
787 3 : PetscFunctionReturn(PETSC_SUCCESS);
788 : }
789 :
790 38 : static PetscErrorCode NEPCISSGetRefinement_CISS(NEP nep,PetscInt *inner,PetscInt *blsize)
791 : {
792 38 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
793 :
794 38 : PetscFunctionBegin;
795 38 : if (inner) *inner = ctx->refine_inner;
796 38 : if (blsize) *blsize = ctx->refine_blocksize;
797 38 : PetscFunctionReturn(PETSC_SUCCESS);
798 : }
799 :
800 : /*@
801 : NEPCISSGetRefinement - Gets the values of various refinement parameters
802 : in the CISS solver.
803 :
804 : Not Collective
805 :
806 : Input Parameter:
807 : . nep - the nonlinear eigensolver context
808 :
809 : Output Parameters:
810 : + inner - number of iterative refinement iterations (inner loop)
811 : - blsize - number of iterative refinement iterations (blocksize loop)
812 :
813 : Level: advanced
814 :
815 : .seealso: NEPCISSSetRefinement()
816 : @*/
817 38 : PetscErrorCode NEPCISSGetRefinement(NEP nep, PetscInt *inner, PetscInt *blsize)
818 : {
819 38 : PetscFunctionBegin;
820 38 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
821 38 : PetscUseMethod(nep,"NEPCISSGetRefinement_C",(NEP,PetscInt*,PetscInt*),(nep,inner,blsize));
822 38 : PetscFunctionReturn(PETSC_SUCCESS);
823 : }
824 :
825 24 : static PetscErrorCode NEPCISSSetExtraction_CISS(NEP nep,NEPCISSExtraction extraction)
826 : {
827 24 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
828 :
829 24 : PetscFunctionBegin;
830 24 : if (ctx->extraction != extraction) {
831 16 : ctx->extraction = extraction;
832 16 : nep->state = NEP_STATE_INITIAL;
833 : }
834 24 : PetscFunctionReturn(PETSC_SUCCESS);
835 : }
836 :
837 : /*@
838 : NEPCISSSetExtraction - Sets the extraction technique used in the CISS solver.
839 :
840 : Logically Collective
841 :
842 : Input Parameters:
843 : + nep - the nonlinear eigensolver context
844 : - extraction - the extraction technique
845 :
846 : Options Database Key:
847 : . -nep_ciss_extraction - Sets the extraction technique (either 'ritz', 'hankel' or 'caa')
848 :
849 : Notes:
850 : By default, the Rayleigh-Ritz extraction is used (NEP_CISS_EXTRACTION_RITZ).
851 :
852 : If the 'hankel' or the 'caa' option is specified (NEP_CISS_EXTRACTION_HANKEL or
853 : NEP_CISS_EXTRACTION_CAA), then the Block Hankel method, or the Communication-avoiding
854 : Arnoldi method, respectively, is used for extracting eigenpairs.
855 :
856 : Level: advanced
857 :
858 : .seealso: NEPCISSGetExtraction(), NEPCISSExtraction
859 : @*/
860 24 : PetscErrorCode NEPCISSSetExtraction(NEP nep,NEPCISSExtraction extraction)
861 : {
862 24 : PetscFunctionBegin;
863 24 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
864 72 : PetscValidLogicalCollectiveEnum(nep,extraction,2);
865 24 : PetscTryMethod(nep,"NEPCISSSetExtraction_C",(NEP,NEPCISSExtraction),(nep,extraction));
866 24 : PetscFunctionReturn(PETSC_SUCCESS);
867 : }
868 :
869 13 : static PetscErrorCode NEPCISSGetExtraction_CISS(NEP nep,NEPCISSExtraction *extraction)
870 : {
871 13 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
872 :
873 13 : PetscFunctionBegin;
874 13 : *extraction = ctx->extraction;
875 13 : PetscFunctionReturn(PETSC_SUCCESS);
876 : }
877 :
878 : /*@
879 : NEPCISSGetExtraction - Gets the extraction technique used in the CISS solver.
880 :
881 : Not Collective
882 :
883 : Input Parameter:
884 : . nep - the nonlinear eigensolver context
885 :
886 : Output Parameters:
887 : . extraction - extraction technique
888 :
889 : Level: advanced
890 :
891 : .seealso: NEPCISSSetExtraction() NEPCISSExtraction
892 : @*/
893 13 : PetscErrorCode NEPCISSGetExtraction(NEP nep,NEPCISSExtraction *extraction)
894 : {
895 13 : PetscFunctionBegin;
896 13 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
897 13 : PetscAssertPointer(extraction,2);
898 13 : PetscUseMethod(nep,"NEPCISSGetExtraction_C",(NEP,NEPCISSExtraction*),(nep,extraction));
899 13 : PetscFunctionReturn(PETSC_SUCCESS);
900 : }
901 :
902 44 : static PetscErrorCode NEPCISSGetKSPs_CISS(NEP nep,PetscInt *nsolve,KSP **ksp)
903 : {
904 44 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
905 44 : SlepcContourData contour;
906 44 : PetscInt i;
907 44 : PC pc;
908 44 : MPI_Comm child;
909 :
910 44 : PetscFunctionBegin;
911 44 : if (!ctx->contour) { /* initialize contour data structure first */
912 44 : PetscCall(RGCanUseConjugates(nep->rg,ctx->isreal,&ctx->useconj));
913 44 : PetscCall(SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)nep,&ctx->contour));
914 : }
915 44 : contour = ctx->contour;
916 44 : if (!contour->ksp) {
917 44 : PetscCall(PetscMalloc1(contour->npoints,&contour->ksp));
918 44 : PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
919 876 : for (i=0;i<contour->npoints;i++) {
920 832 : PetscCall(KSPCreate(child,&contour->ksp[i]));
921 832 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)nep,1));
922 832 : PetscCall(KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)nep)->prefix));
923 832 : PetscCall(KSPAppendOptionsPrefix(contour->ksp[i],"nep_ciss_"));
924 832 : PetscCall(PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)nep)->options));
925 832 : PetscCall(KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE));
926 1248 : PetscCall(KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(nep->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
927 832 : PetscCall(KSPGetPC(contour->ksp[i],&pc));
928 832 : if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
929 96 : PetscCall(KSPSetType(contour->ksp[i],KSPBCGS));
930 96 : PetscCall(PCSetType(pc,PCBJACOBI));
931 : } else {
932 736 : PetscCall(KSPSetType(contour->ksp[i],KSPPREONLY));
933 832 : PetscCall(PCSetType(pc,PCLU));
934 : }
935 : }
936 : }
937 44 : if (nsolve) *nsolve = contour->npoints;
938 44 : if (ksp) *ksp = contour->ksp;
939 44 : PetscFunctionReturn(PETSC_SUCCESS);
940 : }
941 :
942 : /*@C
943 : NEPCISSGetKSPs - Retrieve the array of linear solver objects associated with
944 : the CISS solver.
945 :
946 : Collective
947 :
948 : Input Parameter:
949 : . nep - nonlinear eigenvalue solver
950 :
951 : Output Parameters:
952 : + nsolve - number of solver objects
953 : - ksp - array of linear solver object
954 :
955 : Notes:
956 : The number of KSP solvers is equal to the number of integration points divided by
957 : the number of partitions. This value is halved in the case of real matrices with
958 : a region centered at the real axis.
959 :
960 : Level: advanced
961 :
962 : .seealso: NEPCISSSetSizes()
963 : @*/
964 44 : PetscErrorCode NEPCISSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
965 : {
966 44 : PetscFunctionBegin;
967 44 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
968 44 : PetscUseMethod(nep,"NEPCISSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
969 44 : PetscFunctionReturn(PETSC_SUCCESS);
970 : }
971 :
972 40 : static PetscErrorCode NEPReset_CISS(NEP nep)
973 : {
974 40 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
975 :
976 40 : PetscFunctionBegin;
977 40 : PetscCall(BVDestroy(&ctx->S));
978 40 : PetscCall(BVDestroy(&ctx->V));
979 40 : PetscCall(BVDestroy(&ctx->Y));
980 40 : PetscCall(SlepcContourDataReset(ctx->contour));
981 40 : PetscCall(MatDestroy(&ctx->J));
982 40 : PetscCall(BVDestroy(&ctx->pV));
983 40 : if (ctx->extraction == NEP_CISS_EXTRACTION_RITZ && nep->fui==NEP_USER_INTERFACE_CALLBACK) PetscCall(PetscFree(ctx->dsctxf));
984 40 : PetscFunctionReturn(PETSC_SUCCESS);
985 : }
986 :
987 38 : static PetscErrorCode NEPSetFromOptions_CISS(NEP nep,PetscOptionItems *PetscOptionsObject)
988 : {
989 38 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
990 38 : PetscReal r1,r2;
991 38 : PetscInt i,i1,i2,i3,i4,i5,i6,i7;
992 38 : PetscBool b1,flg,flg2,flg3,flg4,flg5,flg6;
993 38 : NEPCISSExtraction extraction;
994 :
995 38 : PetscFunctionBegin;
996 38 : PetscOptionsHeadBegin(PetscOptionsObject,"NEP CISS Options");
997 :
998 38 : PetscCall(NEPCISSGetSizes(nep,&i1,&i2,&i3,&i4,&i5,&b1));
999 38 : PetscCall(PetscOptionsInt("-nep_ciss_integration_points","Number of integration points","NEPCISSSetSizes",i1,&i1,&flg));
1000 38 : PetscCall(PetscOptionsInt("-nep_ciss_blocksize","Block size","NEPCISSSetSizes",i2,&i2,&flg2));
1001 38 : PetscCall(PetscOptionsInt("-nep_ciss_moments","Moment size","NEPCISSSetSizes",i3,&i3,&flg3));
1002 38 : PetscCall(PetscOptionsInt("-nep_ciss_partitions","Number of partitions","NEPCISSSetSizes",i4,&i4,&flg4));
1003 38 : PetscCall(PetscOptionsInt("-nep_ciss_maxblocksize","Maximum block size","NEPCISSSetSizes",i5,&i5,&flg5));
1004 38 : PetscCall(PetscOptionsBool("-nep_ciss_realmats","True if T(z) is real for real z","NEPCISSSetSizes",b1,&b1,&flg6));
1005 38 : if (flg || flg2 || flg3 || flg4 || flg5 || flg6) PetscCall(NEPCISSSetSizes(nep,i1,i2,i3,i4,i5,b1));
1006 :
1007 38 : PetscCall(NEPCISSGetThreshold(nep,&r1,&r2));
1008 38 : PetscCall(PetscOptionsReal("-nep_ciss_delta","Threshold for numerical rank","NEPCISSSetThreshold",r1,&r1,&flg));
1009 38 : PetscCall(PetscOptionsReal("-nep_ciss_spurious_threshold","Threshold for the spurious eigenpairs","NEPCISSSetThreshold",r2,&r2,&flg2));
1010 38 : if (flg || flg2) PetscCall(NEPCISSSetThreshold(nep,r1,r2));
1011 :
1012 38 : PetscCall(NEPCISSGetRefinement(nep,&i6,&i7));
1013 38 : PetscCall(PetscOptionsInt("-nep_ciss_refine_inner","Number of inner iterative refinement iterations","NEPCISSSetRefinement",i6,&i6,&flg));
1014 38 : PetscCall(PetscOptionsInt("-nep_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","NEPCISSSetRefinement",i7,&i7,&flg2));
1015 38 : if (flg || flg2) PetscCall(NEPCISSSetRefinement(nep,i6,i7));
1016 :
1017 38 : PetscCall(PetscOptionsEnum("-nep_ciss_extraction","Extraction technique","NEPCISSSetExtraction",NEPCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg));
1018 38 : if (flg) PetscCall(NEPCISSSetExtraction(nep,extraction));
1019 :
1020 38 : PetscOptionsHeadEnd();
1021 :
1022 38 : if (!nep->rg) PetscCall(NEPGetRG(nep,&nep->rg));
1023 38 : PetscCall(RGSetFromOptions(nep->rg)); /* this is necessary here to set useconj */
1024 38 : if (!ctx->contour || !ctx->contour->ksp) PetscCall(NEPCISSGetKSPs(nep,NULL,NULL));
1025 38 : PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Something went wrong with NEPCISSGetKSPs()");
1026 774 : for (i=0;i<ctx->contour->npoints;i++) PetscCall(KSPSetFromOptions(ctx->contour->ksp[i]));
1027 38 : PetscCall(PetscSubcommSetFromOptions(ctx->contour->subcomm));
1028 38 : PetscFunctionReturn(PETSC_SUCCESS);
1029 : }
1030 :
1031 38 : static PetscErrorCode NEPDestroy_CISS(NEP nep)
1032 : {
1033 38 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
1034 :
1035 38 : PetscFunctionBegin;
1036 38 : PetscCall(SlepcContourDataDestroy(&ctx->contour));
1037 38 : PetscCall(PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma));
1038 38 : PetscCall(PetscFree(nep->data));
1039 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NULL));
1040 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NULL));
1041 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NULL));
1042 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NULL));
1043 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NULL));
1044 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NULL));
1045 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetExtraction_C",NULL));
1046 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetExtraction_C",NULL));
1047 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NULL));
1048 38 : PetscFunctionReturn(PETSC_SUCCESS);
1049 : }
1050 :
1051 1 : static PetscErrorCode NEPView_CISS(NEP nep,PetscViewer viewer)
1052 : {
1053 1 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
1054 1 : PetscBool isascii;
1055 1 : PetscViewer sviewer;
1056 :
1057 1 : PetscFunctionBegin;
1058 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1059 1 : if (isascii) {
1060 1 : PetscCall(PetscViewerASCIIPrintf(viewer," sizes { integration points: %" PetscInt_FMT ", block size: %" PetscInt_FMT ", moment size: %" PetscInt_FMT ", partitions: %" PetscInt_FMT ", maximum block size: %" PetscInt_FMT " }\n",ctx->N,ctx->L,ctx->M,ctx->npart,ctx->L_max));
1061 1 : if (ctx->isreal) PetscCall(PetscViewerASCIIPrintf(viewer," exploiting symmetry of integration points\n"));
1062 1 : PetscCall(PetscViewerASCIIPrintf(viewer," threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold));
1063 1 : PetscCall(PetscViewerASCIIPrintf(viewer," iterative refinement { inner: %" PetscInt_FMT ", blocksize: %" PetscInt_FMT " }\n",ctx->refine_inner, ctx->refine_blocksize));
1064 1 : PetscCall(PetscViewerASCIIPrintf(viewer," extraction: %s\n",NEPCISSExtractions[ctx->extraction]));
1065 1 : if (!ctx->contour || !ctx->contour->ksp) PetscCall(NEPCISSGetKSPs(nep,NULL,NULL));
1066 1 : PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Something went wrong with NEPCISSGetKSPs()");
1067 1 : PetscCall(PetscViewerASCIIPushTab(viewer));
1068 1 : if (ctx->npart>1 && ctx->contour->subcomm) {
1069 0 : PetscCall(PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
1070 0 : if (!ctx->contour->subcomm->color) PetscCall(KSPView(ctx->contour->ksp[0],sviewer));
1071 0 : PetscCall(PetscViewerFlush(sviewer));
1072 0 : PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
1073 : /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1074 0 : PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1075 1 : } else PetscCall(KSPView(ctx->contour->ksp[0],viewer));
1076 1 : PetscCall(PetscViewerASCIIPopTab(viewer));
1077 : }
1078 1 : PetscFunctionReturn(PETSC_SUCCESS);
1079 : }
1080 :
1081 78 : static PetscErrorCode NEPSetDSType_CISS(NEP nep)
1082 : {
1083 78 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
1084 :
1085 78 : PetscFunctionBegin;
1086 78 : switch (ctx->extraction) {
1087 46 : case NEP_CISS_EXTRACTION_RITZ: PetscCall(DSSetType(nep->ds,DSNEP)); break;
1088 16 : case NEP_CISS_EXTRACTION_HANKEL: PetscCall(DSSetType(nep->ds,DSGNHEP)); break;
1089 16 : case NEP_CISS_EXTRACTION_CAA: PetscCall(DSSetType(nep->ds,DSNHEP)); break;
1090 : }
1091 78 : PetscFunctionReturn(PETSC_SUCCESS);
1092 : }
1093 :
1094 38 : SLEPC_EXTERN PetscErrorCode NEPCreate_CISS(NEP nep)
1095 : {
1096 38 : NEP_CISS *ctx = (NEP_CISS*)nep->data;
1097 :
1098 38 : PetscFunctionBegin;
1099 38 : PetscCall(PetscNew(&ctx));
1100 38 : nep->data = ctx;
1101 : /* set default values of parameters */
1102 38 : ctx->N = 32;
1103 38 : ctx->L = 16;
1104 38 : ctx->M = ctx->N/4;
1105 38 : ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
1106 38 : ctx->L_max = 64;
1107 38 : ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
1108 38 : ctx->isreal = PETSC_FALSE;
1109 38 : ctx->npart = 1;
1110 :
1111 38 : nep->useds = PETSC_TRUE;
1112 :
1113 38 : nep->ops->solve = NEPSolve_CISS;
1114 38 : nep->ops->setup = NEPSetUp_CISS;
1115 38 : nep->ops->setfromoptions = NEPSetFromOptions_CISS;
1116 38 : nep->ops->reset = NEPReset_CISS;
1117 38 : nep->ops->destroy = NEPDestroy_CISS;
1118 38 : nep->ops->view = NEPView_CISS;
1119 38 : nep->ops->setdstype = NEPSetDSType_CISS;
1120 :
1121 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NEPCISSSetSizes_CISS));
1122 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NEPCISSGetSizes_CISS));
1123 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NEPCISSSetThreshold_CISS));
1124 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NEPCISSGetThreshold_CISS));
1125 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NEPCISSSetRefinement_CISS));
1126 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NEPCISSGetRefinement_CISS));
1127 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetExtraction_C",NEPCISSSetExtraction_CISS));
1128 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetExtraction_C",NEPCISSGetExtraction_CISS));
1129 38 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NEPCISSGetKSPs_CISS));
1130 38 : PetscFunctionReturn(PETSC_SUCCESS);
1131 : }
|
