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 : #include <slepc/private/dsimpl.h> /*I "slepcds.h" I*/
12 : #include <slepcblaslapack.h>
13 :
14 : typedef struct {
15 : PetscInt nf; /* number of functions in f[] */
16 : FN f[DS_NUM_EXTRA]; /* functions defining the nonlinear operator */
17 : PetscInt max_mid; /* maximum minimality index */
18 : PetscInt nnod; /* number of nodes for quadrature rules */
19 : PetscInt spls; /* number of sampling columns for quadrature rules */
20 : PetscInt Nit; /* number of refinement iterations */
21 : PetscReal rtol; /* tolerance of Newton refinement */
22 : RG rg; /* region for contour integral */
23 : PetscLayout map; /* used to distribute work among MPI processes */
24 : void *computematrixctx;
25 : DSNEPMatrixFunctionFn *computematrix;
26 : } DS_NEP;
27 :
28 : /*
29 : DSNEPComputeMatrix - Build the matrix associated with a nonlinear operator
30 : T(lambda) or its derivative T'(lambda), given the parameter lambda, where
31 : T(lambda) = sum_i E_i*f_i(lambda). The result is written in mat.
32 : */
33 969 : static PetscErrorCode DSNEPComputeMatrix(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat)
34 : {
35 969 : DS_NEP *ctx = (DS_NEP*)ds->data;
36 969 : PetscScalar *T,alpha;
37 969 : const PetscScalar *E;
38 969 : PetscInt i,ld,n;
39 969 : PetscBLASInt k,inc=1;
40 :
41 969 : PetscFunctionBegin;
42 969 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
43 969 : if (ctx->computematrix) PetscCall((*ctx->computematrix)(ds,lambda,deriv,mat,ctx->computematrixctx));
44 : else {
45 10 : PetscCall(DSGetDimensions(ds,&n,NULL,NULL,NULL));
46 10 : PetscCall(DSGetLeadingDimension(ds,&ld));
47 10 : PetscCall(PetscBLASIntCast(ld*n,&k));
48 10 : PetscCall(MatDenseGetArray(ds->omat[mat],&T));
49 10 : PetscCall(PetscArrayzero(T,k));
50 40 : for (i=0;i<ctx->nf;i++) {
51 30 : if (deriv) PetscCall(FNEvaluateDerivative(ctx->f[i],lambda,&alpha));
52 18 : else PetscCall(FNEvaluateFunction(ctx->f[i],lambda,&alpha));
53 30 : PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&E));
54 30 : PetscCallBLAS("BLASaxpy",BLASaxpy_(&k,&alpha,E,&inc,T,&inc));
55 30 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&E));
56 : }
57 10 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&T));
58 : }
59 969 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
60 969 : PetscFunctionReturn(PETSC_SUCCESS);
61 : }
62 :
63 20 : static PetscErrorCode DSAllocate_NEP(DS ds,PetscInt ld)
64 : {
65 20 : DS_NEP *ctx = (DS_NEP*)ds->data;
66 20 : PetscInt i;
67 :
68 20 : PetscFunctionBegin;
69 20 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
70 80 : for (i=0;i<ctx->nf;i++) PetscCall(DSAllocateMat_Private(ds,DSMatExtra[i]));
71 20 : PetscCall(PetscFree(ds->perm));
72 20 : PetscCall(PetscMalloc1(ld*ctx->max_mid,&ds->perm));
73 20 : PetscFunctionReturn(PETSC_SUCCESS);
74 : }
75 :
76 1 : static PetscErrorCode DSView_NEP(DS ds,PetscViewer viewer)
77 : {
78 1 : DS_NEP *ctx = (DS_NEP*)ds->data;
79 1 : PetscViewerFormat format;
80 1 : PetscInt i;
81 1 : const char *methodname[] = {
82 : "Successive Linear Problems",
83 : "Contour Integral"
84 : };
85 1 : const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
86 :
87 1 : PetscFunctionBegin;
88 1 : PetscCall(PetscViewerGetFormat(viewer,&format));
89 1 : if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
90 1 : if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
91 : #if defined(PETSC_USE_COMPLEX)
92 : if (ds->method==1) { /* contour integral method */
93 : PetscCall(PetscViewerASCIIPrintf(viewer,"number of integration points: %" PetscInt_FMT "\n",ctx->nnod));
94 : PetscCall(PetscViewerASCIIPrintf(viewer,"maximum minimality index: %" PetscInt_FMT "\n",ctx->max_mid));
95 : if (ctx->spls) PetscCall(PetscViewerASCIIPrintf(viewer,"number of sampling columns for quadrature: %" PetscInt_FMT "\n",ctx->spls));
96 : if (ctx->Nit) PetscCall(PetscViewerASCIIPrintf(viewer,"doing iterative refinement (%" PetscInt_FMT " its, tolerance %g)\n",ctx->Nit,(double)ctx->rtol));
97 : PetscCall(RGView(ctx->rg,viewer));
98 : }
99 : #endif
100 1 : if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(PetscViewerASCIIPrintf(viewer,"number of functions: %" PetscInt_FMT "\n",ctx->nf));
101 1 : PetscFunctionReturn(PETSC_SUCCESS);
102 : }
103 0 : for (i=0;i<ctx->nf;i++) {
104 0 : PetscCall(FNView(ctx->f[i],viewer));
105 0 : PetscCall(DSViewMat(ds,viewer,DSMatExtra[i]));
106 : }
107 0 : if (ds->state>DS_STATE_INTERMEDIATE) PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
108 0 : PetscFunctionReturn(PETSC_SUCCESS);
109 : }
110 :
111 2 : static PetscErrorCode DSVectors_NEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
112 : {
113 2 : PetscFunctionBegin;
114 2 : PetscCheck(!rnorm,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
115 2 : switch (mat) {
116 : case DS_MAT_X:
117 2 : break;
118 0 : case DS_MAT_Y:
119 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
120 0 : default:
121 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
122 : }
123 2 : PetscFunctionReturn(PETSC_SUCCESS);
124 : }
125 :
126 2 : static PetscErrorCode DSSort_NEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *dummy)
127 : {
128 2 : DS_NEP *ctx = (DS_NEP*)ds->data;
129 2 : PetscInt n,l,i,*perm,lds;
130 2 : PetscScalar *Q;
131 :
132 2 : PetscFunctionBegin;
133 2 : if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
134 2 : if (!ds->method) PetscFunctionReturn(PETSC_SUCCESS); /* SLP computes just one eigenvalue */
135 0 : n = ds->n*ctx->max_mid;
136 0 : lds = ds->ld*ctx->max_mid;
137 0 : l = ds->l;
138 0 : perm = ds->perm;
139 0 : for (i=0;i<n;i++) perm[i] = i;
140 0 : if (rr) PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
141 0 : else PetscCall(DSSortEigenvalues_Private(ds,wr,NULL,perm,PETSC_FALSE));
142 0 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
143 0 : for (i=l;i<ds->t;i++) Q[i+i*lds] = wr[perm[i]];
144 0 : for (i=l;i<ds->t;i++) wr[i] = Q[i+i*lds];
145 0 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
146 : /* n != ds->n */
147 0 : PetscCall(DSPermuteColumns_Private(ds,0,ds->t,ds->n,DS_MAT_X,perm));
148 0 : PetscFunctionReturn(PETSC_SUCCESS);
149 : }
150 :
151 : #if defined(SLEPC_MISSING_LAPACK_GGEV3)
152 : #define LAPGEEV "ggev"
153 : #else
154 : #define LAPGEEV "ggev3"
155 : #endif
156 :
157 185 : static PetscErrorCode DSSolve_NEP_SLP(DS ds,PetscScalar *wr,PetscScalar *wi)
158 : {
159 185 : PetscScalar *A,*B,*W,*X,*work,*alpha,*beta,a;
160 185 : PetscScalar sigma,lambda,mu,re,re2,sone=1.0,szero=0.0;
161 185 : PetscBLASInt info,n,ld,lwork,one=1,zero=0;
162 185 : PetscInt it,pos,j,maxit=100,result;
163 185 : PetscReal norm,tol,done=1.0;
164 : #if !defined(PETSC_USE_COMPLEX)
165 185 : PetscReal *alphai,im,im2;
166 : #endif
167 :
168 185 : PetscFunctionBegin;
169 185 : PetscCall(PetscBLASIntCast(ds->n,&n));
170 185 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
171 185 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
172 185 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
173 185 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
174 185 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
175 185 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
176 185 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
177 185 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
178 :
179 : /* workspace query and memory allocation */
180 185 : lwork = -1;
181 : #if defined(PETSC_USE_COMPLEX)
182 : PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&n,A,&ld,B,&ld,NULL,NULL,NULL,&ld,W,&ld,&a,&lwork,NULL,&info));
183 : PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
184 : PetscCall(DSAllocateWork_Private(ds,lwork+2*ds->n,8*ds->n,0));
185 : alpha = ds->work;
186 : beta = ds->work + ds->n;
187 : work = ds->work + 2*ds->n;
188 : #else
189 185 : PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&n,A,&ld,B,&ld,NULL,NULL,NULL,NULL,&ld,W,&ld,&a,&lwork,&info));
190 185 : PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
191 185 : PetscCall(DSAllocateWork_Private(ds,lwork+3*ds->n,0,0));
192 185 : alpha = ds->work;
193 185 : beta = ds->work + ds->n;
194 185 : alphai = ds->work + 2*ds->n;
195 185 : work = ds->work + 3*ds->n;
196 : #endif
197 :
198 185 : sigma = 0.0;
199 185 : if (ds->sc->comparison==SlepcCompareTargetMagnitude || ds->sc->comparison==SlepcCompareTargetReal) sigma = *(PetscScalar*)ds->sc->comparisonctx;
200 185 : lambda = sigma;
201 185 : tol = n*PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);
202 :
203 577 : for (it=0;it<maxit;it++) {
204 :
205 : /* evaluate T and T' */
206 577 : PetscCall(DSNEPComputeMatrix(ds,lambda,PETSC_FALSE,DS_MAT_A));
207 577 : if (it) {
208 392 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,A,&ld,X,&one,&szero,X+ld,&one));
209 392 : norm = BLASnrm2_(&n,X+ld,&one);
210 392 : if (norm/PetscAbsScalar(lambda)<=tol) break;
211 : }
212 392 : PetscCall(DSNEPComputeMatrix(ds,lambda,PETSC_TRUE,DS_MAT_B));
213 :
214 : /* compute eigenvalue correction mu and eigenvector u */
215 : #if defined(PETSC_USE_COMPLEX)
216 : PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&n,A,&ld,B,&ld,alpha,beta,NULL,&ld,W,&ld,work,&lwork,ds->rwork,&info));
217 : #else
218 392 : PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&n,A,&ld,B,&ld,alpha,alphai,beta,NULL,&ld,W,&ld,work,&lwork,&info));
219 : #endif
220 392 : SlepcCheckLapackInfo(LAPGEEV,info);
221 :
222 : /* find smallest eigenvalue */
223 392 : j = 0;
224 392 : if (beta[j]==0.0) re = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
225 392 : else re = alpha[j]/beta[j];
226 : #if !defined(PETSC_USE_COMPLEX)
227 392 : if (beta[j]==0.0) im = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
228 392 : else im = alphai[j]/beta[j];
229 : #endif
230 : pos = 0;
231 1861 : for (j=1;j<n;j++) {
232 1469 : if (beta[j]==0.0) re2 = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
233 1469 : else re2 = alpha[j]/beta[j];
234 : #if !defined(PETSC_USE_COMPLEX)
235 1469 : if (beta[j]==0.0) im2 = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
236 1469 : else im2 = alphai[j]/beta[j];
237 1469 : PetscCall(SlepcCompareSmallestMagnitude(re,im,re2,im2,&result,NULL));
238 : #else
239 : PetscCall(SlepcCompareSmallestMagnitude(re,0.0,re2,0.0,&result,NULL));
240 : #endif
241 1469 : if (result > 0) {
242 938 : re = re2;
243 : #if !defined(PETSC_USE_COMPLEX)
244 938 : im = im2;
245 : #endif
246 938 : pos = j;
247 : }
248 : }
249 :
250 : #if !defined(PETSC_USE_COMPLEX)
251 392 : PetscCheck(im==0.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"DSNEP found a complex eigenvalue; try rerunning with complex scalars");
252 : #endif
253 392 : mu = alpha[pos]/beta[pos];
254 392 : PetscCall(PetscArraycpy(X,W+pos*ld,n));
255 392 : norm = BLASnrm2_(&n,X,&one);
256 392 : PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&one,X,&n,&info));
257 392 : SlepcCheckLapackInfo("lascl",info);
258 :
259 : /* correct eigenvalue approximation */
260 392 : lambda = lambda - mu;
261 : }
262 185 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
263 185 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
264 185 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
265 185 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
266 :
267 185 : PetscCheck(it<maxit,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"DSNEP did not converge");
268 185 : ds->t = 1;
269 185 : wr[0] = lambda;
270 185 : if (wi) wi[0] = 0.0;
271 185 : PetscFunctionReturn(PETSC_SUCCESS);
272 : }
273 :
274 : #if defined(PETSC_USE_COMPLEX)
275 : /*
276 : Newton refinement for eigenpairs computed with contour integral.
277 : k - number of eigenpairs to refine
278 : wr - eigenvalues (eigenvectors are stored in DS_MAT_X)
279 : */
280 : static PetscErrorCode DSNEPNewtonRefine(DS ds,PetscInt k,PetscScalar *wr)
281 : {
282 : DS_NEP *ctx = (DS_NEP*)ds->data;
283 : PetscScalar *X,*W,*U,*R,sone=1.0,szero=0.0;
284 : PetscReal norm;
285 : PetscInt i,j,ii,nwu=0,*p,jstart=0,jend=k;
286 : const PetscInt *range;
287 : PetscBLASInt n,*perm,info,ld,one=1,n1;
288 : PetscMPIInt len,size,root;
289 : PetscLayout map;
290 :
291 : PetscFunctionBegin;
292 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
293 : PetscCall(PetscBLASIntCast(ds->n,&n));
294 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
295 : n1 = n+1;
296 : p = ds->perm;
297 : PetscCall(PetscArrayzero(p,k));
298 : PetscCall(DSAllocateWork_Private(ds,(n+2)*(n+1),0,n+1));
299 : U = ds->work+nwu; nwu += (n+1)*(n+1);
300 : R = ds->work+nwu; /*nwu += n+1;*/
301 : perm = ds->iwork;
302 : if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {
303 : PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)ds),PETSC_DECIDE,k,1,&map));
304 : PetscCall(PetscLayoutGetRange(map,&jstart,&jend));
305 : }
306 : for (ii=0;ii<ctx->Nit;ii++) {
307 : for (j=jstart;j<jend;j++) {
308 : if (p[j]<2) {
309 : PetscCall(DSNEPComputeMatrix(ds,wr[j],PETSC_FALSE,DS_MAT_W));
310 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
311 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,W,&ld,X+ld*j,&one,&szero,R,&one));
312 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
313 : norm = BLASnrm2_(&n,R,&one);
314 : if (norm/PetscAbsScalar(wr[j]) > ctx->rtol) {
315 : PetscCall(PetscInfo(NULL,"Refining eigenpair %" PetscInt_FMT ", residual=%g\n",j,(double)(norm/PetscAbsScalar(wr[j]))));
316 : p[j] = 1;
317 : R[n] = 0.0;
318 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
319 : for (i=0;i<n;i++) {
320 : PetscCall(PetscArraycpy(U+i*n1,W+i*ld,n));
321 : U[n+i*n1] = PetscConj(X[j*ld+i]);
322 : }
323 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
324 : U[n+n*n1] = 0.0;
325 : PetscCall(DSNEPComputeMatrix(ds,wr[j],PETSC_TRUE,DS_MAT_W));
326 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
327 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,W,&ld,X+ld*j,&one,&szero,U+n*(n+1),&one));
328 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
329 : /* solve system */
330 : PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n1,&n1,U,&n1,perm,&info));
331 : SlepcCheckLapackInfo("getrf",info);
332 : PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n1,&one,U,&n1,perm,R,&n1,&info));
333 : SlepcCheckLapackInfo("getrs",info);
334 : wr[j] -= R[n];
335 : for (i=0;i<n;i++) X[j*ld+i] -= R[i];
336 : /* normalization */
337 : norm = BLASnrm2_(&n,X+ld*j,&one);
338 : for (i=0;i<n;i++) X[ld*j+i] /= norm;
339 : } else p[j] = 2;
340 : }
341 : }
342 : }
343 : if (ds->pmode==DS_PARALLEL_DISTRIBUTED) { /* communicate results */
344 : PetscCall(PetscMPIIntCast(k,&len));
345 : PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,p,len,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)ds)));
346 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
347 : PetscCall(PetscLayoutGetRanges(map,&range));
348 : for (j=0;j<k;j++) {
349 : if (p[j]) { /* j-th eigenpair has been refined */
350 : for (root=0;root<size;root++) if (range[root+1]>j) break;
351 : PetscCall(PetscMPIIntCast(1,&len));
352 : PetscCallMPI(MPI_Bcast(wr+j,len,MPIU_SCALAR,root,PetscObjectComm((PetscObject)ds)));
353 : PetscCall(PetscMPIIntCast(n,&len));
354 : PetscCallMPI(MPI_Bcast(X+ld*j,len,MPIU_SCALAR,root,PetscObjectComm((PetscObject)ds)));
355 : }
356 : }
357 : PetscCall(PetscLayoutDestroy(&map));
358 : }
359 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
360 : PetscFunctionReturn(PETSC_SUCCESS);
361 : }
362 :
363 : PetscErrorCode DSSolve_NEP_Contour(DS ds,PetscScalar *wr,PetscScalar *wi)
364 : {
365 : DS_NEP *ctx = (DS_NEP*)ds->data;
366 : PetscScalar *alpha,*beta,*Q,*Z,*X,*U,*V,*W,*work,*Rc,*R,*w,*z,*zn,*S;
367 : PetscScalar sone=1.0,szero=0.0,center,a;
368 : PetscReal *rwork,norm,radius,vscale,rgscale,*sigma;
369 : PetscBLASInt info,n,*perm,p,pp,ld,lwork,k_,rk_,colA,rowA,one=1;
370 : PetscInt mid,lds,nnod=ctx->nnod,k,i,ii,jj,j,s,off,rk,nwu=0,nw,lrwork,*inside,kstart=0,kend=nnod;
371 : PetscMPIInt len;
372 : PetscBool isellipse;
373 : PetscRandom rand;
374 :
375 : PetscFunctionBegin;
376 : PetscCheck(ctx->rg,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"The contour solver requires a region passed with DSNEPSetRG()");
377 : /* Contour parameters */
378 : PetscCall(PetscObjectTypeCompare((PetscObject)ctx->rg,RGELLIPSE,&isellipse));
379 : PetscCheck(isellipse,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Region must be Ellipse");
380 : PetscCall(RGEllipseGetParameters(ctx->rg,¢er,&radius,&vscale));
381 : PetscCall(RGGetScale(ctx->rg,&rgscale));
382 : if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {
383 : if (!ctx->map) PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)ds),PETSC_DECIDE,ctx->nnod,1,&ctx->map));
384 : PetscCall(PetscLayoutGetRange(ctx->map,&kstart,&kend));
385 : }
386 :
387 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_W)); /* size n */
388 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q)); /* size mid*n */
389 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_Z)); /* size mid*n */
390 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_U)); /* size mid*n */
391 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_V)); /* size mid*n */
392 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
393 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Z],&Z));
394 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
395 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
396 : mid = ctx->max_mid;
397 : PetscCall(PetscBLASIntCast(ds->n,&n));
398 : p = n; /* maximum number of columns for the probing matrix */
399 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
400 : PetscCall(PetscBLASIntCast(mid*n,&rowA));
401 : nw = 2*n*(p+mid)+3*nnod+2*mid*n*p;
402 : lrwork = 9*mid*n;
403 :
404 : /* workspace query and memory allocation */
405 : lwork = -1;
406 : PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&rowA,Q,&rowA,Z,&rowA,NULL,NULL,NULL,&ld,V,&rowA,&a,&lwork,NULL,&info));
407 : PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
408 : PetscCall(DSAllocateWork_Private(ds,lwork+nw,lrwork,n+1));
409 :
410 : sigma = ds->rwork;
411 : rwork = ds->rwork+mid*n;
412 : perm = ds->iwork;
413 : z = ds->work+nwu; nwu += nnod; /* quadrature points */
414 : zn = ds->work+nwu; nwu += nnod; /* normalized quadrature points */
415 : w = ds->work+nwu; nwu += nnod; /* quadrature weights */
416 : Rc = ds->work+nwu; nwu += n*p;
417 : R = ds->work+nwu; nwu += n*p;
418 : alpha = ds->work+nwu; nwu += mid*n;
419 : beta = ds->work+nwu; nwu += mid*n;
420 : S = ds->work+nwu; nwu += 2*mid*n*p;
421 : work = ds->work+nwu;
422 :
423 : /* Compute quadrature parameters */
424 : PetscCall(RGComputeQuadrature(ctx->rg,RG_QUADRULE_TRAPEZOIDAL,nnod,z,zn,w));
425 :
426 : /* Set random matrix */
427 : PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)ds),&rand));
428 : PetscCall(PetscRandomSetSeed(rand,0x12345678));
429 : PetscCall(PetscRandomSeed(rand));
430 : for (j=0;j<p;j++)
431 : for (i=0;i<n;i++) PetscCall(PetscRandomGetValue(rand,Rc+i+j*n));
432 : PetscCall(PetscArrayzero(S,2*mid*n*p));
433 : /* Loop of integration points */
434 : for (k=kstart;k<kend;k++) {
435 : PetscCall(PetscInfo(NULL,"Solving integration point %" PetscInt_FMT "\n",k));
436 : PetscCall(PetscArraycpy(R,Rc,p*n));
437 : PetscCall(DSNEPComputeMatrix(ds,z[k],PETSC_FALSE,DS_MAT_W));
438 :
439 : /* LU factorization */
440 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
441 : PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n,&n,W,&ld,perm,&info));
442 : SlepcCheckLapackInfo("getrf",info);
443 : PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&p,W,&ld,perm,R,&n,&info));
444 : SlepcCheckLapackInfo("getrs",info);
445 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
446 :
447 : /* Moments computation */
448 : for (s=0;s<2*ctx->max_mid;s++) {
449 : off = s*n*p;
450 : for (j=0;j<p;j++)
451 : for (i=0;i<n;i++) S[off+i+j*n] += w[k]*R[j*n+i];
452 : w[k] *= zn[k];
453 : }
454 : }
455 :
456 : if (ds->pmode==DS_PARALLEL_DISTRIBUTED) { /* compute final S via reduction */
457 : PetscCall(PetscMPIIntCast(2*mid*n*p,&len));
458 : PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,S,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ds)));
459 : }
460 : PetscCall(PetscBLASIntCast(ctx->spls?PetscMin(ctx->spls,n):n,&p));
461 : pp = p;
462 : do {
463 : p = pp;
464 : PetscCall(PetscBLASIntCast(mid*p,&colA));
465 :
466 : PetscCall(PetscInfo(ds,"Computing SVD of size %" PetscBLASInt_FMT "x%" PetscBLASInt_FMT "\n",rowA,colA));
467 : for (jj=0;jj<mid;jj++) {
468 : for (ii=0;ii<mid;ii++) {
469 : off = jj*p*rowA+ii*n;
470 : for (j=0;j<p;j++)
471 : for (i=0;i<n;i++) Q[off+j*rowA+i] = S[((jj+ii)*n+j)*n+i];
472 : }
473 : }
474 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&rowA,&colA,Q,&rowA,sigma,U,&rowA,V,&colA,work,&lwork,rwork,&info));
475 : SlepcCheckLapackInfo("gesvd",info);
476 :
477 : rk = colA;
478 : for (i=1;i<colA;i++) if (sigma[i]/sigma[0]<PETSC_MACHINE_EPSILON*1e4) {rk = i; break;}
479 : if (rk<colA || p==n) break;
480 : pp *= 2;
481 : } while (pp<=n);
482 : PetscCall(PetscInfo(ds,"Solving generalized eigenproblem of size %" PetscInt_FMT "\n",rk));
483 : for (jj=0;jj<mid;jj++) {
484 : for (ii=0;ii<mid;ii++) {
485 : off = jj*p*rowA+ii*n;
486 : for (j=0;j<p;j++)
487 : for (i=0;i<n;i++) Q[off+j*rowA+i] = S[((jj+ii+1)*n+j)*n+i];
488 : }
489 : }
490 : PetscCall(PetscBLASIntCast(rk,&rk_));
491 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&rowA,&rk_,&colA,&sone,Q,&rowA,V,&colA,&szero,Z,&rowA));
492 : PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&rk_,&rk_,&rowA,&sone,U,&rowA,Z,&rowA,&szero,Q,&rk_));
493 : PetscCall(PetscArrayzero(Z,n*mid*n*mid));
494 : for (j=0;j<rk;j++) Z[j+j*rk_] = sigma[j];
495 : PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&rk_,Q,&rk_,Z,&rk_,alpha,beta,NULL,&ld,V,&rk_,work,&lwork,rwork,&info));
496 : SlepcCheckLapackInfo(LAPGEEV,info);
497 : for (i=0;i<rk;i++) wr[i] = (center+radius*alpha[i]/beta[i])*rgscale;
498 : PetscCall(PetscMalloc1(rk,&inside));
499 : PetscCall(RGCheckInside(ctx->rg,rk,wr,wi,inside));
500 : k=0;
501 : for (i=0;i<rk;i++)
502 : if (inside[i]==1) inside[k++] = i;
503 : /* Discard values outside region */
504 : lds = ld*mid;
505 : PetscCall(PetscArrayzero(Q,lds*lds));
506 : PetscCall(PetscArrayzero(Z,lds*lds));
507 : for (i=0;i<k;i++) Q[i+i*lds] = (center*beta[inside[i]]+radius*alpha[inside[i]])*rgscale;
508 : for (i=0;i<k;i++) Z[i+i*lds] = beta[inside[i]];
509 : for (i=0;i<k;i++) wr[i] = Q[i+i*lds]/Z[i+i*lds];
510 : for (j=0;j<k;j++) for (i=0;i<rk;i++) V[j*rk+i] = sigma[i]*V[inside[j]*rk+i];
511 : PetscCall(PetscBLASIntCast(k,&k_));
512 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
513 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k_,&rk_,&sone,U,&rowA,V,&rk_,&szero,X,&ld));
514 : /* Normalize */
515 : for (j=0;j<k;j++) {
516 : norm = BLASnrm2_(&n,X+ld*j,&one);
517 : for (i=0;i<n;i++) X[ld*j+i] /= norm;
518 : }
519 : PetscCall(PetscFree(inside));
520 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
521 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
522 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Z],&Z));
523 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
524 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
525 :
526 : /* Newton refinement */
527 : if (ctx->Nit) PetscCall(DSNEPNewtonRefine(ds,k,wr));
528 : ds->t = k;
529 : PetscCall(PetscRandomDestroy(&rand));
530 : PetscFunctionReturn(PETSC_SUCCESS);
531 : }
532 : #endif
533 :
534 : #if !defined(PETSC_HAVE_MPIUNI)
535 44 : static PetscErrorCode DSSynchronize_NEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
536 : {
537 44 : DS_NEP *ctx = (DS_NEP*)ds->data;
538 44 : PetscInt ld=ds->ld,k=0;
539 44 : PetscMPIInt n,n2,rank,size,off=0;
540 44 : PetscScalar *X;
541 :
542 44 : PetscFunctionBegin;
543 44 : if (!ds->method) { /* SLP */
544 44 : if (ds->state>=DS_STATE_CONDENSED) k += ds->n;
545 44 : if (eigr) k += 1;
546 44 : if (eigi) k += 1;
547 44 : PetscCall(PetscMPIIntCast(1,&n));
548 44 : PetscCall(PetscMPIIntCast(ds->n,&n2));
549 : } else { /* Contour */
550 0 : if (ds->state>=DS_STATE_CONDENSED) k += ctx->max_mid*ds->n*ld;
551 0 : if (eigr) k += ctx->max_mid*ds->n;
552 0 : if (eigi) k += ctx->max_mid*ds->n;
553 0 : PetscCall(PetscMPIIntCast(ctx->max_mid*ds->n,&n));
554 0 : PetscCall(PetscMPIIntCast(ctx->max_mid*ds->n*ld,&n2));
555 : }
556 44 : PetscCall(DSAllocateWork_Private(ds,k,0,0));
557 44 : PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
558 44 : if (ds->state>=DS_STATE_CONDENSED) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
559 44 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
560 44 : if (!rank) {
561 22 : if (ds->state>=DS_STATE_CONDENSED) PetscCallMPI(MPI_Pack(X,n2,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
562 22 : if (eigr) PetscCallMPI(MPI_Pack(eigr,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
563 : #if !defined(PETSC_USE_COMPLEX)
564 22 : if (eigi) PetscCallMPI(MPI_Pack(eigi,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
565 : #endif
566 : }
567 88 : PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
568 44 : if (rank) {
569 22 : if (ds->state>=DS_STATE_CONDENSED) PetscCallMPI(MPI_Unpack(ds->work,size,&off,X,n2,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
570 22 : if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
571 : #if !defined(PETSC_USE_COMPLEX)
572 22 : if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
573 : #endif
574 : }
575 44 : if (ds->state>=DS_STATE_CONDENSED) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
576 44 : PetscFunctionReturn(PETSC_SUCCESS);
577 : }
578 : #endif
579 :
580 20 : static PetscErrorCode DSNEPSetFN_NEP(DS ds,PetscInt n,FN fn[])
581 : {
582 20 : DS_NEP *ctx = (DS_NEP*)ds->data;
583 20 : PetscInt i;
584 :
585 20 : PetscFunctionBegin;
586 20 : PetscCheck(n>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more functions, you have %" PetscInt_FMT,n);
587 20 : PetscCheck(n<=DS_NUM_EXTRA,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Too many functions, you specified %" PetscInt_FMT " but the limit is %d",n,DS_NUM_EXTRA);
588 20 : if (ds->ld) PetscCall(PetscInfo(ds,"DSNEPSetFN() called after DSAllocate()\n"));
589 80 : for (i=0;i<n;i++) PetscCall(PetscObjectReference((PetscObject)fn[i]));
590 20 : for (i=0;i<ctx->nf;i++) PetscCall(FNDestroy(&ctx->f[i]));
591 80 : for (i=0;i<n;i++) ctx->f[i] = fn[i];
592 20 : ctx->nf = n;
593 20 : PetscFunctionReturn(PETSC_SUCCESS);
594 : }
595 :
596 : /*@
597 : DSNEPSetFN - Sets a number of functions that define the nonlinear
598 : eigenproblem.
599 :
600 : Collective
601 :
602 : Input Parameters:
603 : + ds - the direct solver context
604 : . n - number of functions
605 : - fn - array of functions
606 :
607 : Notes:
608 : The nonlinear eigenproblem is defined in terms of the split nonlinear
609 : operator T(lambda) = sum_i A_i*f_i(lambda).
610 :
611 : This function must be called before DSAllocate(). Then DSAllocate()
612 : will allocate an extra matrix A_i per each function, that can be
613 : filled in the usual way.
614 :
615 : Level: advanced
616 :
617 : .seealso: DSNEPGetFN(), DSAllocate()
618 : @*/
619 20 : PetscErrorCode DSNEPSetFN(DS ds,PetscInt n,FN fn[])
620 : {
621 20 : PetscInt i;
622 :
623 20 : PetscFunctionBegin;
624 20 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
625 60 : PetscValidLogicalCollectiveInt(ds,n,2);
626 20 : PetscAssertPointer(fn,3);
627 80 : for (i=0;i<n;i++) {
628 60 : PetscValidHeaderSpecific(fn[i],FN_CLASSID,3);
629 60 : PetscCheckSameComm(ds,1,fn[i],3);
630 : }
631 20 : PetscTryMethod(ds,"DSNEPSetFN_C",(DS,PetscInt,FN[]),(ds,n,fn));
632 20 : PetscFunctionReturn(PETSC_SUCCESS);
633 : }
634 :
635 3 : static PetscErrorCode DSNEPGetFN_NEP(DS ds,PetscInt k,FN *fn)
636 : {
637 3 : DS_NEP *ctx = (DS_NEP*)ds->data;
638 :
639 3 : PetscFunctionBegin;
640 3 : PetscCheck(k>=0 && k<ctx->nf,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,ctx->nf-1);
641 3 : *fn = ctx->f[k];
642 3 : PetscFunctionReturn(PETSC_SUCCESS);
643 : }
644 :
645 : /*@
646 : DSNEPGetFN - Gets the functions associated with the nonlinear DS.
647 :
648 : Not Collective
649 :
650 : Input Parameters:
651 : + ds - the direct solver context
652 : - k - the index of the requested function (starting in 0)
653 :
654 : Output Parameter:
655 : . fn - the function
656 :
657 : Level: advanced
658 :
659 : .seealso: DSNEPSetFN()
660 : @*/
661 3 : PetscErrorCode DSNEPGetFN(DS ds,PetscInt k,FN *fn)
662 : {
663 3 : PetscFunctionBegin;
664 3 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
665 3 : PetscAssertPointer(fn,3);
666 3 : PetscUseMethod(ds,"DSNEPGetFN_C",(DS,PetscInt,FN*),(ds,k,fn));
667 3 : PetscFunctionReturn(PETSC_SUCCESS);
668 : }
669 :
670 2 : static PetscErrorCode DSNEPGetNumFN_NEP(DS ds,PetscInt *n)
671 : {
672 2 : DS_NEP *ctx = (DS_NEP*)ds->data;
673 :
674 2 : PetscFunctionBegin;
675 2 : *n = ctx->nf;
676 2 : PetscFunctionReturn(PETSC_SUCCESS);
677 : }
678 :
679 : /*@
680 : DSNEPGetNumFN - Returns the number of functions stored internally by
681 : the DS.
682 :
683 : Not Collective
684 :
685 : Input Parameter:
686 : . ds - the direct solver context
687 :
688 : Output Parameters:
689 : . n - the number of functions passed in DSNEPSetFN()
690 :
691 : Level: advanced
692 :
693 : .seealso: DSNEPSetFN()
694 : @*/
695 2 : PetscErrorCode DSNEPGetNumFN(DS ds,PetscInt *n)
696 : {
697 2 : PetscFunctionBegin;
698 2 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
699 2 : PetscAssertPointer(n,2);
700 2 : PetscUseMethod(ds,"DSNEPGetNumFN_C",(DS,PetscInt*),(ds,n));
701 2 : PetscFunctionReturn(PETSC_SUCCESS);
702 : }
703 :
704 1 : static PetscErrorCode DSNEPSetMinimality_NEP(DS ds,PetscInt n)
705 : {
706 1 : DS_NEP *ctx = (DS_NEP*)ds->data;
707 :
708 1 : PetscFunctionBegin;
709 1 : if (n == PETSC_DECIDE || n == PETSC_DEFAULT) ctx->max_mid = 4;
710 : else {
711 1 : PetscCheck(n>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The minimality value must be > 0");
712 1 : ctx->max_mid = n;
713 : }
714 1 : PetscFunctionReturn(PETSC_SUCCESS);
715 : }
716 :
717 : /*@
718 : DSNEPSetMinimality - Sets the maximum minimality index used internally by
719 : the DSNEP.
720 :
721 : Logically Collective
722 :
723 : Input Parameters:
724 : + ds - the direct solver context
725 : - n - the maximum minimality index
726 :
727 : Options Database Key:
728 : . -ds_nep_minimality <n> - sets the maximum minimality index
729 :
730 : Notes:
731 : The maximum minimality index is used only in the contour integral method,
732 : and is related to the highest momemts used in the method. The default
733 : value is 1, an larger value might give better accuracy in some cases, but
734 : at a higher cost.
735 :
736 : Level: advanced
737 :
738 : .seealso: DSNEPGetMinimality()
739 : @*/
740 1 : PetscErrorCode DSNEPSetMinimality(DS ds,PetscInt n)
741 : {
742 1 : PetscFunctionBegin;
743 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
744 3 : PetscValidLogicalCollectiveInt(ds,n,2);
745 1 : PetscTryMethod(ds,"DSNEPSetMinimality_C",(DS,PetscInt),(ds,n));
746 1 : PetscFunctionReturn(PETSC_SUCCESS);
747 : }
748 :
749 21 : static PetscErrorCode DSNEPGetMinimality_NEP(DS ds,PetscInt *n)
750 : {
751 21 : DS_NEP *ctx = (DS_NEP*)ds->data;
752 :
753 21 : PetscFunctionBegin;
754 21 : *n = ctx->max_mid;
755 21 : PetscFunctionReturn(PETSC_SUCCESS);
756 : }
757 :
758 : /*@
759 : DSNEPGetMinimality - Returns the maximum minimality index used internally by
760 : the DSNEP.
761 :
762 : Not Collective
763 :
764 : Input Parameter:
765 : . ds - the direct solver context
766 :
767 : Output Parameters:
768 : . n - the maximum minimality index passed in DSNEPSetMinimality()
769 :
770 : Level: advanced
771 :
772 : .seealso: DSNEPSetMinimality()
773 : @*/
774 21 : PetscErrorCode DSNEPGetMinimality(DS ds,PetscInt *n)
775 : {
776 21 : PetscFunctionBegin;
777 21 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
778 21 : PetscAssertPointer(n,2);
779 21 : PetscUseMethod(ds,"DSNEPGetMinimality_C",(DS,PetscInt*),(ds,n));
780 21 : PetscFunctionReturn(PETSC_SUCCESS);
781 : }
782 :
783 2 : static PetscErrorCode DSNEPSetRefine_NEP(DS ds,PetscReal tol,PetscInt its)
784 : {
785 2 : DS_NEP *ctx = (DS_NEP*)ds->data;
786 :
787 2 : PetscFunctionBegin;
788 2 : if (tol == (PetscReal)PETSC_DETERMINE) {
789 0 : ctx->rtol = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);
790 2 : } else if (tol != (PetscReal)PETSC_CURRENT) {
791 1 : PetscCheck(tol>0.0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The tolerance must be > 0");
792 1 : ctx->rtol = tol;
793 : }
794 2 : if (its == PETSC_DETERMINE) {
795 1 : ctx->Nit = 3;
796 1 : } else if (its != PETSC_CURRENT) {
797 1 : PetscCheck(its>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The number of iterations must be >= 0");
798 1 : ctx->Nit = its;
799 : }
800 2 : PetscFunctionReturn(PETSC_SUCCESS);
801 : }
802 :
803 : /*@
804 : DSNEPSetRefine - Sets the tolerance and the number of iterations of Newton iterative
805 : refinement for eigenpairs.
806 :
807 : Logically Collective
808 :
809 : Input Parameters:
810 : + ds - the direct solver context
811 : . tol - the tolerance
812 : - its - the number of iterations
813 :
814 : Options Database Key:
815 : + -ds_nep_refine_tol <tol> - sets the tolerance
816 : - -ds_nep_refine_its <its> - sets the number of Newton iterations
817 :
818 : Notes:
819 : Iterative refinement of eigenpairs is currently used only in the contour
820 : integral method.
821 :
822 : Use PETSC_CURRENT to retain the current value of any of the parameters.
823 : Use PETSC_DETERMINE for either argument to assign a default value computed
824 : internally.
825 :
826 : Level: advanced
827 :
828 : .seealso: DSNEPGetRefine()
829 : @*/
830 2 : PetscErrorCode DSNEPSetRefine(DS ds,PetscReal tol,PetscInt its)
831 : {
832 2 : PetscFunctionBegin;
833 2 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
834 6 : PetscValidLogicalCollectiveReal(ds,tol,2);
835 6 : PetscValidLogicalCollectiveInt(ds,its,3);
836 2 : PetscTryMethod(ds,"DSNEPSetRefine_C",(DS,PetscReal,PetscInt),(ds,tol,its));
837 2 : PetscFunctionReturn(PETSC_SUCCESS);
838 : }
839 :
840 1 : static PetscErrorCode DSNEPGetRefine_NEP(DS ds,PetscReal *tol,PetscInt *its)
841 : {
842 1 : DS_NEP *ctx = (DS_NEP*)ds->data;
843 :
844 1 : PetscFunctionBegin;
845 1 : if (tol) *tol = ctx->rtol;
846 1 : if (its) *its = ctx->Nit;
847 1 : PetscFunctionReturn(PETSC_SUCCESS);
848 : }
849 :
850 : /*@
851 : DSNEPGetRefine - Returns the tolerance and the number of iterations of Newton iterative
852 : refinement for eigenpairs.
853 :
854 : Not Collective
855 :
856 : Input Parameter:
857 : . ds - the direct solver context
858 :
859 : Output Parameters:
860 : + tol - the tolerance
861 : - its - the number of iterations
862 :
863 : Level: advanced
864 :
865 : .seealso: DSNEPSetRefine()
866 : @*/
867 1 : PetscErrorCode DSNEPGetRefine(DS ds,PetscReal *tol,PetscInt *its)
868 : {
869 1 : PetscFunctionBegin;
870 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
871 1 : PetscUseMethod(ds,"DSNEPGetRefine_C",(DS,PetscReal*,PetscInt*),(ds,tol,its));
872 1 : PetscFunctionReturn(PETSC_SUCCESS);
873 : }
874 :
875 1 : static PetscErrorCode DSNEPSetIntegrationPoints_NEP(DS ds,PetscInt ip)
876 : {
877 1 : DS_NEP *ctx = (DS_NEP*)ds->data;
878 :
879 1 : PetscFunctionBegin;
880 1 : if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) ctx->nnod = 64;
881 : else {
882 1 : PetscCheck(ip>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The number of integration points must be > 0");
883 1 : ctx->nnod = ip;
884 : }
885 1 : PetscCall(PetscLayoutDestroy(&ctx->map)); /* need to redistribute at next solve */
886 1 : PetscFunctionReturn(PETSC_SUCCESS);
887 : }
888 :
889 : /*@
890 : DSNEPSetIntegrationPoints - Sets the number of integration points to be
891 : used in the contour integral method.
892 :
893 : Logically Collective
894 :
895 : Input Parameters:
896 : + ds - the direct solver context
897 : - ip - the number of integration points
898 :
899 : Options Database Key:
900 : . -ds_nep_integration_points <ip> - sets the number of integration points
901 :
902 : Notes:
903 : This parameter is relevant only in the contour integral method.
904 :
905 : Level: advanced
906 :
907 : .seealso: DSNEPGetIntegrationPoints()
908 : @*/
909 1 : PetscErrorCode DSNEPSetIntegrationPoints(DS ds,PetscInt ip)
910 : {
911 1 : PetscFunctionBegin;
912 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
913 3 : PetscValidLogicalCollectiveInt(ds,ip,2);
914 1 : PetscTryMethod(ds,"DSNEPSetIntegrationPoints_C",(DS,PetscInt),(ds,ip));
915 1 : PetscFunctionReturn(PETSC_SUCCESS);
916 : }
917 :
918 1 : static PetscErrorCode DSNEPGetIntegrationPoints_NEP(DS ds,PetscInt *ip)
919 : {
920 1 : DS_NEP *ctx = (DS_NEP*)ds->data;
921 :
922 1 : PetscFunctionBegin;
923 1 : *ip = ctx->nnod;
924 1 : PetscFunctionReturn(PETSC_SUCCESS);
925 : }
926 :
927 : /*@
928 : DSNEPGetIntegrationPoints - Returns the number of integration points used
929 : in the contour integral method.
930 :
931 : Not Collective
932 :
933 : Input Parameter:
934 : . ds - the direct solver context
935 :
936 : Output Parameters:
937 : . ip - the number of integration points
938 :
939 : Level: advanced
940 :
941 : .seealso: DSNEPSetIntegrationPoints()
942 : @*/
943 1 : PetscErrorCode DSNEPGetIntegrationPoints(DS ds,PetscInt *ip)
944 : {
945 1 : PetscFunctionBegin;
946 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
947 1 : PetscAssertPointer(ip,2);
948 1 : PetscUseMethod(ds,"DSNEPGetIntegrationPoints_C",(DS,PetscInt*),(ds,ip));
949 1 : PetscFunctionReturn(PETSC_SUCCESS);
950 : }
951 :
952 1 : static PetscErrorCode DSNEPSetSamplingSize_NEP(DS ds,PetscInt p)
953 : {
954 1 : DS_NEP *ctx = (DS_NEP*)ds->data;
955 :
956 1 : PetscFunctionBegin;
957 1 : if (p == PETSC_DECIDE || p == PETSC_DEFAULT) ctx->spls = 0;
958 : else {
959 1 : PetscCheck(p>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The sample size must be > 0");
960 1 : PetscCheck(p>=20,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The sample size cannot be smaller than 20");
961 1 : ctx->spls = p;
962 : }
963 1 : PetscFunctionReturn(PETSC_SUCCESS);
964 : }
965 :
966 : /*@
967 : DSNEPSetSamplingSize - Sets the number of sampling columns to be
968 : used in the contour integral method.
969 :
970 : Logically Collective
971 :
972 : Input Parameters:
973 : + ds - the direct solver context
974 : - p - the number of columns for the sampling matrix
975 :
976 : Options Database Key:
977 : . -ds_nep_sampling_size <p> - sets the number of sampling columns
978 :
979 : Notes:
980 : This parameter is relevant only in the contour integral method.
981 :
982 : Level: advanced
983 :
984 : .seealso: DSNEPGetSamplingSize()
985 : @*/
986 1 : PetscErrorCode DSNEPSetSamplingSize(DS ds,PetscInt p)
987 : {
988 1 : PetscFunctionBegin;
989 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
990 3 : PetscValidLogicalCollectiveInt(ds,p,2);
991 1 : PetscTryMethod(ds,"DSNEPSetSamplingSize_C",(DS,PetscInt),(ds,p));
992 1 : PetscFunctionReturn(PETSC_SUCCESS);
993 : }
994 :
995 1 : static PetscErrorCode DSNEPGetSamplingSize_NEP(DS ds,PetscInt *p)
996 : {
997 1 : DS_NEP *ctx = (DS_NEP*)ds->data;
998 :
999 1 : PetscFunctionBegin;
1000 1 : *p = ctx->spls;
1001 1 : PetscFunctionReturn(PETSC_SUCCESS);
1002 : }
1003 :
1004 : /*@
1005 : DSNEPGetSamplingSize - Returns the number of sampling columns used
1006 : in the contour integral method.
1007 :
1008 : Not Collective
1009 :
1010 : Input Parameter:
1011 : . ds - the direct solver context
1012 :
1013 : Output Parameters:
1014 : . p - the number of columns for the sampling matrix
1015 :
1016 : Level: advanced
1017 :
1018 : .seealso: DSNEPSetSamplingSize()
1019 : @*/
1020 1 : PetscErrorCode DSNEPGetSamplingSize(DS ds,PetscInt *p)
1021 : {
1022 1 : PetscFunctionBegin;
1023 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
1024 1 : PetscAssertPointer(p,2);
1025 1 : PetscUseMethod(ds,"DSNEPGetSamplingSize_C",(DS,PetscInt*),(ds,p));
1026 1 : PetscFunctionReturn(PETSC_SUCCESS);
1027 : }
1028 :
1029 17 : static PetscErrorCode DSNEPSetComputeMatrixFunction_NEP(DS ds,DSNEPMatrixFunctionFn *fun,void *ctx)
1030 : {
1031 17 : DS_NEP *dsctx = (DS_NEP*)ds->data;
1032 :
1033 17 : PetscFunctionBegin;
1034 17 : dsctx->computematrix = fun;
1035 17 : dsctx->computematrixctx = ctx;
1036 17 : PetscFunctionReturn(PETSC_SUCCESS);
1037 : }
1038 :
1039 : /*@C
1040 : DSNEPSetComputeMatrixFunction - Sets a user-provided subroutine to compute
1041 : the matrices T(lambda) or T'(lambda).
1042 :
1043 : Logically Collective
1044 :
1045 : Input Parameters:
1046 : + ds - the direct solver context
1047 : . fun - matrix function evaluation routine, see DSNEPMatrixFunctionFn for the calling sequence
1048 : - ctx - a context pointer (the last parameter to the user function)
1049 :
1050 : Note:
1051 : The result is computed as T(lambda) = sum_i E_i*f_i(lambda), and similarly
1052 : for the derivative.
1053 :
1054 : Level: developer
1055 :
1056 : .seealso: DSNEPGetComputeMatrixFunction()
1057 : @*/
1058 17 : PetscErrorCode DSNEPSetComputeMatrixFunction(DS ds,DSNEPMatrixFunctionFn *fun,void *ctx)
1059 : {
1060 17 : PetscFunctionBegin;
1061 17 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
1062 17 : PetscTryMethod(ds,"DSNEPSetComputeMatrixFunction_C",(DS,DSNEPMatrixFunctionFn*,void*),(ds,fun,ctx));
1063 17 : PetscFunctionReturn(PETSC_SUCCESS);
1064 : }
1065 :
1066 0 : static PetscErrorCode DSNEPGetComputeMatrixFunction_NEP(DS ds,DSNEPMatrixFunctionFn **fun,void **ctx)
1067 : {
1068 0 : DS_NEP *dsctx = (DS_NEP*)ds->data;
1069 :
1070 0 : PetscFunctionBegin;
1071 0 : if (fun) *fun = dsctx->computematrix;
1072 0 : if (ctx) *ctx = dsctx->computematrixctx;
1073 0 : PetscFunctionReturn(PETSC_SUCCESS);
1074 : }
1075 :
1076 : /*@C
1077 : DSNEPGetComputeMatrixFunction - Returns the user-provided callback function
1078 : set in DSNEPSetComputeMatrixFunction().
1079 :
1080 : Not Collective
1081 :
1082 : Input Parameter:
1083 : . ds - the direct solver context
1084 :
1085 : Output Parameters:
1086 : + fun - the pointer to the user function
1087 : - ctx - the context pointer
1088 :
1089 : Level: developer
1090 :
1091 : .seealso: DSNEPSetComputeMatrixFunction()
1092 : @*/
1093 0 : PetscErrorCode DSNEPGetComputeMatrixFunction(DS ds,DSNEPMatrixFunctionFn **fun,void **ctx)
1094 : {
1095 0 : PetscFunctionBegin;
1096 0 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
1097 0 : PetscUseMethod(ds,"DSNEPGetComputeMatrixFunction_C",(DS,DSNEPMatrixFunctionFn**,void**),(ds,fun,ctx));
1098 0 : PetscFunctionReturn(PETSC_SUCCESS);
1099 : }
1100 :
1101 1 : static PetscErrorCode DSNEPSetRG_NEP(DS ds,RG rg)
1102 : {
1103 1 : DS_NEP *dsctx = (DS_NEP*)ds->data;
1104 :
1105 1 : PetscFunctionBegin;
1106 1 : PetscCall(PetscObjectReference((PetscObject)rg));
1107 1 : PetscCall(RGDestroy(&dsctx->rg));
1108 1 : dsctx->rg = rg;
1109 1 : PetscFunctionReturn(PETSC_SUCCESS);
1110 : }
1111 :
1112 : /*@
1113 : DSNEPSetRG - Associates a region object to the DSNEP solver.
1114 :
1115 : Collective
1116 :
1117 : Input Parameters:
1118 : + ds - the direct solver context
1119 : - rg - the region context
1120 :
1121 : Notes:
1122 : The region is used only in the contour integral method, and
1123 : should enclose the wanted eigenvalues.
1124 :
1125 : Level: developer
1126 :
1127 : .seealso: DSNEPGetRG()
1128 : @*/
1129 1 : PetscErrorCode DSNEPSetRG(DS ds,RG rg)
1130 : {
1131 1 : PetscFunctionBegin;
1132 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
1133 1 : if (rg) {
1134 1 : PetscValidHeaderSpecific(rg,RG_CLASSID,2);
1135 1 : PetscCheckSameComm(ds,1,rg,2);
1136 : }
1137 1 : PetscTryMethod(ds,"DSNEPSetRG_C",(DS,RG),(ds,rg));
1138 1 : PetscFunctionReturn(PETSC_SUCCESS);
1139 : }
1140 :
1141 1 : static PetscErrorCode DSNEPGetRG_NEP(DS ds,RG *rg)
1142 : {
1143 1 : DS_NEP *ctx = (DS_NEP*)ds->data;
1144 :
1145 1 : PetscFunctionBegin;
1146 1 : if (!ctx->rg) {
1147 1 : PetscCall(RGCreate(PetscObjectComm((PetscObject)ds),&ctx->rg));
1148 1 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->rg,(PetscObject)ds,1));
1149 1 : PetscCall(RGSetOptionsPrefix(ctx->rg,((PetscObject)ds)->prefix));
1150 1 : PetscCall(RGAppendOptionsPrefix(ctx->rg,"ds_nep_"));
1151 1 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->rg,((PetscObject)ds)->options));
1152 : }
1153 1 : *rg = ctx->rg;
1154 1 : PetscFunctionReturn(PETSC_SUCCESS);
1155 : }
1156 :
1157 : /*@
1158 : DSNEPGetRG - Obtain the region object associated to the DSNEP solver.
1159 :
1160 : Collective
1161 :
1162 : Input Parameter:
1163 : . ds - the direct solver context
1164 :
1165 : Output Parameter:
1166 : . rg - the region context
1167 :
1168 : Level: developer
1169 :
1170 : .seealso: DSNEPSetRG()
1171 : @*/
1172 1 : PetscErrorCode DSNEPGetRG(DS ds,RG *rg)
1173 : {
1174 1 : PetscFunctionBegin;
1175 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
1176 1 : PetscAssertPointer(rg,2);
1177 1 : PetscUseMethod(ds,"DSNEPGetRG_C",(DS,RG*),(ds,rg));
1178 1 : PetscFunctionReturn(PETSC_SUCCESS);
1179 : }
1180 :
1181 19 : static PetscErrorCode DSSetFromOptions_NEP(DS ds,PetscOptionItems *PetscOptionsObject)
1182 : {
1183 19 : PetscInt k;
1184 19 : PetscBool flg;
1185 : #if defined(PETSC_USE_COMPLEX)
1186 : PetscReal r;
1187 : PetscBool flg1;
1188 : DS_NEP *ctx = (DS_NEP*)ds->data;
1189 : #endif
1190 :
1191 19 : PetscFunctionBegin;
1192 19 : PetscOptionsHeadBegin(PetscOptionsObject,"DS NEP Options");
1193 :
1194 19 : PetscCall(PetscOptionsInt("-ds_nep_minimality","Maximum minimality index","DSNEPSetMinimality",4,&k,&flg));
1195 19 : if (flg) PetscCall(DSNEPSetMinimality(ds,k));
1196 :
1197 19 : PetscCall(PetscOptionsInt("-ds_nep_integration_points","Number of integration points","DSNEPSetIntegrationPoints",64,&k,&flg));
1198 19 : if (flg) PetscCall(DSNEPSetIntegrationPoints(ds,k));
1199 :
1200 19 : PetscCall(PetscOptionsInt("-ds_nep_sampling_size","Number of sampling columns","DSNEPSetSamplingSize",0,&k,&flg));
1201 19 : if (flg) PetscCall(DSNEPSetSamplingSize(ds,k));
1202 :
1203 : #if defined(PETSC_USE_COMPLEX)
1204 : r = ctx->rtol;
1205 : PetscCall(PetscOptionsReal("-ds_nep_refine_tol","Refinement tolerance","DSNEPSetRefine",ctx->rtol,&r,&flg1));
1206 : k = ctx->Nit;
1207 : PetscCall(PetscOptionsInt("-ds_nep_refine_its","Number of iterative refinement iterations","DSNEPSetRefine",ctx->Nit,&k,&flg));
1208 : if (flg1||flg) PetscCall(DSNEPSetRefine(ds,r,k));
1209 :
1210 : if (ds->method==1) {
1211 : if (!ctx->rg) PetscCall(DSNEPGetRG(ds,&ctx->rg));
1212 : PetscCall(RGSetFromOptions(ctx->rg));
1213 : }
1214 : #endif
1215 :
1216 19 : PetscOptionsHeadEnd();
1217 19 : PetscFunctionReturn(PETSC_SUCCESS);
1218 : }
1219 :
1220 20 : static PetscErrorCode DSDestroy_NEP(DS ds)
1221 : {
1222 20 : DS_NEP *ctx = (DS_NEP*)ds->data;
1223 20 : PetscInt i;
1224 :
1225 20 : PetscFunctionBegin;
1226 80 : for (i=0;i<ctx->nf;i++) PetscCall(FNDestroy(&ctx->f[i]));
1227 20 : PetscCall(RGDestroy(&ctx->rg));
1228 20 : PetscCall(PetscLayoutDestroy(&ctx->map));
1229 20 : PetscCall(PetscFree(ds->data));
1230 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",NULL));
1231 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",NULL));
1232 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",NULL));
1233 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",NULL));
1234 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",NULL));
1235 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",NULL));
1236 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",NULL));
1237 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",NULL));
1238 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",NULL));
1239 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",NULL));
1240 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",NULL));
1241 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",NULL));
1242 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",NULL));
1243 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",NULL));
1244 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",NULL));
1245 20 : PetscFunctionReturn(PETSC_SUCCESS);
1246 : }
1247 :
1248 4388 : static PetscErrorCode DSMatGetSize_NEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
1249 : {
1250 4388 : DS_NEP *ctx = (DS_NEP*)ds->data;
1251 :
1252 4388 : PetscFunctionBegin;
1253 4388 : *rows = ds->n;
1254 4388 : if (t==DS_MAT_Q || t==DS_MAT_Z || t==DS_MAT_U || t==DS_MAT_V) *rows *= ctx->max_mid;
1255 4388 : *cols = ds->n;
1256 4388 : if (t==DS_MAT_Q || t==DS_MAT_Z || t==DS_MAT_U || t==DS_MAT_V || t==DS_MAT_X || t==DS_MAT_Y) *cols *= ctx->max_mid;
1257 4388 : PetscFunctionReturn(PETSC_SUCCESS);
1258 : }
1259 :
1260 : /*MC
1261 : DSNEP - Dense Nonlinear Eigenvalue Problem.
1262 :
1263 : Level: beginner
1264 :
1265 : Notes:
1266 : The problem is expressed as T(lambda)*x = 0, where T(lambda) is a
1267 : parameter-dependent matrix written as T(lambda) = sum_i E_i*f_i(lambda).
1268 : The eigenvalues lambda are the arguments returned by DSSolve()..
1269 :
1270 : The coefficient matrices E_i are the extra matrices of the DS, and
1271 : the scalar functions f_i are passed via DSNEPSetFN(). Optionally, a
1272 : callback function to fill the E_i matrices can be set with
1273 : DSNEPSetComputeMatrixFunction().
1274 :
1275 : Used DS matrices:
1276 : + DS_MAT_Ex - coefficient matrices of the split form of T(lambda)
1277 : . DS_MAT_X - eigenvectors
1278 : . DS_MAT_A - (workspace) T(lambda) evaluated at a given lambda (SLP only)
1279 : . DS_MAT_B - (workspace) T'(lambda) evaluated at a given lambda (SLP only)
1280 : . DS_MAT_Q - (workspace) left Hankel matrix (contour only)
1281 : . DS_MAT_Z - (workspace) right Hankel matrix (contour only)
1282 : . DS_MAT_U - (workspace) left singular vectors (contour only)
1283 : . DS_MAT_V - (workspace) right singular vectors (contour only)
1284 : - DS_MAT_W - (workspace) auxiliary matrix of size nxn
1285 :
1286 : Implemented methods:
1287 : + 0 - Successive Linear Problems (SLP), computes just one eigenpair
1288 : - 1 - Contour integral, computes all eigenvalues inside a region
1289 :
1290 : .seealso: DSCreate(), DSSetType(), DSType, DSNEPSetFN(), DSNEPSetComputeMatrixFunction()
1291 : M*/
1292 20 : SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS ds)
1293 : {
1294 20 : DS_NEP *ctx;
1295 :
1296 20 : PetscFunctionBegin;
1297 20 : PetscCall(PetscNew(&ctx));
1298 20 : ds->data = (void*)ctx;
1299 20 : ctx->max_mid = 4;
1300 20 : ctx->nnod = 64;
1301 20 : ctx->Nit = 3;
1302 20 : ctx->rtol = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);
1303 :
1304 20 : ds->ops->allocate = DSAllocate_NEP;
1305 20 : ds->ops->setfromoptions = DSSetFromOptions_NEP;
1306 20 : ds->ops->view = DSView_NEP;
1307 20 : ds->ops->vectors = DSVectors_NEP;
1308 20 : ds->ops->solve[0] = DSSolve_NEP_SLP;
1309 : #if defined(PETSC_USE_COMPLEX)
1310 : ds->ops->solve[1] = DSSolve_NEP_Contour;
1311 : #endif
1312 20 : ds->ops->sort = DSSort_NEP;
1313 : #if !defined(PETSC_HAVE_MPIUNI)
1314 20 : ds->ops->synchronize = DSSynchronize_NEP;
1315 : #endif
1316 20 : ds->ops->destroy = DSDestroy_NEP;
1317 20 : ds->ops->matgetsize = DSMatGetSize_NEP;
1318 :
1319 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",DSNEPSetFN_NEP));
1320 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",DSNEPGetFN_NEP));
1321 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",DSNEPGetNumFN_NEP));
1322 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",DSNEPGetMinimality_NEP));
1323 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",DSNEPSetMinimality_NEP));
1324 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",DSNEPGetRefine_NEP));
1325 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",DSNEPSetRefine_NEP));
1326 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",DSNEPGetIntegrationPoints_NEP));
1327 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",DSNEPSetIntegrationPoints_NEP));
1328 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",DSNEPGetSamplingSize_NEP));
1329 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",DSNEPSetSamplingSize_NEP));
1330 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",DSNEPSetRG_NEP));
1331 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",DSNEPGetRG_NEP));
1332 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",DSNEPSetComputeMatrixFunction_NEP));
1333 20 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",DSNEPGetComputeMatrixFunction_NEP));
1334 20 : PetscFunctionReturn(PETSC_SUCCESS);
1335 : }
|
