Actual source code: fnrational.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: */
10: /*
11: Rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials
12: */
14: #include <slepc/private/fnimpl.h> /*I "slepcfn.h" I*/
15: #include <slepcblaslapack.h>
17: typedef struct {
18: PetscScalar *pcoeff; /* numerator coefficients */
19: PetscInt np; /* length of array pcoeff, p(x) has degree np-1 */
20: PetscScalar *qcoeff; /* denominator coefficients */
21: PetscInt nq; /* length of array qcoeff, q(x) has degree nq-1 */
22: } FN_RATIONAL;
24: PetscErrorCode FNEvaluateFunction_Rational(FN fn,PetscScalar x,PetscScalar *y)
25: {
26: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
27: PetscInt i;
28: PetscScalar p,q;
31: if (!ctx->np) p = 1.0;
32: else {
33: p = ctx->pcoeff[0];
34: for (i=1;i<ctx->np;i++)
35: p = ctx->pcoeff[i]+x*p;
36: }
37: if (!ctx->nq) *y = p;
38: else {
39: q = ctx->qcoeff[0];
40: for (i=1;i<ctx->nq;i++)
41: q = ctx->qcoeff[i]+x*q;
42: if (q==0.0) SETERRQ(PETSC_COMM_SELF,1,"Function not defined in the requested value");
43: *y = p/q;
44: }
45: return(0);
46: }
48: static PetscErrorCode FNEvaluateFunctionMat_Rational_Private(FN fn,PetscScalar *Aa,PetscScalar *Ba,PetscInt m,PetscBool firstonly)
49: {
50: #if defined(PETSC_MISSING_LAPACK_GESV)
52: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESV - Lapack routines are unavailable");
53: #else
55: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
56: PetscBLASInt n,k,ld,*ipiv,info;
57: PetscInt i,j;
58: PetscScalar *W,*P,*Q,one=1.0,zero=0.0;
61: PetscBLASIntCast(m,&n);
62: ld = n;
63: k = firstonly? 1: n;
64: if (Aa==Ba) {
65: PetscMalloc4(m*m,&P,m*m,&Q,m*m,&W,ld,&ipiv);
66: } else {
67: P = Ba;
68: PetscMalloc3(m*m,&Q,m*m,&W,ld,&ipiv);
69: }
70: PetscMemzero(P,m*m*sizeof(PetscScalar));
71: if (!ctx->np) {
72: for (i=0;i<m;i++) P[i+i*ld] = 1.0;
73: } else {
74: for (i=0;i<m;i++) P[i+i*ld] = ctx->pcoeff[0];
75: for (j=1;j<ctx->np;j++) {
76: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,P,&ld,Aa,&ld,&zero,W,&ld));
77: PetscMemcpy(P,W,m*m*sizeof(PetscScalar));
78: for (i=0;i<m;i++) P[i+i*ld] += ctx->pcoeff[j];
79: }
80: PetscLogFlops(2.0*n*n*n*(ctx->np-1));
81: }
82: if (ctx->nq) {
83: PetscMemzero(Q,m*m*sizeof(PetscScalar));
84: for (i=0;i<m;i++) Q[i+i*ld] = ctx->qcoeff[0];
85: for (j=1;j<ctx->nq;j++) {
86: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,Q,&ld,Aa,&ld,&zero,W,&ld));
87: PetscMemcpy(Q,W,m*m*sizeof(PetscScalar));
88: for (i=0;i<m;i++) Q[i+i*ld] += ctx->qcoeff[j];
89: }
90: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&k,Q,&ld,ipiv,P,&ld,&info));
91: SlepcCheckLapackInfo("gesv",info);
92: PetscLogFlops(2.0*n*n*n*(ctx->nq-1)+2.0*n*n*n/3.0+2.0*n*n*k);
93: }
94: if (Aa==Ba) {
95: PetscMemcpy(Aa,P,m*k*sizeof(PetscScalar));
96: PetscFree4(P,Q,W,ipiv);
97: } else {
98: PetscFree3(Q,W,ipiv);
99: }
100: return(0);
101: #endif
102: }
104: PetscErrorCode FNEvaluateFunctionMat_Rational(FN fn,Mat A,Mat B)
105: {
107: PetscInt m;
108: PetscScalar *Aa,*Ba;
111: MatDenseGetArray(A,&Aa);
112: MatDenseGetArray(B,&Ba);
113: MatGetSize(A,&m,NULL);
114: FNEvaluateFunctionMat_Rational_Private(fn,Aa,Ba,m,PETSC_FALSE);
115: MatDenseRestoreArray(A,&Aa);
116: MatDenseRestoreArray(B,&Ba);
117: return(0);
118: }
120: PetscErrorCode FNEvaluateFunctionMatVec_Rational(FN fn,Mat A,Vec v)
121: {
123: PetscInt m;
124: PetscScalar *Aa,*Ba;
125: Mat B;
128: FN_AllocateWorkMat(fn,A,&B);
129: MatDenseGetArray(A,&Aa);
130: MatDenseGetArray(B,&Ba);
131: MatGetSize(A,&m,NULL);
132: FNEvaluateFunctionMat_Rational_Private(fn,Aa,Ba,m,PETSC_TRUE);
133: MatDenseRestoreArray(A,&Aa);
134: MatDenseRestoreArray(B,&Ba);
135: MatGetColumnVector(B,v,0);
136: FN_FreeWorkMat(fn,&B);
137: return(0);
138: }
140: PetscErrorCode FNEvaluateDerivative_Rational(FN fn,PetscScalar x,PetscScalar *yp)
141: {
142: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
143: PetscInt i;
144: PetscScalar p,q,pp,qp;
147: if (!ctx->np) {
148: p = 1.0;
149: pp = 0.0;
150: } else {
151: p = ctx->pcoeff[0];
152: pp = 0.0;
153: for (i=1;i<ctx->np;i++) {
154: pp = p+x*pp;
155: p = ctx->pcoeff[i]+x*p;
156: }
157: }
158: if (!ctx->nq) *yp = pp;
159: else {
160: q = ctx->qcoeff[0];
161: qp = 0.0;
162: for (i=1;i<ctx->nq;i++) {
163: qp = q+x*qp;
164: q = ctx->qcoeff[i]+x*q;
165: }
166: if (q==0.0) SETERRQ(PETSC_COMM_SELF,1,"Derivative not defined in the requested value");
167: *yp = (pp*q-p*qp)/(q*q);
168: }
169: return(0);
170: }
172: PetscErrorCode FNView_Rational(FN fn,PetscViewer viewer)
173: {
175: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
176: PetscBool isascii;
177: PetscInt i;
178: char str[50];
181: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
182: if (isascii) {
183: if (fn->alpha!=(PetscScalar)1.0 || fn->beta!=(PetscScalar)1.0) {
184: SlepcSNPrintfScalar(str,50,fn->alpha,PETSC_FALSE);
185: PetscViewerASCIIPrintf(viewer," Scale factors: alpha=%s,",str);
186: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
187: SlepcSNPrintfScalar(str,50,fn->beta,PETSC_FALSE);
188: PetscViewerASCIIPrintf(viewer," beta=%s\n",str);
189: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
190: }
191: if (!ctx->nq) {
192: if (!ctx->np) {
193: PetscViewerASCIIPrintf(viewer," Constant: 1.0\n");
194: } else if (ctx->np==1) {
195: SlepcSNPrintfScalar(str,50,ctx->pcoeff[0],PETSC_FALSE);
196: PetscViewerASCIIPrintf(viewer," Constant: %s\n",str);
197: } else {
198: PetscViewerASCIIPrintf(viewer," Polynomial: ");
199: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
200: for (i=0;i<ctx->np-1;i++) {
201: SlepcSNPrintfScalar(str,50,ctx->pcoeff[i],PETSC_TRUE);
202: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->np-i-1);
203: }
204: SlepcSNPrintfScalar(str,50,ctx->pcoeff[ctx->np-1],PETSC_TRUE);
205: PetscViewerASCIIPrintf(viewer,"%s\n",str);
206: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
207: }
208: } else if (!ctx->np) {
209: PetscViewerASCIIPrintf(viewer," Inverse polinomial: 1 / (");
210: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
211: for (i=0;i<ctx->nq-1;i++) {
212: SlepcSNPrintfScalar(str,50,ctx->qcoeff[i],PETSC_TRUE);
213: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->nq-i-1);
214: }
215: SlepcSNPrintfScalar(str,50,ctx->qcoeff[ctx->nq-1],PETSC_TRUE);
216: PetscViewerASCIIPrintf(viewer,"%s)\n",str);
217: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
218: } else {
219: PetscViewerASCIIPrintf(viewer," Rational function: (");
220: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
221: for (i=0;i<ctx->np-1;i++) {
222: SlepcSNPrintfScalar(str,50,ctx->pcoeff[i],PETSC_TRUE);
223: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->np-i-1);
224: }
225: SlepcSNPrintfScalar(str,50,ctx->pcoeff[ctx->np-1],PETSC_TRUE);
226: PetscViewerASCIIPrintf(viewer,"%s) / (",str);
227: for (i=0;i<ctx->nq-1;i++) {
228: SlepcSNPrintfScalar(str,50,ctx->qcoeff[i],PETSC_TRUE);
229: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->nq-i-1);
230: }
231: SlepcSNPrintfScalar(str,50,ctx->qcoeff[ctx->nq-1],PETSC_TRUE);
232: PetscViewerASCIIPrintf(viewer,"%s)\n",str);
233: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
234: }
235: }
236: return(0);
237: }
239: static PetscErrorCode FNRationalSetNumerator_Rational(FN fn,PetscInt np,PetscScalar *pcoeff)
240: {
242: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
243: PetscInt i;
246: if (np<0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Argument np cannot be negative");
247: ctx->np = np;
248: PetscFree(ctx->pcoeff);
249: if (np) {
250: PetscMalloc1(np,&ctx->pcoeff);
251: PetscLogObjectMemory((PetscObject)fn,np*sizeof(PetscScalar));
252: for (i=0;i<np;i++) ctx->pcoeff[i] = pcoeff[i];
253: }
254: return(0);
255: }
257: /*@C
258: FNRationalSetNumerator - Sets the parameters defining the numerator of the
259: rational function.
261: Logically Collective on FN
263: Input Parameters:
264: + fn - the math function context
265: . np - number of coefficients
266: - pcoeff - coefficients (array of scalar values)
268: Notes:
269: Let the rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials.
270: This function provides the coefficients of the numerator p(x).
271: Hence, p(x) is of degree np-1.
272: If np is zero, then the numerator is assumed to be p(x)=1.
274: In polynomials, high order coefficients are stored in the first positions
275: of the array, e.g. to represent x^2-3 use {1,0,-3}.
277: Level: intermediate
279: .seealso: FNRationalSetDenominator(), FNRationalGetNumerator()
280: @*/
281: PetscErrorCode FNRationalSetNumerator(FN fn,PetscInt np,PetscScalar *pcoeff)
282: {
289: PetscTryMethod(fn,"FNRationalSetNumerator_C",(FN,PetscInt,PetscScalar*),(fn,np,pcoeff));
290: return(0);
291: }
293: static PetscErrorCode FNRationalGetNumerator_Rational(FN fn,PetscInt *np,PetscScalar *pcoeff[])
294: {
296: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
297: PetscInt i;
300: if (np) *np = ctx->np;
301: if (pcoeff) {
302: if (!ctx->np) *pcoeff = NULL;
303: else {
304: PetscMalloc1(ctx->np,pcoeff);
305: for (i=0;i<ctx->np;i++) (*pcoeff)[i] = ctx->pcoeff[i];
306: }
307: }
308: return(0);
309: }
311: /*@C
312: FNRationalGetNumerator - Gets the parameters that define the numerator of the
313: rational function.
315: Not Collective
317: Input Parameter:
318: . fn - the math function context
320: Output Parameters:
321: + np - number of coefficients
322: - pcoeff - coefficients (array of scalar values, length nq)
324: Notes:
325: The values passed by user with FNRationalSetNumerator() are returned (or null
326: pointers otherwise).
327: The pcoeff array should be freed by the user when no longer needed.
329: Level: intermediate
331: .seealso: FNRationalSetNumerator()
332: @*/
333: PetscErrorCode FNRationalGetNumerator(FN fn,PetscInt *np,PetscScalar *pcoeff[])
334: {
339: PetscUseMethod(fn,"FNRationalGetNumerator_C",(FN,PetscInt*,PetscScalar**),(fn,np,pcoeff));
340: return(0);
341: }
343: static PetscErrorCode FNRationalSetDenominator_Rational(FN fn,PetscInt nq,PetscScalar *qcoeff)
344: {
346: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
347: PetscInt i;
350: if (nq<0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Argument nq cannot be negative");
351: ctx->nq = nq;
352: PetscFree(ctx->qcoeff);
353: if (nq) {
354: PetscMalloc1(nq,&ctx->qcoeff);
355: PetscLogObjectMemory((PetscObject)fn,nq*sizeof(PetscScalar));
356: for (i=0;i<nq;i++) ctx->qcoeff[i] = qcoeff[i];
357: }
358: return(0);
359: }
361: /*@C
362: FNRationalSetDenominator - Sets the parameters defining the denominator of the
363: rational function.
365: Logically Collective on FN
367: Input Parameters:
368: + fn - the math function context
369: . nq - number of coefficients
370: - qcoeff - coefficients (array of scalar values)
372: Notes:
373: Let the rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials.
374: This function provides the coefficients of the denominator q(x).
375: Hence, q(x) is of degree nq-1.
376: If nq is zero, then the function is assumed to be polynomial, r(x) = p(x).
378: In polynomials, high order coefficients are stored in the first positions
379: of the array, e.g. to represent x^2-3 use {1,0,-3}.
381: Level: intermediate
383: .seealso: FNRationalSetNumerator(), FNRationalGetDenominator()
384: @*/
385: PetscErrorCode FNRationalSetDenominator(FN fn,PetscInt nq,PetscScalar *qcoeff)
386: {
393: PetscTryMethod(fn,"FNRationalSetDenominator_C",(FN,PetscInt,PetscScalar*),(fn,nq,qcoeff));
394: return(0);
395: }
397: static PetscErrorCode FNRationalGetDenominator_Rational(FN fn,PetscInt *nq,PetscScalar *qcoeff[])
398: {
400: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
401: PetscInt i;
404: if (nq) *nq = ctx->nq;
405: if (qcoeff) {
406: if (!ctx->nq) *qcoeff = NULL;
407: else {
408: PetscMalloc1(ctx->nq,qcoeff);
409: for (i=0;i<ctx->nq;i++) (*qcoeff)[i] = ctx->qcoeff[i];
410: }
411: }
412: return(0);
413: }
415: /*@C
416: FNRationalGetDenominator - Gets the parameters that define the denominator of the
417: rational function.
419: Not Collective
421: Input Parameter:
422: . fn - the math function context
424: Output Parameters:
425: + nq - number of coefficients
426: - qcoeff - coefficients (array of scalar values, length nq)
428: Notes:
429: The values passed by user with FNRationalSetDenominator() are returned (or null
430: pointers otherwise).
431: The qcoeff array should be freed by the user when no longer needed.
433: Level: intermediate
435: .seealso: FNRationalSetDenominator()
436: @*/
437: PetscErrorCode FNRationalGetDenominator(FN fn,PetscInt *nq,PetscScalar *qcoeff[])
438: {
443: PetscUseMethod(fn,"FNRationalGetDenominator_C",(FN,PetscInt*,PetscScalar**),(fn,nq,qcoeff));
444: return(0);
445: }
447: PetscErrorCode FNSetFromOptions_Rational(PetscOptionItems *PetscOptionsObject,FN fn)
448: {
450: #define PARMAX 10
451: PetscScalar array[PARMAX];
452: PetscInt i,k;
453: PetscBool flg;
456: PetscOptionsHead(PetscOptionsObject,"FN Rational Options");
458: k = PARMAX;
459: for (i=0;i<k;i++) array[i] = 0;
460: PetscOptionsScalarArray("-fn_rational_numerator","Numerator coefficients (one or more scalar values separated with a comma without spaces)","FNRationalSetNumerator",array,&k,&flg);
461: if (flg) { FNRationalSetNumerator(fn,k,array); }
463: k = PARMAX;
464: for (i=0;i<k;i++) array[i] = 0;
465: PetscOptionsScalarArray("-fn_rational_denominator","Denominator coefficients (one or more scalar values separated with a comma without spaces)","FNRationalSetDenominator",array,&k,&flg);
466: if (flg) { FNRationalSetDenominator(fn,k,array); }
468: PetscOptionsTail();
469: return(0);
470: }
472: PetscErrorCode FNDuplicate_Rational(FN fn,MPI_Comm comm,FN *newfn)
473: {
475: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data,*ctx2 = (FN_RATIONAL*)(*newfn)->data;
476: PetscInt i;
479: ctx2->np = ctx->np;
480: if (ctx->np) {
481: PetscMalloc1(ctx->np,&ctx2->pcoeff);
482: PetscLogObjectMemory((PetscObject)(*newfn),ctx->np*sizeof(PetscScalar));
483: for (i=0;i<ctx->np;i++) ctx2->pcoeff[i] = ctx->pcoeff[i];
484: }
485: ctx2->nq = ctx->nq;
486: if (ctx->nq) {
487: PetscMalloc1(ctx->nq,&ctx2->qcoeff);
488: PetscLogObjectMemory((PetscObject)(*newfn),ctx->nq*sizeof(PetscScalar));
489: for (i=0;i<ctx->nq;i++) ctx2->qcoeff[i] = ctx->qcoeff[i];
490: }
491: return(0);
492: }
494: PetscErrorCode FNDestroy_Rational(FN fn)
495: {
497: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
500: PetscFree(ctx->pcoeff);
501: PetscFree(ctx->qcoeff);
502: PetscFree(fn->data);
503: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetNumerator_C",NULL);
504: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetNumerator_C",NULL);
505: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetDenominator_C",NULL);
506: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetDenominator_C",NULL);
507: return(0);
508: }
510: SLEPC_EXTERN PetscErrorCode FNCreate_Rational(FN fn)
511: {
513: FN_RATIONAL *ctx;
516: PetscNewLog(fn,&ctx);
517: fn->data = (void*)ctx;
519: fn->ops->evaluatefunction = FNEvaluateFunction_Rational;
520: fn->ops->evaluatederivative = FNEvaluateDerivative_Rational;
521: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Rational;
522: fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Rational;
523: fn->ops->setfromoptions = FNSetFromOptions_Rational;
524: fn->ops->view = FNView_Rational;
525: fn->ops->duplicate = FNDuplicate_Rational;
526: fn->ops->destroy = FNDestroy_Rational;
527: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetNumerator_C",FNRationalSetNumerator_Rational);
528: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetNumerator_C",FNRationalGetNumerator_Rational);
529: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetDenominator_C",FNRationalSetDenominator_Rational);
530: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetDenominator_C",FNRationalGetDenominator_Rational);
531: return(0);
532: }