Actual source code: fnrational.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }