Actual source code: fnrational.c
slepc-3.13.2 2020-05-12
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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,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: {
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(Aa,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: {
102: PetscInt m;
103: PetscScalar *Aa,*Ba;
106: MatDenseGetArray(A,&Aa);
107: MatDenseGetArray(B,&Ba);
108: MatGetSize(A,&m,NULL);
109: FNEvaluateFunctionMat_Rational_Private(fn,Aa,Ba,m,PETSC_FALSE);
110: MatDenseRestoreArray(A,&Aa);
111: MatDenseRestoreArray(B,&Ba);
112: return(0);
113: }
115: PetscErrorCode FNEvaluateFunctionMatVec_Rational(FN fn,Mat A,Vec v)
116: {
118: PetscInt m;
119: PetscScalar *Aa,*Ba;
120: Mat B;
123: FN_AllocateWorkMat(fn,A,&B);
124: MatDenseGetArray(A,&Aa);
125: MatDenseGetArray(B,&Ba);
126: MatGetSize(A,&m,NULL);
127: FNEvaluateFunctionMat_Rational_Private(fn,Aa,Ba,m,PETSC_TRUE);
128: MatDenseRestoreArray(A,&Aa);
129: MatDenseRestoreArray(B,&Ba);
130: MatGetColumnVector(B,v,0);
131: FN_FreeWorkMat(fn,&B);
132: return(0);
133: }
135: PetscErrorCode FNEvaluateDerivative_Rational(FN fn,PetscScalar x,PetscScalar *yp)
136: {
137: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
138: PetscInt i;
139: PetscScalar p,q,pp,qp;
142: if (!ctx->np) {
143: p = 1.0;
144: pp = 0.0;
145: } else {
146: p = ctx->pcoeff[0];
147: pp = 0.0;
148: for (i=1;i<ctx->np;i++) {
149: pp = p+x*pp;
150: p = ctx->pcoeff[i]+x*p;
151: }
152: }
153: if (!ctx->nq) *yp = pp;
154: else {
155: q = ctx->qcoeff[0];
156: qp = 0.0;
157: for (i=1;i<ctx->nq;i++) {
158: qp = q+x*qp;
159: q = ctx->qcoeff[i]+x*q;
160: }
161: if (q==0.0) SETERRQ(PETSC_COMM_SELF,1,"Derivative not defined in the requested value");
162: *yp = (pp*q-p*qp)/(q*q);
163: }
164: return(0);
165: }
167: PetscErrorCode FNView_Rational(FN fn,PetscViewer viewer)
168: {
170: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
171: PetscBool isascii;
172: PetscInt i;
173: char str[50];
176: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
177: if (isascii) {
178: if (fn->alpha!=(PetscScalar)1.0 || fn->beta!=(PetscScalar)1.0) {
179: SlepcSNPrintfScalar(str,50,fn->alpha,PETSC_FALSE);
180: PetscViewerASCIIPrintf(viewer," Scale factors: alpha=%s,",str);
181: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
182: SlepcSNPrintfScalar(str,50,fn->beta,PETSC_FALSE);
183: PetscViewerASCIIPrintf(viewer," beta=%s\n",str);
184: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
185: }
186: if (!ctx->nq) {
187: if (!ctx->np) {
188: PetscViewerASCIIPrintf(viewer," Constant: 1.0\n");
189: } else if (ctx->np==1) {
190: SlepcSNPrintfScalar(str,50,ctx->pcoeff[0],PETSC_FALSE);
191: PetscViewerASCIIPrintf(viewer," Constant: %s\n",str);
192: } else {
193: PetscViewerASCIIPrintf(viewer," Polynomial: ");
194: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
195: for (i=0;i<ctx->np-1;i++) {
196: SlepcSNPrintfScalar(str,50,ctx->pcoeff[i],PETSC_TRUE);
197: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->np-i-1);
198: }
199: SlepcSNPrintfScalar(str,50,ctx->pcoeff[ctx->np-1],PETSC_TRUE);
200: PetscViewerASCIIPrintf(viewer,"%s\n",str);
201: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
202: }
203: } else if (!ctx->np) {
204: PetscViewerASCIIPrintf(viewer," Inverse polinomial: 1 / (");
205: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
206: for (i=0;i<ctx->nq-1;i++) {
207: SlepcSNPrintfScalar(str,50,ctx->qcoeff[i],PETSC_TRUE);
208: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->nq-i-1);
209: }
210: SlepcSNPrintfScalar(str,50,ctx->qcoeff[ctx->nq-1],PETSC_TRUE);
211: PetscViewerASCIIPrintf(viewer,"%s)\n",str);
212: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
213: } else {
214: PetscViewerASCIIPrintf(viewer," Rational function: (");
215: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
216: for (i=0;i<ctx->np-1;i++) {
217: SlepcSNPrintfScalar(str,50,ctx->pcoeff[i],PETSC_TRUE);
218: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->np-i-1);
219: }
220: SlepcSNPrintfScalar(str,50,ctx->pcoeff[ctx->np-1],PETSC_TRUE);
221: PetscViewerASCIIPrintf(viewer,"%s) / (",str);
222: for (i=0;i<ctx->nq-1;i++) {
223: SlepcSNPrintfScalar(str,50,ctx->qcoeff[i],PETSC_TRUE);
224: PetscViewerASCIIPrintf(viewer,"%s*x^%1D",str,ctx->nq-i-1);
225: }
226: SlepcSNPrintfScalar(str,50,ctx->qcoeff[ctx->nq-1],PETSC_TRUE);
227: PetscViewerASCIIPrintf(viewer,"%s)\n",str);
228: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
229: }
230: }
231: return(0);
232: }
234: static PetscErrorCode FNRationalSetNumerator_Rational(FN fn,PetscInt np,PetscScalar *pcoeff)
235: {
237: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
238: PetscInt i;
241: if (np<0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Argument np cannot be negative");
242: ctx->np = np;
243: PetscFree(ctx->pcoeff);
244: if (np) {
245: PetscMalloc1(np,&ctx->pcoeff);
246: PetscLogObjectMemory((PetscObject)fn,np*sizeof(PetscScalar));
247: for (i=0;i<np;i++) ctx->pcoeff[i] = pcoeff[i];
248: }
249: return(0);
250: }
252: /*@C
253: FNRationalSetNumerator - Sets the parameters defining the numerator of the
254: rational function.
256: Logically Collective on fn
258: Input Parameters:
259: + fn - the math function context
260: . np - number of coefficients
261: - pcoeff - coefficients (array of scalar values)
263: Notes:
264: Let the rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials.
265: This function provides the coefficients of the numerator p(x).
266: Hence, p(x) is of degree np-1.
267: If np is zero, then the numerator is assumed to be p(x)=1.
269: In polynomials, high order coefficients are stored in the first positions
270: of the array, e.g. to represent x^2-3 use {1,0,-3}.
272: Level: intermediate
274: .seealso: FNRationalSetDenominator(), FNRationalGetNumerator()
275: @*/
276: PetscErrorCode FNRationalSetNumerator(FN fn,PetscInt np,PetscScalar *pcoeff)
277: {
284: PetscTryMethod(fn,"FNRationalSetNumerator_C",(FN,PetscInt,PetscScalar*),(fn,np,pcoeff));
285: return(0);
286: }
288: static PetscErrorCode FNRationalGetNumerator_Rational(FN fn,PetscInt *np,PetscScalar *pcoeff[])
289: {
291: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
292: PetscInt i;
295: if (np) *np = ctx->np;
296: if (pcoeff) {
297: if (!ctx->np) *pcoeff = NULL;
298: else {
299: PetscMalloc1(ctx->np,pcoeff);
300: for (i=0;i<ctx->np;i++) (*pcoeff)[i] = ctx->pcoeff[i];
301: }
302: }
303: return(0);
304: }
306: /*@C
307: FNRationalGetNumerator - Gets the parameters that define the numerator of the
308: rational function.
310: Not Collective
312: Input Parameter:
313: . fn - the math function context
315: Output Parameters:
316: + np - number of coefficients
317: - pcoeff - coefficients (array of scalar values, length nq)
319: Notes:
320: The values passed by user with FNRationalSetNumerator() are returned (or null
321: pointers otherwise).
322: The pcoeff array should be freed by the user when no longer needed.
324: Level: intermediate
326: .seealso: FNRationalSetNumerator()
327: @*/
328: PetscErrorCode FNRationalGetNumerator(FN fn,PetscInt *np,PetscScalar *pcoeff[])
329: {
334: PetscUseMethod(fn,"FNRationalGetNumerator_C",(FN,PetscInt*,PetscScalar**),(fn,np,pcoeff));
335: return(0);
336: }
338: static PetscErrorCode FNRationalSetDenominator_Rational(FN fn,PetscInt nq,PetscScalar *qcoeff)
339: {
341: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
342: PetscInt i;
345: if (nq<0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Argument nq cannot be negative");
346: ctx->nq = nq;
347: PetscFree(ctx->qcoeff);
348: if (nq) {
349: PetscMalloc1(nq,&ctx->qcoeff);
350: PetscLogObjectMemory((PetscObject)fn,nq*sizeof(PetscScalar));
351: for (i=0;i<nq;i++) ctx->qcoeff[i] = qcoeff[i];
352: }
353: return(0);
354: }
356: /*@C
357: FNRationalSetDenominator - Sets the parameters defining the denominator of the
358: rational function.
360: Logically Collective on fn
362: Input Parameters:
363: + fn - the math function context
364: . nq - number of coefficients
365: - qcoeff - coefficients (array of scalar values)
367: Notes:
368: Let the rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials.
369: This function provides the coefficients of the denominator q(x).
370: Hence, q(x) is of degree nq-1.
371: If nq is zero, then the function is assumed to be polynomial, r(x) = p(x).
373: In polynomials, high order coefficients are stored in the first positions
374: of the array, e.g. to represent x^2-3 use {1,0,-3}.
376: Level: intermediate
378: .seealso: FNRationalSetNumerator(), FNRationalGetDenominator()
379: @*/
380: PetscErrorCode FNRationalSetDenominator(FN fn,PetscInt nq,PetscScalar *qcoeff)
381: {
388: PetscTryMethod(fn,"FNRationalSetDenominator_C",(FN,PetscInt,PetscScalar*),(fn,nq,qcoeff));
389: return(0);
390: }
392: static PetscErrorCode FNRationalGetDenominator_Rational(FN fn,PetscInt *nq,PetscScalar *qcoeff[])
393: {
395: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
396: PetscInt i;
399: if (nq) *nq = ctx->nq;
400: if (qcoeff) {
401: if (!ctx->nq) *qcoeff = NULL;
402: else {
403: PetscMalloc1(ctx->nq,qcoeff);
404: for (i=0;i<ctx->nq;i++) (*qcoeff)[i] = ctx->qcoeff[i];
405: }
406: }
407: return(0);
408: }
410: /*@C
411: FNRationalGetDenominator - Gets the parameters that define the denominator of the
412: rational function.
414: Not Collective
416: Input Parameter:
417: . fn - the math function context
419: Output Parameters:
420: + nq - number of coefficients
421: - qcoeff - coefficients (array of scalar values, length nq)
423: Notes:
424: The values passed by user with FNRationalSetDenominator() are returned (or a null
425: pointer otherwise).
426: The qcoeff array should be freed by the user when no longer needed.
428: Level: intermediate
430: .seealso: FNRationalSetDenominator()
431: @*/
432: PetscErrorCode FNRationalGetDenominator(FN fn,PetscInt *nq,PetscScalar *qcoeff[])
433: {
438: PetscUseMethod(fn,"FNRationalGetDenominator_C",(FN,PetscInt*,PetscScalar**),(fn,nq,qcoeff));
439: return(0);
440: }
442: PetscErrorCode FNSetFromOptions_Rational(PetscOptionItems *PetscOptionsObject,FN fn)
443: {
445: #define PARMAX 10
446: PetscScalar array[PARMAX];
447: PetscInt i,k;
448: PetscBool flg;
451: PetscOptionsHead(PetscOptionsObject,"FN Rational Options");
453: k = PARMAX;
454: for (i=0;i<k;i++) array[i] = 0;
455: PetscOptionsScalarArray("-fn_rational_numerator","Numerator coefficients (one or more scalar values separated with a comma without spaces)","FNRationalSetNumerator",array,&k,&flg);
456: if (flg) { FNRationalSetNumerator(fn,k,array); }
458: k = PARMAX;
459: for (i=0;i<k;i++) array[i] = 0;
460: PetscOptionsScalarArray("-fn_rational_denominator","Denominator coefficients (one or more scalar values separated with a comma without spaces)","FNRationalSetDenominator",array,&k,&flg);
461: if (flg) { FNRationalSetDenominator(fn,k,array); }
463: PetscOptionsTail();
464: return(0);
465: }
467: PetscErrorCode FNDuplicate_Rational(FN fn,MPI_Comm comm,FN *newfn)
468: {
470: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data,*ctx2 = (FN_RATIONAL*)(*newfn)->data;
471: PetscInt i;
474: ctx2->np = ctx->np;
475: if (ctx->np) {
476: PetscMalloc1(ctx->np,&ctx2->pcoeff);
477: PetscLogObjectMemory((PetscObject)(*newfn),ctx->np*sizeof(PetscScalar));
478: for (i=0;i<ctx->np;i++) ctx2->pcoeff[i] = ctx->pcoeff[i];
479: }
480: ctx2->nq = ctx->nq;
481: if (ctx->nq) {
482: PetscMalloc1(ctx->nq,&ctx2->qcoeff);
483: PetscLogObjectMemory((PetscObject)(*newfn),ctx->nq*sizeof(PetscScalar));
484: for (i=0;i<ctx->nq;i++) ctx2->qcoeff[i] = ctx->qcoeff[i];
485: }
486: return(0);
487: }
489: PetscErrorCode FNDestroy_Rational(FN fn)
490: {
492: FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
495: PetscFree(ctx->pcoeff);
496: PetscFree(ctx->qcoeff);
497: PetscFree(fn->data);
498: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetNumerator_C",NULL);
499: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetNumerator_C",NULL);
500: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetDenominator_C",NULL);
501: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetDenominator_C",NULL);
502: return(0);
503: }
505: SLEPC_EXTERN PetscErrorCode FNCreate_Rational(FN fn)
506: {
508: FN_RATIONAL *ctx;
511: PetscNewLog(fn,&ctx);
512: fn->data = (void*)ctx;
514: fn->ops->evaluatefunction = FNEvaluateFunction_Rational;
515: fn->ops->evaluatederivative = FNEvaluateDerivative_Rational;
516: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Rational;
517: fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Rational;
518: fn->ops->setfromoptions = FNSetFromOptions_Rational;
519: fn->ops->view = FNView_Rational;
520: fn->ops->duplicate = FNDuplicate_Rational;
521: fn->ops->destroy = FNDestroy_Rational;
522: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetNumerator_C",FNRationalSetNumerator_Rational);
523: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetNumerator_C",FNRationalGetNumerator_Rational);
524: PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetDenominator_C",FNRationalSetDenominator_Rational);
525: PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetDenominator_C",FNRationalGetDenominator_Rational);
526: return(0);
527: }