Actual source code: fnrational.c
slepc-3.16.2 2022-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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>
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,PETSC_ERR_ARG_OUTOFRANGE,"Function not defined in the requested value");
43: *y = p/q;
44: }
45: return(0);
46: }
48: static PetscErrorCode FNEvaluateFunctionMat_Rational_Private(FN fn,const PetscScalar *Aa,PetscScalar *Ba,PetscInt m,PetscBool firstonly)
49: {
51: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
52: PetscBLASInt n,k,ld,*ipiv,info;
53: PetscInt i,j;
54: PetscScalar *W,*P,*Q,one=1.0,zero=0.0;
57: PetscBLASIntCast(m,&n);
58: ld = n;
59: k = firstonly? 1: n;
60: if (Aa==Ba) {
61: PetscMalloc4(m*m,&P,m*m,&Q,m*m,&W,ld,&ipiv);
62: } else {
63: P = Ba;
64: PetscMalloc3(m*m,&Q,m*m,&W,ld,&ipiv);
65: }
66: PetscArrayzero(P,m*m);
67: if (!ctx->np) {
68: for (i=0;i<m;i++) P[i+i*ld] = 1.0;
69: } else {
70: for (i=0;i<m;i++) P[i+i*ld] = ctx->pcoeff[0];
71: for (j=1;j<ctx->np;j++) {
72: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,P,&ld,Aa,&ld,&zero,W,&ld));
73: PetscArraycpy(P,W,m*m);
74: for (i=0;i<m;i++) P[i+i*ld] += ctx->pcoeff[j];
75: }
76: PetscLogFlops(2.0*n*n*n*(ctx->np-1));
77: }
78: if (ctx->nq) {
79: PetscArrayzero(Q,m*m);
80: for (i=0;i<m;i++) Q[i+i*ld] = ctx->qcoeff[0];
81: for (j=1;j<ctx->nq;j++) {
82: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,Q,&ld,Aa,&ld,&zero,W,&ld));
83: PetscArraycpy(Q,W,m*m);
84: for (i=0;i<m;i++) Q[i+i*ld] += ctx->qcoeff[j];
85: }
86: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&k,Q,&ld,ipiv,P,&ld,&info));
87: SlepcCheckLapackInfo("gesv",info);
88: PetscLogFlops(2.0*n*n*n*(ctx->nq-1)+2.0*n*n*n/3.0+2.0*n*n*k);
89: }
90: if (Aa==Ba) {
91: PetscArraycpy(Ba,P,m*k);
92: PetscFree4(P,Q,W,ipiv);
93: } else {
94: PetscFree3(Q,W,ipiv);
95: }
96: return(0);
97: }
99: PetscErrorCode FNEvaluateFunctionMat_Rational(FN fn,Mat A,Mat B)
100: {
101: PetscErrorCode ierr;
102: PetscInt m;
103: const PetscScalar *Aa;
104: PetscScalar *Ba;
107: MatDenseGetArrayRead(A,&Aa);
108: MatDenseGetArray(B,&Ba);
109: MatGetSize(A,&m,NULL);
110: FNEvaluateFunctionMat_Rational_Private(fn,Aa,Ba,m,PETSC_FALSE);
111: MatDenseRestoreArrayRead(A,&Aa);
112: MatDenseRestoreArray(B,&Ba);
113: return(0);
114: }
116: PetscErrorCode FNEvaluateFunctionMatVec_Rational(FN fn,Mat A,Vec v)
117: {
118: PetscErrorCode ierr;
119: PetscInt m;
120: const PetscScalar *Aa;
121: PetscScalar *Ba;
122: Mat B;
125: FN_AllocateWorkMat(fn,A,&B);
126: MatDenseGetArrayRead(A,&Aa);
127: MatDenseGetArray(B,&Ba);
128: MatGetSize(A,&m,NULL);
129: FNEvaluateFunctionMat_Rational_Private(fn,Aa,Ba,m,PETSC_TRUE);
130: MatDenseRestoreArrayRead(A,&Aa);
131: MatDenseRestoreArray(B,&Ba);
132: MatGetColumnVector(B,v,0);
133: FN_FreeWorkMat(fn,&B);
134: return(0);
135: }
137: PetscErrorCode FNEvaluateDerivative_Rational(FN fn,PetscScalar x,PetscScalar *yp)
138: {
139: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
140: PetscInt i;
141: PetscScalar p,q,pp,qp;
144: if (!ctx->np) {
145: p = 1.0;
146: pp = 0.0;
147: } else {
148: p = ctx->pcoeff[0];
149: pp = 0.0;
150: for (i=1;i<ctx->np;i++) {
151: pp = p+x*pp;
152: p = ctx->pcoeff[i]+x*p;
153: }
154: }
155: if (!ctx->nq) *yp = pp;
156: else {
157: q = ctx->qcoeff[0];
158: qp = 0.0;
159: for (i=1;i<ctx->nq;i++) {
160: qp = q+x*qp;
161: q = ctx->qcoeff[i]+x*q;
162: }
163: if (q==0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
164: *yp = (pp*q-p*qp)/(q*q);
165: }
166: return(0);
167: }
169: PetscErrorCode FNView_Rational(FN fn,PetscViewer viewer)
170: {
172: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
173: PetscBool isascii;
174: PetscInt i;
175: char str[50];
178: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
179: if (isascii) {
180: if (fn->alpha!=(PetscScalar)1.0 || fn->beta!=(PetscScalar)1.0) {
181: SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_FALSE);
182: PetscViewerASCIIPrintf(viewer," Scale factors: alpha=%s,",str);
183: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
184: SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_FALSE);
185: PetscViewerASCIIPrintf(viewer," beta=%s\n",str);
186: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
187: }
188: if (!ctx->nq) {
189: if (!ctx->np) {
190: PetscViewerASCIIPrintf(viewer," Constant: 1.0\n");
191: } else if (ctx->np==1) {
192: SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[0],PETSC_FALSE);
193: PetscViewerASCIIPrintf(viewer," Constant: %s\n",str);
194: } else {
195: PetscViewerASCIIPrintf(viewer," Polynomial: ");
196: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
197: for (i=0;i<ctx->np-1;i++) {
198: SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[i],PETSC_TRUE);
199: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->np-i-1);
200: }
201: SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[ctx->np-1],PETSC_TRUE);
202: PetscViewerASCIIPrintf(viewer,"%s\n",str);
203: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
204: }
205: } else if (!ctx->np) {
206: PetscViewerASCIIPrintf(viewer," Inverse polinomial: 1 / (");
207: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
208: for (i=0;i<ctx->nq-1;i++) {
209: SlepcSNPrintfScalar(str,sizeof(str),ctx->qcoeff[i],PETSC_TRUE);
210: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->nq-i-1);
211: }
212: SlepcSNPrintfScalar(str,sizeof(str),ctx->qcoeff[ctx->nq-1],PETSC_TRUE);
213: PetscViewerASCIIPrintf(viewer,"%s)\n",str);
214: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
215: } else {
216: PetscViewerASCIIPrintf(viewer," Rational function: (");
217: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
218: for (i=0;i<ctx->np-1;i++) {
219: SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[i],PETSC_TRUE);
220: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->np-i-1);
221: }
222: SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[ctx->np-1],PETSC_TRUE);
223: PetscViewerASCIIPrintf(viewer,"%s) / (",str);
224: for (i=0;i<ctx->nq-1;i++) {
225: SlepcSNPrintfScalar(str,sizeof(str),ctx->qcoeff[i],PETSC_TRUE);
226: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->nq-i-1);
227: }
228: SlepcSNPrintfScalar(str,sizeof(str),ctx->qcoeff[ctx->nq-1],PETSC_TRUE);
229: PetscViewerASCIIPrintf(viewer,"%s)\n",str);
230: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
231: }
232: }
233: return(0);
234: }
236: static PetscErrorCode FNRationalSetNumerator_Rational(FN fn,PetscInt np,PetscScalar *pcoeff)
237: {
239: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
240: PetscInt i;
243: if (np<0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Argument np cannot be negative");
244: ctx->np = np;
245: PetscFree(ctx->pcoeff);
246: if (np) {
247: PetscMalloc1(np,&ctx->pcoeff);
248: PetscLogObjectMemory((PetscObject)fn,np*sizeof(PetscScalar));
249: for (i=0;i<np;i++) ctx->pcoeff[i] = pcoeff[i];
250: }
251: return(0);
252: }
254: /*@C
255: FNRationalSetNumerator - Sets the parameters defining the numerator of the
256: rational function.
258: Logically Collective on fn
260: Input Parameters:
261: + fn - the math function context
262: . np - number of coefficients
263: - pcoeff - coefficients (array of scalar values)
265: Notes:
266: Let the rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials.
267: This function provides the coefficients of the numerator p(x).
268: Hence, p(x) is of degree np-1.
269: If np is zero, then the numerator is assumed to be p(x)=1.
271: In polynomials, high order coefficients are stored in the first positions
272: of the array, e.g. to represent x^2-3 use {1,0,-3}.
274: Level: intermediate
276: .seealso: FNRationalSetDenominator(), FNRationalGetNumerator()
277: @*/
278: PetscErrorCode FNRationalSetNumerator(FN fn,PetscInt np,PetscScalar *pcoeff)
279: {
286: PetscTryMethod(fn,"FNRationalSetNumerator_C",(FN,PetscInt,PetscScalar*),(fn,np,pcoeff));
287: return(0);
288: }
290: static PetscErrorCode FNRationalGetNumerator_Rational(FN fn,PetscInt *np,PetscScalar *pcoeff[])
291: {
293: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
294: PetscInt i;
297: if (np) *np = ctx->np;
298: if (pcoeff) {
299: if (!ctx->np) *pcoeff = NULL;
300: else {
301: PetscMalloc1(ctx->np,pcoeff);
302: for (i=0;i<ctx->np;i++) (*pcoeff)[i] = ctx->pcoeff[i];
303: }
304: }
305: return(0);
306: }
308: /*@C
309: FNRationalGetNumerator - Gets the parameters that define the numerator of the
310: rational function.
312: Not Collective
314: Input Parameter:
315: . fn - the math function context
317: Output Parameters:
318: + np - number of coefficients
319: - pcoeff - coefficients (array of scalar values, length nq)
321: Notes:
322: The values passed by user with FNRationalSetNumerator() are returned (or null
323: pointers otherwise).
324: The pcoeff array should be freed by the user when no longer needed.
326: Level: intermediate
328: .seealso: FNRationalSetNumerator()
329: @*/
330: PetscErrorCode FNRationalGetNumerator(FN fn,PetscInt *np,PetscScalar *pcoeff[])
331: {
336: PetscUseMethod(fn,"FNRationalGetNumerator_C",(FN,PetscInt*,PetscScalar**),(fn,np,pcoeff));
337: return(0);
338: }
340: static PetscErrorCode FNRationalSetDenominator_Rational(FN fn,PetscInt nq,PetscScalar *qcoeff)
341: {
343: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
344: PetscInt i;
347: if (nq<0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Argument nq cannot be negative");
348: ctx->nq = nq;
349: PetscFree(ctx->qcoeff);
350: if (nq) {
351: PetscMalloc1(nq,&ctx->qcoeff);
352: PetscLogObjectMemory((PetscObject)fn,nq*sizeof(PetscScalar));
353: for (i=0;i<nq;i++) ctx->qcoeff[i] = qcoeff[i];
354: }
355: return(0);
356: }
358: /*@C
359: FNRationalSetDenominator - Sets the parameters defining the denominator of the
360: rational function.
362: Logically Collective on fn
364: Input Parameters:
365: + fn - the math function context
366: . nq - number of coefficients
367: - qcoeff - coefficients (array of scalar values)
369: Notes:
370: Let the rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials.
371: This function provides the coefficients of the denominator q(x).
372: Hence, q(x) is of degree nq-1.
373: If nq is zero, then the function is assumed to be polynomial, r(x) = p(x).
375: In polynomials, high order coefficients are stored in the first positions
376: of the array, e.g. to represent x^2-3 use {1,0,-3}.
378: Level: intermediate
380: .seealso: FNRationalSetNumerator(), FNRationalGetDenominator()
381: @*/
382: PetscErrorCode FNRationalSetDenominator(FN fn,PetscInt nq,PetscScalar *qcoeff)
383: {
390: PetscTryMethod(fn,"FNRationalSetDenominator_C",(FN,PetscInt,PetscScalar*),(fn,nq,qcoeff));
391: return(0);
392: }
394: static PetscErrorCode FNRationalGetDenominator_Rational(FN fn,PetscInt *nq,PetscScalar *qcoeff[])
395: {
397: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
398: PetscInt i;
401: if (nq) *nq = ctx->nq;
402: if (qcoeff) {
403: if (!ctx->nq) *qcoeff = NULL;
404: else {
405: PetscMalloc1(ctx->nq,qcoeff);
406: for (i=0;i<ctx->nq;i++) (*qcoeff)[i] = ctx->qcoeff[i];
407: }
408: }
409: return(0);
410: }
412: /*@C
413: FNRationalGetDenominator - Gets the parameters that define the denominator of the
414: rational function.
416: Not Collective
418: Input Parameter:
419: . fn - the math function context
421: Output Parameters:
422: + nq - number of coefficients
423: - qcoeff - coefficients (array of scalar values, length nq)
425: Notes:
426: The values passed by user with FNRationalSetDenominator() are returned (or a null
427: pointer otherwise).
428: The qcoeff array should be freed by the user when no longer needed.
430: Level: intermediate
432: .seealso: FNRationalSetDenominator()
433: @*/
434: PetscErrorCode FNRationalGetDenominator(FN fn,PetscInt *nq,PetscScalar *qcoeff[])
435: {
440: PetscUseMethod(fn,"FNRationalGetDenominator_C",(FN,PetscInt*,PetscScalar**),(fn,nq,qcoeff));
441: return(0);
442: }
444: PetscErrorCode FNSetFromOptions_Rational(PetscOptionItems *PetscOptionsObject,FN fn)
445: {
447: #define PARMAX 10
448: PetscScalar array[PARMAX];
449: PetscInt i,k;
450: PetscBool flg;
453: PetscOptionsHead(PetscOptionsObject,"FN Rational Options");
455: k = PARMAX;
456: for (i=0;i<k;i++) array[i] = 0;
457: PetscOptionsScalarArray("-fn_rational_numerator","Numerator coefficients (one or more scalar values separated with a comma without spaces)","FNRationalSetNumerator",array,&k,&flg);
458: if (flg) { FNRationalSetNumerator(fn,k,array); }
460: k = PARMAX;
461: for (i=0;i<k;i++) array[i] = 0;
462: PetscOptionsScalarArray("-fn_rational_denominator","Denominator coefficients (one or more scalar values separated with a comma without spaces)","FNRationalSetDenominator",array,&k,&flg);
463: if (flg) { FNRationalSetDenominator(fn,k,array); }
465: PetscOptionsTail();
466: return(0);
467: }
469: PetscErrorCode FNDuplicate_Rational(FN fn,MPI_Comm comm,FN *newfn)
470: {
472: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data,*ctx2 = (FN_RATIONAL*)(*newfn)->data;
473: PetscInt i;
476: ctx2->np = ctx->np;
477: if (ctx->np) {
478: PetscMalloc1(ctx->np,&ctx2->pcoeff);
479: PetscLogObjectMemory((PetscObject)(*newfn),ctx->np*sizeof(PetscScalar));
480: for (i=0;i<ctx->np;i++) ctx2->pcoeff[i] = ctx->pcoeff[i];
481: }
482: ctx2->nq = ctx->nq;
483: if (ctx->nq) {
484: PetscMalloc1(ctx->nq,&ctx2->qcoeff);
485: PetscLogObjectMemory((PetscObject)(*newfn),ctx->nq*sizeof(PetscScalar));
486: for (i=0;i<ctx->nq;i++) ctx2->qcoeff[i] = ctx->qcoeff[i];
487: }
488: return(0);
489: }
491: PetscErrorCode FNDestroy_Rational(FN fn)
492: {
494: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
497: PetscFree(ctx->pcoeff);
498: PetscFree(ctx->qcoeff);
499: PetscFree(fn->data);
500: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetNumerator_C",NULL);
501: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetNumerator_C",NULL);
502: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetDenominator_C",NULL);
503: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetDenominator_C",NULL);
504: return(0);
505: }
507: SLEPC_EXTERN PetscErrorCode FNCreate_Rational(FN fn)
508: {
510: FN_RATIONAL *ctx;
513: PetscNewLog(fn,&ctx);
514: fn->data = (void*)ctx;
516: fn->ops->evaluatefunction = FNEvaluateFunction_Rational;
517: fn->ops->evaluatederivative = FNEvaluateDerivative_Rational;
518: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Rational;
519: fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Rational;
520: fn->ops->setfromoptions = FNSetFromOptions_Rational;
521: fn->ops->view = FNView_Rational;
522: fn->ops->duplicate = FNDuplicate_Rational;
523: fn->ops->destroy = FNDestroy_Rational;
524: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetNumerator_C",FNRationalSetNumerator_Rational);
525: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetNumerator_C",FNRationalGetNumerator_Rational);
526: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetDenominator_C",FNRationalSetDenominator_Rational);
527: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetDenominator_C",FNRationalGetDenominator_Rational);
528: return(0);
529: }