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