Actual source code: dsnep.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: void *computematrixctx;
19: PetscErrorCode (*computematrix)(DS,PetscScalar,PetscBool,DSMatType,void*);
20: } DS_NEP;
22: /*
23: DSNEPComputeMatrix - Build the matrix associated with a nonlinear operator
24: T(lambda) or its derivative T'(lambda), given the parameter lambda, where
25: T(lambda) = sum_i E_i*f_i(lambda). The result is written in mat.
26: */
27: static PetscErrorCode DSNEPComputeMatrix(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat)
28: {
30: DS_NEP *ctx = (DS_NEP*)ds->data;
31: PetscScalar *T,*E,alpha;
32: PetscInt i,ld,n;
33: PetscBLASInt k,inc=1;
36: PetscLogEventBegin(DS_Other,ds,0,0,0);
37: if (ctx->computematrix) {
38: (*ctx->computematrix)(ds,lambda,deriv,mat,ctx->computematrixctx);
39: } else {
40: DSGetDimensions(ds,&n,NULL,NULL,NULL,NULL);
41: DSGetLeadingDimension(ds,&ld);
42: PetscBLASIntCast(ld*n,&k);
43: DSGetArray(ds,mat,&T);
44: PetscMemzero(T,k*sizeof(PetscScalar));
45: for (i=0;i<ctx->nf;i++) {
46: if (deriv) {
47: FNEvaluateDerivative(ctx->f[i],lambda,&alpha);
48: } else {
49: FNEvaluateFunction(ctx->f[i],lambda,&alpha);
50: }
51: E = ds->mat[DSMatExtra[i]];
52: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k,&alpha,E,&inc,T,&inc));
53: }
54: DSRestoreArray(ds,mat,&T);
55: }
56: PetscLogEventEnd(DS_Other,ds,0,0,0);
57: return(0);
58: }
60: PetscErrorCode DSAllocate_NEP(DS ds,PetscInt ld)
61: {
63: DS_NEP *ctx = (DS_NEP*)ds->data;
64: PetscInt i;
67: if (!ctx->nf) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSNEP requires passing some functions via DSSetFN()");
68: DSAllocateMat_Private(ds,DS_MAT_X);
69: for (i=0;i<ctx->nf;i++) {
70: DSAllocateMat_Private(ds,DSMatExtra[i]);
71: }
72: PetscFree(ds->perm);
73: PetscMalloc1(ld,&ds->perm);
74: PetscLogObjectMemory((PetscObject)ds,ld*sizeof(PetscInt));
75: return(0);
76: }
78: PetscErrorCode DSView_NEP(DS ds,PetscViewer viewer)
79: {
80: PetscErrorCode ierr;
81: DS_NEP *ctx = (DS_NEP*)ds->data;
82: PetscViewerFormat format;
83: PetscInt i;
86: PetscViewerGetFormat(viewer,&format);
87: PetscViewerASCIIPrintf(viewer,"number of functions: %D\n",ctx->nf);
88: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
89: for (i=0;i<ctx->nf;i++) {
90: FNView(ctx->f[i],viewer);
91: DSViewMat(ds,viewer,DSMatExtra[i]);
92: }
93: if (ds->state>DS_STATE_INTERMEDIATE) {
94: DSViewMat(ds,viewer,DS_MAT_X);
95: }
96: return(0);
97: }
99: PetscErrorCode DSVectors_NEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
100: {
102: if (rnorm) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
103: switch (mat) {
104: case DS_MAT_X:
105: break;
106: case DS_MAT_Y:
107: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
108: break;
109: default:
110: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
111: }
112: return(0);
113: }
115: PetscErrorCode DSSort_NEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *dummy)
116: {
118: DS_NEP *ctx = (DS_NEP*)ds->data;
119: PetscInt n,l,i,j,k,p,*perm,told,ld=ds->ld;
120: PetscScalar *A,*X,rtmp;
123: if (!ds->sc) return(0);
124: n = ds->n;
125: l = ds->l;
126: A = ds->mat[DS_MAT_A];
127: perm = ds->perm;
128: for (i=0;i<ctx->neig;i++) perm[i] = i;
129: told = ds->t;
130: ds->t = ctx->neig; /* force the sorting routines to consider ctx->neig eigenvalues */
131: if (rr) {
132: DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
133: } else {
134: DSSortEigenvalues_Private(ds,wr,NULL,perm,PETSC_FALSE);
135: }
136: ds->t = told; /* restore value of t */
137: for (i=l;i<n;i++) A[i+i*ld] = wr[perm[i]];
138: for (i=l;i<n;i++) wr[i] = A[i+i*ld];
139: /* cannot use DSPermuteColumns_Private() since not all columns are filled */
140: X = ds->mat[DS_MAT_X];
141: for (i=0;i<ctx->neig;i++) {
142: p = perm[i];
143: if (p != i) {
144: j = i + 1;
145: while (perm[j] != i) j++;
146: perm[j] = p; perm[i] = i;
147: /* swap columns i and j */
148: for (k=0;k<n;k++) {
149: rtmp = X[k+p*ld]; X[k+p*ld] = X[k+i*ld]; X[k+i*ld] = rtmp;
150: }
151: }
152: }
153: return(0);
154: }
156: PetscErrorCode DSSolve_NEP_SLP(DS ds,PetscScalar *wr,PetscScalar *wi)
157: {
158: #if defined(SLEPC_MISSING_LAPACK_GGEV)
160: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GGEV - Lapack routine is unavailable");
161: #else
163: DS_NEP *ctx = (DS_NEP*)ds->data;
164: PetscScalar *A,*B,*W,*X,*work,*alpha,*beta;
165: PetscScalar norm,sigma,lambda,mu,re,re2,sone=1.0,zero=0.0;
166: PetscBLASInt info,n,ld,lrwork=0,lwork,one=1;
167: PetscInt it,pos,j,maxit=100,result;
168: PetscReal tol;
169: #if defined(PETSC_USE_COMPLEX)
170: PetscReal *rwork;
171: #else
172: PetscReal *alphai,im,im2;
173: #endif
176: if (!ds->mat[DS_MAT_A]) {
177: DSAllocateMat_Private(ds,DS_MAT_A);
178: }
179: if (!ds->mat[DS_MAT_B]) {
180: DSAllocateMat_Private(ds,DS_MAT_B);
181: }
182: if (!ds->mat[DS_MAT_W]) {
183: DSAllocateMat_Private(ds,DS_MAT_W);
184: }
185: PetscBLASIntCast(ds->n,&n);
186: PetscBLASIntCast(ds->ld,&ld);
187: #if defined(PETSC_USE_COMPLEX)
188: PetscBLASIntCast(2*ds->n+2*ds->n,&lwork);
189: PetscBLASIntCast(8*ds->n,&lrwork);
190: #else
191: PetscBLASIntCast(3*ds->n+8*ds->n,&lwork);
192: #endif
193: DSAllocateWork_Private(ds,lwork,lrwork,0);
194: alpha = ds->work;
195: beta = ds->work + ds->n;
196: #if defined(PETSC_USE_COMPLEX)
197: work = ds->work + 2*ds->n;
198: lwork -= 2*ds->n;
199: #else
200: alphai = ds->work + 2*ds->n;
201: work = ds->work + 3*ds->n;
202: lwork -= 3*ds->n;
203: #endif
204: A = ds->mat[DS_MAT_A];
205: B = ds->mat[DS_MAT_B];
206: W = ds->mat[DS_MAT_W];
207: X = ds->mat[DS_MAT_X];
209: sigma = 0.0;
210: if (ds->sc->comparison==SlepcCompareTargetMagnitude || ds->sc->comparison==SlepcCompareTargetReal) sigma = *(PetscScalar*)ds->sc->comparisonctx;
211: lambda = sigma;
212: tol = 1000*n*PETSC_MACHINE_EPSILON;
214: for (it=0;it<maxit;it++) {
216: /* evaluate T and T' */
217: DSNEPComputeMatrix(ds,lambda,PETSC_FALSE,DS_MAT_A);
218: if (it) {
219: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,A,&ld,X,&one,&zero,X+ld,&one));
220: norm = BLASnrm2_(&n,X+ld,&one);
221: if (PetscRealPart(norm)/PetscAbsScalar(lambda)<=tol) break;
222: }
223: DSNEPComputeMatrix(ds,lambda,PETSC_TRUE,DS_MAT_B);
225: /* compute eigenvalue correction mu and eigenvector u */
226: #if defined(PETSC_USE_COMPLEX)
227: rwork = ds->rwork;
228: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,beta,NULL,&ld,W,&ld,work,&lwork,rwork,&info));
229: #else
230: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,alphai,beta,NULL,&ld,W,&ld,work,&lwork,&info));
231: #endif
232: SlepcCheckLapackInfo("ggev",info);
234: /* find smallest eigenvalue */
235: j = 0;
236: if (beta[j]==0.0) re = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
237: else re = alpha[j]/beta[j];
238: #if !defined(PETSC_USE_COMPLEX)
239: if (beta[j]==0.0) im = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
240: else im = alphai[j]/beta[j];
241: #endif
242: pos = 0;
243: for (j=1;j<n;j++) {
244: if (beta[j]==0.0) re2 = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
245: else re2 = alpha[j]/beta[j];
246: #if !defined(PETSC_USE_COMPLEX)
247: if (beta[j]==0.0) im2 = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
248: else im2 = alphai[j]/beta[j];
249: SlepcCompareSmallestMagnitude(re,im,re2,im2,&result,NULL);
250: #else
251: SlepcCompareSmallestMagnitude(re,0.0,re2,0.0,&result,NULL);
252: #endif
253: if (result > 0) {
254: re = re2;
255: #if !defined(PETSC_USE_COMPLEX)
256: im = im2;
257: #endif
258: pos = j;
259: }
260: }
262: #if !defined(PETSC_USE_COMPLEX)
263: if (im!=0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DSNEP found a complex eigenvalue; try rerunning with complex scalars");
264: #endif
265: mu = alpha[pos]/beta[pos];
266: PetscMemcpy(X,W+pos*ld,n*sizeof(PetscScalar));
267: norm = BLASnrm2_(&n,X,&one);
268: norm = 1.0/norm;
269: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X,&one));
271: /* correct eigenvalue approximation */
272: lambda = lambda - mu;
273: }
275: if (it==maxit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"DSNEP did not converge");
276: ctx->neig = 1;
277: wr[0] = lambda;
278: if (wi) wi[0] = 0.0;
279: return(0);
280: #endif
281: }
283: PetscErrorCode DSSynchronize_NEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
284: {
286: PetscInt k=0;
287: PetscMPIInt n,rank,size,off=0;
290: if (ds->state>=DS_STATE_CONDENSED) k += ds->n;
291: if (eigr) k += 1;
292: if (eigi) k += 1;
293: DSAllocateWork_Private(ds,k,0,0);
294: PetscMPIIntCast(k*sizeof(PetscScalar),&size);
295: PetscMPIIntCast(ds->n,&n);
296: MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
297: if (!rank) {
298: if (ds->state>=DS_STATE_CONDENSED) {
299: MPI_Pack(ds->mat[DS_MAT_X],n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
300: }
301: if (eigr) {
302: MPI_Pack(eigr,1,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
303: }
304: if (eigi) {
305: MPI_Pack(eigi,1,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
306: }
307: }
308: MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
309: if (rank) {
310: if (ds->state>=DS_STATE_CONDENSED) {
311: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_X],n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
312: }
313: if (eigr) {
314: MPI_Unpack(ds->work,size,&off,eigr,1,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
315: }
316: if (eigi) {
317: MPI_Unpack(ds->work,size,&off,eigi,1,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
318: }
319: }
320: return(0);
321: }
323: static PetscErrorCode DSNEPSetFN_NEP(DS ds,PetscInt n,FN fn[])
324: {
326: DS_NEP *ctx = (DS_NEP*)ds->data;
327: PetscInt i;
330: if (n<=0) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more functions, you have %D",n);
331: 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);
332: if (ds->ld) { PetscInfo(ds,"DSNEPSetFN() called after DSAllocate()\n"); }
333: for (i=0;i<n;i++) {
334: PetscObjectReference((PetscObject)fn[i]);
335: }
336: for (i=0;i<ctx->nf;i++) {
337: FNDestroy(&ctx->f[i]);
338: }
339: for (i=0;i<n;i++) ctx->f[i] = fn[i];
340: ctx->nf = n;
341: return(0);
342: }
344: /*@
345: DSNEPSetFN - Sets a number of functions that define the nonlinear
346: eigenproblem.
348: Collective on DS and FN
350: Input Parameters:
351: + ds - the direct solver context
352: . n - number of functions
353: - fn - array of functions
355: Notes:
356: The nonlinear eigenproblem is defined in terms of the split nonlinear
357: operator T(lambda) = sum_i A_i*f_i(lambda).
359: This function must be called before DSAllocate(). Then DSAllocate()
360: will allocate an extra matrix A_i per each function, that can be
361: filled in the usual way.
363: Level: advanced
365: .seealso: DSNEPGetFN(), DSAllocate()
366: @*/
367: PetscErrorCode DSNEPSetFN(DS ds,PetscInt n,FN fn[])
368: {
369: PetscInt i;
376: for (i=0;i<n;i++) {
379: }
380: PetscTryMethod(ds,"DSNEPSetFN_C",(DS,PetscInt,FN[]),(ds,n,fn));
381: return(0);
382: }
384: static PetscErrorCode DSNEPGetFN_NEP(DS ds,PetscInt k,FN *fn)
385: {
386: DS_NEP *ctx = (DS_NEP*)ds->data;
389: if (k<0 || k>=ctx->nf) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",ctx->nf-1);
390: *fn = ctx->f[k];
391: return(0);
392: }
394: /*@
395: DSNEPGetFN - Gets the functions associated with the nonlinear DS.
397: Not collective, though parallel FNs are returned if the DS is parallel
399: Input Parameter:
400: + ds - the direct solver context
401: - k - the index of the requested function (starting in 0)
403: Output Parameter:
404: . fn - the function
406: Level: advanced
408: .seealso: DSNEPSetFN()
409: @*/
410: PetscErrorCode DSNEPGetFN(DS ds,PetscInt k,FN *fn)
411: {
417: PetscUseMethod(ds,"DSNEPGetFN_C",(DS,PetscInt,FN*),(ds,k,fn));
418: return(0);
419: }
421: static PetscErrorCode DSNEPGetNumFN_NEP(DS ds,PetscInt *n)
422: {
423: DS_NEP *ctx = (DS_NEP*)ds->data;
426: *n = ctx->nf;
427: return(0);
428: }
430: /*@
431: DSNEPGetNumFN - Returns the number of functions stored internally by
432: the DS.
434: Not collective
436: Input Parameter:
437: . ds - the direct solver context
439: Output Parameters:
440: . n - the number of functions passed in DSNEPSetFN()
442: Level: advanced
444: .seealso: DSNEPSetFN()
445: @*/
446: PetscErrorCode DSNEPGetNumFN(DS ds,PetscInt *n)
447: {
453: PetscUseMethod(ds,"DSNEPGetNumFN_C",(DS,PetscInt*),(ds,n));
454: return(0);
455: }
457: static PetscErrorCode DSNEPSetComputeMatrixFunction_NEP(DS ds,PetscErrorCode (*fun)(DS,PetscScalar,PetscBool,DSMatType,void*),void *ctx)
458: {
459: DS_NEP *dsctx = (DS_NEP*)ds->data;
462: dsctx->computematrix = fun;
463: dsctx->computematrixctx = ctx;
464: return(0);
465: }
467: /*@
468: DSNEPSetComputeMatrixFunction - Sets a user-provided subroutine to compute
469: the matrices T(lambda) or T'(lambda).
471: Logically Collective on DS
473: Input Parameters:
474: + ds - the direct solver context
475: . fun - a pointer to the user function
476: - ctx - a context pointer (the last parameter to the user function)
478: Calling Sequence of fun:
479: $ fun(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat,void *ctx)
481: + ds - the direct solver object
482: . lambda - point where T(lambda) or T'(lambda) must be evaluated
483: . deriv - if true compute T'(lambda), otherwise compute T(lambda)
484: . mat - the DS matrix where the result must be stored
485: - ctx - optional context, as set by DSNEPSetComputeMatrixFunction()
487: Note:
488: The result is computed as T(lambda) = sum_i E_i*f_i(lambda), and similarly
489: for the derivative.
491: Level: developer
492: @*/
493: PetscErrorCode DSNEPSetComputeMatrixFunction(DS ds,PetscErrorCode (*fun)(DS,PetscScalar,PetscBool,DSMatType,void*),void *ctx)
494: {
499: PetscTryMethod(ds,"DSNEPSetComputeMatrixFunction_C",(DS,PetscErrorCode (*)(DS,PetscScalar,PetscBool,DSMatType,void*),void*),(ds,fun,ctx));
500: return(0);
501: }
503: PetscErrorCode DSDestroy_NEP(DS ds)
504: {
506: DS_NEP *ctx = (DS_NEP*)ds->data;
507: PetscInt i;
510: for (i=0;i<ctx->nf;i++) {
511: FNDestroy(&ctx->f[i]);
512: }
513: PetscFree(ds->data);
514: PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",NULL);
515: PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",NULL);
516: PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",NULL);
517: PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",NULL);
518: return(0);
519: }
521: SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS ds)
522: {
523: DS_NEP *ctx;
527: PetscNewLog(ds,&ctx);
528: ds->data = (void*)ctx;
530: ds->ops->allocate = DSAllocate_NEP;
531: ds->ops->view = DSView_NEP;
532: ds->ops->vectors = DSVectors_NEP;
533: ds->ops->solve[0] = DSSolve_NEP_SLP;
534: ds->ops->sort = DSSort_NEP;
535: ds->ops->synchronize = DSSynchronize_NEP;
536: ds->ops->destroy = DSDestroy_NEP;
537: PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",DSNEPSetFN_NEP);
538: PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",DSNEPGetFN_NEP);
539: PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",DSNEPGetNumFN_NEP);
540: PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",DSNEPSetComputeMatrixFunction_NEP);
541: return(0);
542: }