Actual source code: dsnep.c
slepc-3.9.2 2018-07-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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> /*I "slepcds.h" I*/
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 neig; /* number of available eigenpairs */
18: } DS_NEP;
20: /*
21: DSNEPComputeMatrix - Build the matrix associated with a nonlinear operator
22: T(lambda) or its derivative T'(lambda), given the parameter lambda, where
23: T(lambda) = sum_i E_i*f_i(lambda). The result is written in mat.
24: */
25: static PetscErrorCode DSNEPComputeMatrix(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat)
26: {
28: DS_NEP *ctx = (DS_NEP*)ds->data;
29: PetscScalar *T,*E,alpha;
30: PetscInt i,ld,n;
31: PetscBLASInt k,inc=1;
34: DSGetDimensions(ds,&n,NULL,NULL,NULL,NULL);
35: DSGetLeadingDimension(ds,&ld);
36: PetscBLASIntCast(ld*n,&k);
37: PetscLogEventBegin(DS_Other,ds,0,0,0);
38: DSGetArray(ds,mat,&T);
39: PetscMemzero(T,k*sizeof(PetscScalar));
40: for (i=0;i<ctx->nf;i++) {
41: if (deriv) {
42: FNEvaluateDerivative(ctx->f[i],lambda,&alpha);
43: } else {
44: FNEvaluateFunction(ctx->f[i],lambda,&alpha);
45: }
46: E = ds->mat[DSMatExtra[i]];
47: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k,&alpha,E,&inc,T,&inc));
48: }
49: DSRestoreArray(ds,mat,&T);
50: PetscLogEventEnd(DS_Other,ds,0,0,0);
51: return(0);
52: }
54: PetscErrorCode DSAllocate_NEP(DS ds,PetscInt ld)
55: {
57: DS_NEP *ctx = (DS_NEP*)ds->data;
58: PetscInt i;
61: if (!ctx->nf) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSNEP requires passing some functions via DSSetFN()");
62: DSAllocateMat_Private(ds,DS_MAT_X);
63: for (i=0;i<ctx->nf;i++) {
64: DSAllocateMat_Private(ds,DSMatExtra[i]);
65: }
66: PetscFree(ds->perm);
67: PetscMalloc1(ld,&ds->perm);
68: PetscLogObjectMemory((PetscObject)ds,ld*sizeof(PetscInt));
69: return(0);
70: }
72: PetscErrorCode DSView_NEP(DS ds,PetscViewer viewer)
73: {
74: PetscErrorCode ierr;
75: DS_NEP *ctx = (DS_NEP*)ds->data;
76: PetscViewerFormat format;
77: PetscInt i;
80: PetscViewerGetFormat(viewer,&format);
81: PetscViewerASCIIPrintf(viewer,"number of functions: %D\n",ctx->nf);
82: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
83: for (i=0;i<ctx->nf;i++) {
84: FNView(ctx->f[i],viewer);
85: DSViewMat(ds,viewer,DSMatExtra[i]);
86: }
87: if (ds->state>DS_STATE_INTERMEDIATE) {
88: DSViewMat(ds,viewer,DS_MAT_X);
89: }
90: return(0);
91: }
93: PetscErrorCode DSVectors_NEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
94: {
96: if (rnorm) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
97: switch (mat) {
98: case DS_MAT_X:
99: break;
100: case DS_MAT_Y:
101: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
102: break;
103: default:
104: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
105: }
106: return(0);
107: }
109: PetscErrorCode DSSort_NEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *dummy)
110: {
112: DS_NEP *ctx = (DS_NEP*)ds->data;
113: PetscInt n,l,i,j,k,p,*perm,told,ld=ds->ld;
114: PetscScalar *A,*X,rtmp;
117: if (!ds->sc) return(0);
118: n = ds->n;
119: l = ds->l;
120: A = ds->mat[DS_MAT_A];
121: perm = ds->perm;
122: for (i=0;i<ctx->neig;i++) perm[i] = i;
123: told = ds->t;
124: ds->t = ctx->neig; /* force the sorting routines to consider ctx->neig eigenvalues */
125: if (rr) {
126: DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
127: } else {
128: DSSortEigenvalues_Private(ds,wr,NULL,perm,PETSC_FALSE);
129: }
130: ds->t = told; /* restore value of t */
131: for (i=l;i<n;i++) A[i+i*ld] = wr[perm[i]];
132: for (i=l;i<n;i++) wr[i] = A[i+i*ld];
133: /* cannot use DSPermuteColumns_Private() since not all columns are filled */
134: X = ds->mat[DS_MAT_X];
135: for (i=0;i<ctx->neig;i++) {
136: p = perm[i];
137: if (p != i) {
138: j = i + 1;
139: while (perm[j] != i) j++;
140: perm[j] = p; perm[i] = i;
141: /* swap columns i and j */
142: for (k=0;k<n;k++) {
143: rtmp = X[k+p*ld]; X[k+p*ld] = X[k+i*ld]; X[k+i*ld] = rtmp;
144: }
145: }
146: }
147: return(0);
148: }
150: PetscErrorCode DSSolve_NEP_SLP(DS ds,PetscScalar *wr,PetscScalar *wi)
151: {
152: #if defined(SLEPC_MISSING_LAPACK_GGEV)
154: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GGEV - Lapack routine is unavailable");
155: #else
157: DS_NEP *ctx = (DS_NEP*)ds->data;
158: PetscScalar *A,*B,*W,*X,*work,*alpha,*beta;
159: PetscScalar norm,sigma,lambda,mu,re,re2;
160: PetscBLASInt info,n,ld,lrwork=0,lwork,one=1;
161: PetscInt it,pos,j,maxit=100,result;
162: PetscReal tol;
163: #if defined(PETSC_USE_COMPLEX)
164: PetscReal *rwork;
165: #else
166: PetscReal *alphai,im,im2;
167: #endif
170: if (!ds->mat[DS_MAT_A]) {
171: DSAllocateMat_Private(ds,DS_MAT_A);
172: }
173: if (!ds->mat[DS_MAT_B]) {
174: DSAllocateMat_Private(ds,DS_MAT_B);
175: }
176: if (!ds->mat[DS_MAT_W]) {
177: DSAllocateMat_Private(ds,DS_MAT_W);
178: }
179: PetscBLASIntCast(ds->n,&n);
180: PetscBLASIntCast(ds->ld,&ld);
181: #if defined(PETSC_USE_COMPLEX)
182: PetscBLASIntCast(2*ds->n+2*ds->n,&lwork);
183: PetscBLASIntCast(8*ds->n,&lrwork);
184: #else
185: PetscBLASIntCast(3*ds->n+8*ds->n,&lwork);
186: #endif
187: DSAllocateWork_Private(ds,lwork,lrwork,0);
188: alpha = ds->work;
189: beta = ds->work + ds->n;
190: #if defined(PETSC_USE_COMPLEX)
191: work = ds->work + 2*ds->n;
192: lwork -= 2*ds->n;
193: #else
194: alphai = ds->work + 2*ds->n;
195: work = ds->work + 3*ds->n;
196: lwork -= 3*ds->n;
197: #endif
198: A = ds->mat[DS_MAT_A];
199: B = ds->mat[DS_MAT_B];
200: W = ds->mat[DS_MAT_W];
201: X = ds->mat[DS_MAT_X];
203: sigma = 0.0;
204: if (ds->sc->comparison==SlepcCompareTargetMagnitude || ds->sc->comparison==SlepcCompareTargetReal) sigma = *(PetscScalar*)ds->sc->comparisonctx;
205: lambda = sigma;
206: tol = 1000*n*PETSC_MACHINE_EPSILON;
208: for (it=0;it<maxit;it++) {
210: /* evaluate T and T' */
211: DSNEPComputeMatrix(ds,lambda,PETSC_FALSE,DS_MAT_A);
212: DSNEPComputeMatrix(ds,lambda,PETSC_TRUE,DS_MAT_B);
214: /* compute eigenvalue correction mu and eigenvector u */
215: #if defined(PETSC_USE_COMPLEX)
216: rwork = ds->rwork;
217: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,beta,NULL,&ld,W,&ld,work,&lwork,rwork,&info));
218: #else
219: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,alphai,beta,NULL,&ld,W,&ld,work,&lwork,&info));
220: #endif
221: SlepcCheckLapackInfo("ggev",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: SlepcCompareSmallestMagnitude(re,im,re2,im2,&result,NULL);
239: #else
240: 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: if (im!=0.0) SETERRQ(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: PetscMemcpy(X,W+pos*ld,n*sizeof(PetscScalar));
256: norm = BLASnrm2_(&n,X,&one);
257: norm = 1.0/norm;
258: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X,&one));
260: /* correct eigenvalue approximation */
261: lambda = lambda - mu;
262: if (PetscAbsScalar(mu)<=tol) break;
263: }
265: if (it==maxit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"DSNEP did not converge");
266: ctx->neig = 1;
267: wr[0] = lambda;
268: if (wi) wi[0] = 0.0;
269: return(0);
270: #endif
271: }
273: PetscErrorCode DSSynchronize_NEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
274: {
276: PetscInt k=0;
277: PetscMPIInt n,rank,size,off=0;
280: if (ds->state>=DS_STATE_CONDENSED) k += ds->n;
281: if (eigr) k += 1;
282: if (eigi) k += 1;
283: DSAllocateWork_Private(ds,k,0,0);
284: PetscMPIIntCast(k*sizeof(PetscScalar),&size);
285: PetscMPIIntCast(ds->n,&n);
286: MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
287: if (!rank) {
288: if (ds->state>=DS_STATE_CONDENSED) {
289: MPI_Pack(ds->mat[DS_MAT_X],n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
290: }
291: if (eigr) {
292: MPI_Pack(eigr,1,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
293: }
294: if (eigi) {
295: MPI_Pack(eigi,1,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
296: }
297: }
298: MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
299: if (rank) {
300: if (ds->state>=DS_STATE_CONDENSED) {
301: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_X],n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
302: }
303: if (eigr) {
304: MPI_Unpack(ds->work,size,&off,eigr,1,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
305: }
306: if (eigi) {
307: MPI_Unpack(ds->work,size,&off,eigi,1,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
308: }
309: }
310: return(0);
311: }
313: static PetscErrorCode DSNEPSetFN_NEP(DS ds,PetscInt n,FN fn[])
314: {
316: DS_NEP *ctx = (DS_NEP*)ds->data;
317: PetscInt i;
320: if (n<=0) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more functions, you have %D",n);
321: if (n>DS_NUM_EXTRA) SETERRQ2(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Too many functions, you specified %D but the limit is %D",n,DS_NUM_EXTRA);
322: if (ds->ld) { PetscInfo(ds,"DSNEPSetFN() called after DSAllocate()\n"); }
323: for (i=0;i<ctx->nf;i++) {
324: FNDestroy(&ctx->f[i]);
325: }
326: for (i=0;i<n;i++) {
327: PetscObjectReference((PetscObject)fn[i]);
328: ctx->f[i] = fn[i];
329: }
330: ctx->nf = n;
331: return(0);
332: }
334: /*@
335: DSNEPSetFN - Sets a number of functions that define the nonlinear
336: eigenproblem.
338: Collective on DS and FN
340: Input Parameters:
341: + ds - the direct solver context
342: . n - number of functions
343: - fn - array of functions
345: Notes:
346: The nonlinear eigenproblem is defined in terms of the split nonlinear
347: operator T(lambda) = sum_i A_i*f_i(lambda).
349: This function must be called before DSAllocate(). Then DSAllocate()
350: will allocate an extra matrix A_i per each function, that can be
351: filled in the usual way.
353: Level: advanced
355: .seealso: DSNEPGetFN(), DSAllocate()
356: @*/
357: PetscErrorCode DSNEPSetFN(DS ds,PetscInt n,FN fn[])
358: {
359: PetscInt i;
366: for (i=0;i<n;i++) {
369: }
370: PetscTryMethod(ds,"DSNEPSetFN_C",(DS,PetscInt,FN[]),(ds,n,fn));
371: return(0);
372: }
374: static PetscErrorCode DSNEPGetFN_NEP(DS ds,PetscInt k,FN *fn)
375: {
376: DS_NEP *ctx = (DS_NEP*)ds->data;
379: if (k<0 || k>=ctx->nf) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",ctx->nf-1);
380: *fn = ctx->f[k];
381: return(0);
382: }
384: /*@
385: DSNEPGetFN - Gets the functions associated with the nonlinear DS.
387: Not collective, though parallel FNs are returned if the DS is parallel
389: Input Parameter:
390: + ds - the direct solver context
391: - k - the index of the requested function (starting in 0)
393: Output Parameter:
394: . fn - the function
396: Level: advanced
398: .seealso: DSNEPSetFN()
399: @*/
400: PetscErrorCode DSNEPGetFN(DS ds,PetscInt k,FN *fn)
401: {
407: PetscUseMethod(ds,"DSNEPGetFN_C",(DS,PetscInt,FN*),(ds,k,fn));
408: return(0);
409: }
411: static PetscErrorCode DSNEPGetNumFN_NEP(DS ds,PetscInt *n)
412: {
413: DS_NEP *ctx = (DS_NEP*)ds->data;
416: *n = ctx->nf;
417: return(0);
418: }
420: /*@
421: DSNEPGetNumFN - Returns the number of functions stored internally by
422: the DS.
424: Not collective
426: Input Parameter:
427: . ds - the direct solver context
429: Output Parameters:
430: . n - the number of functions passed in DSNEPSetFN()
432: Level: advanced
434: .seealso: DSNEPSetFN()
435: @*/
436: PetscErrorCode DSNEPGetNumFN(DS ds,PetscInt *n)
437: {
443: PetscUseMethod(ds,"DSNEPGetNumFN_C",(DS,PetscInt*),(ds,n));
444: return(0);
445: }
447: PetscErrorCode DSDestroy_NEP(DS ds)
448: {
450: DS_NEP *ctx = (DS_NEP*)ds->data;
451: PetscInt i;
454: for (i=0;i<ctx->nf;i++) {
455: FNDestroy(&ctx->f[i]);
456: }
457: PetscFree(ds->data);
458: PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",NULL);
459: PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",NULL);
460: PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",NULL);
461: return(0);
462: }
464: PETSC_EXTERN PetscErrorCode DSCreate_NEP(DS ds)
465: {
466: DS_NEP *ctx;
470: PetscNewLog(ds,&ctx);
471: ds->data = (void*)ctx;
473: ds->ops->allocate = DSAllocate_NEP;
474: ds->ops->view = DSView_NEP;
475: ds->ops->vectors = DSVectors_NEP;
476: ds->ops->solve[0] = DSSolve_NEP_SLP;
477: ds->ops->sort = DSSort_NEP;
478: ds->ops->synchronize = DSSynchronize_NEP;
479: ds->ops->destroy = DSDestroy_NEP;
480: PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",DSNEPSetFN_NEP);
481: PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",DSNEPGetFN_NEP);
482: PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",DSNEPGetNumFN_NEP);
483: return(0);
484: }