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