Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : /*
      11             :    Phi functions
      12             :       phi_0(x) = exp(x)
      13             :       phi_1(x) = (exp(x)-1)/x
      14             :       phi_k(x) = (phi_{k-1}(x)-1/(k-1)!)/x
      15             : */
      16             : 
      17             : #include <slepc/private/fnimpl.h>      /*I "slepcfn.h" I*/
      18             : 
      19             : typedef struct {
      20             :   PetscInt k;    /* index of the phi-function, defaults to k=1 */
      21             :   Mat      H;    /* auxiliary matrix of order m+k */
      22             :   Mat      F;    /* auxiliary matrix to store exp(H) */
      23             : } FN_PHI;
      24             : 
      25             : #define MAX_INDEX 10
      26             : 
      27             : static const PetscReal rfactorial[MAX_INDEX+2] = { 1, 1, 0.5, 1.0/6, 1.0/24, 1.0/120, 1.0/720, 1.0/5040, 1.0/40320, 1.0/362880, 1.0/3628800, 1.0/39916800 };
      28             : 
      29          27 : static PetscErrorCode FNEvaluateFunction_Phi(FN fn,PetscScalar x,PetscScalar *y)
      30             : {
      31          27 :   FN_PHI      *ctx = (FN_PHI*)fn->data;
      32          27 :   PetscInt    i;
      33          27 :   PetscScalar phi[MAX_INDEX+1];
      34             : 
      35          27 :   PetscFunctionBegin;
      36          27 :   if (x==0.0) *y = rfactorial[ctx->k];
      37             :   else {
      38          27 :     phi[0] = PetscExpScalar(x);
      39          63 :     for (i=1;i<=ctx->k;i++) phi[i] = (phi[i-1]-rfactorial[i-1])/x;
      40          27 :     *y = phi[ctx->k];
      41             :   }
      42          27 :   PetscFunctionReturn(PETSC_SUCCESS);
      43             : }
      44             : 
      45           3 : static PetscErrorCode FNEvaluateDerivative_Phi(FN fn,PetscScalar x,PetscScalar *y)
      46             : {
      47           3 :   FN_PHI      *ctx = (FN_PHI*)fn->data;
      48           3 :   PetscInt    i;
      49           3 :   PetscScalar phi[MAX_INDEX+2];
      50             : 
      51           3 :   PetscFunctionBegin;
      52           3 :   if (x==0.0) *y = rfactorial[ctx->k+1];
      53             :   else {
      54           3 :     phi[0] = PetscExpScalar(x);
      55          10 :     for (i=1;i<=ctx->k+1;i++) phi[i] = (phi[i-1]-rfactorial[i-1])/x;
      56           3 :     *y = phi[ctx->k] - phi[ctx->k+1]*(PetscReal)ctx->k;
      57             :   }
      58           3 :   PetscFunctionReturn(PETSC_SUCCESS);
      59             : }
      60             : 
      61          23 : static PetscErrorCode FNEvaluateFunctionMatVec_Phi(FN fn,Mat A,Vec v)
      62             : {
      63          23 :   FN_PHI            *ctx = (FN_PHI*)fn->data;
      64          23 :   PetscInt          i,j,m,n,nh;
      65          23 :   PetscScalar       *Ha,*va,sfactor=1.0;
      66          23 :   const PetscScalar *Aa,*Fa;
      67             : 
      68          23 :   PetscFunctionBegin;
      69          23 :   PetscCall(MatGetSize(A,&m,NULL));
      70          23 :   n = m+ctx->k;
      71          23 :   if (ctx->H) {
      72          19 :     PetscCall(MatGetSize(ctx->H,&nh,NULL));
      73          19 :     if (n!=nh) {
      74          19 :       PetscCall(MatDestroy(&ctx->H));
      75          19 :       PetscCall(MatDestroy(&ctx->F));
      76             :     }
      77             :   }
      78          23 :   if (!ctx->H) {
      79          23 :     PetscCall(MatCreateDense(PETSC_COMM_SELF,n,n,n,n,NULL,&ctx->H));
      80          23 :     PetscCall(MatCreateDense(PETSC_COMM_SELF,n,n,n,n,NULL,&ctx->F));
      81             :   }
      82          23 :   PetscCall(MatDenseGetArray(ctx->H,&Ha));
      83          23 :   PetscCall(MatDenseGetArrayRead(A,&Aa));
      84        1547 :   for (j=0;j<m;j++) PetscCall(PetscArraycpy(Ha+j*n,Aa+j*m,m));
      85          23 :   PetscCall(MatDenseRestoreArrayRead(A,&Aa));
      86          23 :   if (ctx->k) {
      87        3070 :     for (j=0;j<m;j++) for (i=m;i<n;i++) Ha[i+j*n] = 0.0;
      88        1608 :     for (j=m;j<n;j++) for (i=0;i<n;i++) Ha[i+j*n] = 0.0;
      89          22 :     Ha[0+m*n] = fn->alpha;
      90          24 :     for (j=m+1;j<n;j++) Ha[j-1+j*n] = fn->alpha;
      91             :   }
      92          23 :   PetscCall(MatDenseRestoreArray(ctx->H,&Ha));
      93             : 
      94          23 :   PetscCall(FNEvaluateFunctionMat_Exp_Higham(fn,ctx->H,ctx->F));
      95             : 
      96          23 :   PetscCall(MatDenseGetArrayRead(ctx->F,&Fa));
      97          23 :   PetscCall(VecGetArray(v,&va));
      98          23 :   if (ctx->k) {
      99          22 :     sfactor = PetscPowScalarInt(fn->alpha,-ctx->k);
     100        1560 :     for (i=0;i<m;i++) va[i] = sfactor*Fa[i+(n-1)*n];
     101             :   } else {
     102           9 :     for (i=0;i<m;i++) va[i] = sfactor*Fa[i+0*n];
     103             :   }
     104          23 :   PetscCall(VecRestoreArray(v,&va));
     105          23 :   PetscCall(MatDenseRestoreArrayRead(ctx->F,&Fa));
     106          23 :   PetscFunctionReturn(PETSC_SUCCESS);
     107             : }
     108             : 
     109           3 : static PetscErrorCode FNPhiSetIndex_Phi(FN fn,PetscInt k)
     110             : {
     111           3 :   FN_PHI         *ctx = (FN_PHI*)fn->data;
     112             : 
     113           3 :   PetscFunctionBegin;
     114           3 :   PetscCheck(k>=0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Index cannot be negative");
     115           3 :   PetscCheck(k<=MAX_INDEX,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Phi functions only implemented for k<=%d",MAX_INDEX);
     116           3 :   if (k!=ctx->k) {
     117           2 :     ctx->k = k;
     118           2 :     PetscCall(MatDestroy(&ctx->H));
     119           2 :     PetscCall(MatDestroy(&ctx->F));
     120             :   }
     121           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     122             : }
     123             : 
     124             : /*@
     125             :    FNPhiSetIndex - Sets the index of the phi-function.
     126             : 
     127             :    Logically Collective
     128             : 
     129             :    Input Parameters:
     130             : +  fn - the math function context
     131             : -  k  - the index
     132             : 
     133             :    Notes:
     134             :    The phi-functions are defined as follows. The default is k=1.
     135             : .vb
     136             :       phi_0(x) = exp(x)
     137             :       phi_1(x) = (exp(x)-1)/x
     138             :       phi_k(x) = (phi_{k-1}(x)-1/(k-1)!)/x
     139             : .ve
     140             : 
     141             :    Level: intermediate
     142             : 
     143             : .seealso: FNPhiGetIndex()
     144             : @*/
     145           3 : PetscErrorCode FNPhiSetIndex(FN fn,PetscInt k)
     146             : {
     147           3 :   PetscFunctionBegin;
     148           3 :   PetscValidHeaderSpecific(fn,FN_CLASSID,1);
     149           9 :   PetscValidLogicalCollectiveInt(fn,k,2);
     150           3 :   PetscTryMethod(fn,"FNPhiSetIndex_C",(FN,PetscInt),(fn,k));
     151           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     152             : }
     153             : 
     154           1 : static PetscErrorCode FNPhiGetIndex_Phi(FN fn,PetscInt *k)
     155             : {
     156           1 :   FN_PHI *ctx = (FN_PHI*)fn->data;
     157             : 
     158           1 :   PetscFunctionBegin;
     159           1 :   *k = ctx->k;
     160           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     161             : }
     162             : 
     163             : /*@
     164             :    FNPhiGetIndex - Gets the index of the phi-function.
     165             : 
     166             :    Not Collective
     167             : 
     168             :    Input Parameter:
     169             : .  fn - the math function context
     170             : 
     171             :    Output Parameter:
     172             : .  k  - the index
     173             : 
     174             :    Level: intermediate
     175             : 
     176             : .seealso: FNPhiSetIndex()
     177             : @*/
     178           1 : PetscErrorCode FNPhiGetIndex(FN fn,PetscInt *k)
     179             : {
     180           1 :   PetscFunctionBegin;
     181           1 :   PetscValidHeaderSpecific(fn,FN_CLASSID,1);
     182           1 :   PetscAssertPointer(k,2);
     183           1 :   PetscUseMethod(fn,"FNPhiGetIndex_C",(FN,PetscInt*),(fn,k));
     184           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     185             : }
     186             : 
     187           3 : static PetscErrorCode FNView_Phi(FN fn,PetscViewer viewer)
     188             : {
     189           3 :   FN_PHI         *ctx = (FN_PHI*)fn->data;
     190           3 :   PetscBool      isascii;
     191           3 :   char           str[50],strx[50];
     192             : 
     193           3 :   PetscFunctionBegin;
     194           3 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     195           3 :   if (isascii) {
     196           3 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  phi_%" PetscInt_FMT ": ",ctx->k));
     197           3 :     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
     198           3 :     if (fn->beta!=(PetscScalar)1.0) {
     199           1 :       PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_TRUE));
     200           1 :       PetscCall(PetscViewerASCIIPrintf(viewer,"%s*",str));
     201             :     }
     202           3 :     if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscSNPrintf(strx,sizeof(strx),"x"));
     203             :     else {
     204           1 :       PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
     205           1 :       PetscCall(PetscSNPrintf(strx,sizeof(strx),"(%s*x)",str));
     206             :     }
     207           3 :     if (!ctx->k) PetscCall(PetscViewerASCIIPrintf(viewer,"exp(%s)\n",strx));
     208           2 :     else if (ctx->k==1) PetscCall(PetscViewerASCIIPrintf(viewer,"(exp(%s)-1)/%s\n",strx,strx));
     209           1 :     else PetscCall(PetscViewerASCIIPrintf(viewer,"(phi_%" PetscInt_FMT "(%s)-1/%" PetscInt_FMT "!)/%s\n",ctx->k-1,strx,ctx->k-1,strx));
     210           3 :     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
     211             :   }
     212           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     213             : }
     214             : 
     215           2 : static PetscErrorCode FNSetFromOptions_Phi(FN fn,PetscOptionItems *PetscOptionsObject)
     216             : {
     217           2 :   FN_PHI         *ctx = (FN_PHI*)fn->data;
     218           2 :   PetscInt       k;
     219           2 :   PetscBool      flag;
     220             : 
     221           2 :   PetscFunctionBegin;
     222           2 :   PetscOptionsHeadBegin(PetscOptionsObject,"FN Phi Options");
     223             : 
     224           2 :     PetscCall(PetscOptionsInt("-fn_phi_index","Index of the phi-function","FNPhiSetIndex",ctx->k,&k,&flag));
     225           2 :     if (flag) PetscCall(FNPhiSetIndex(fn,k));
     226             : 
     227           2 :   PetscOptionsHeadEnd();
     228           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     229             : }
     230             : 
     231           1 : static PetscErrorCode FNDuplicate_Phi(FN fn,MPI_Comm comm,FN *newfn)
     232             : {
     233           1 :   FN_PHI *ctx = (FN_PHI*)fn->data,*ctx2 = (FN_PHI*)(*newfn)->data;
     234             : 
     235           1 :   PetscFunctionBegin;
     236           1 :   ctx2->k = ctx->k;
     237           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     238             : }
     239             : 
     240           5 : static PetscErrorCode FNDestroy_Phi(FN fn)
     241             : {
     242           5 :   FN_PHI         *ctx = (FN_PHI*)fn->data;
     243             : 
     244           5 :   PetscFunctionBegin;
     245           5 :   PetscCall(MatDestroy(&ctx->H));
     246           5 :   PetscCall(MatDestroy(&ctx->F));
     247           5 :   PetscCall(PetscFree(fn->data));
     248           5 :   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNPhiSetIndex_C",NULL));
     249           5 :   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNPhiGetIndex_C",NULL));
     250           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     251             : }
     252             : 
     253           5 : SLEPC_EXTERN PetscErrorCode FNCreate_Phi(FN fn)
     254             : {
     255           5 :   FN_PHI         *ctx;
     256             : 
     257           5 :   PetscFunctionBegin;
     258           5 :   PetscCall(PetscNew(&ctx));
     259           5 :   fn->data = (void*)ctx;
     260           5 :   ctx->k   = 1;
     261             : 
     262           5 :   fn->ops->evaluatefunction          = FNEvaluateFunction_Phi;
     263           5 :   fn->ops->evaluatederivative        = FNEvaluateDerivative_Phi;
     264           5 :   fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Phi;
     265           5 :   fn->ops->setfromoptions            = FNSetFromOptions_Phi;
     266           5 :   fn->ops->view                      = FNView_Phi;
     267           5 :   fn->ops->duplicate                 = FNDuplicate_Phi;
     268           5 :   fn->ops->destroy                   = FNDestroy_Phi;
     269           5 :   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNPhiSetIndex_C",FNPhiSetIndex_Phi));
     270           5 :   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNPhiGetIndex_C",FNPhiGetIndex_Phi));
     271           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     272             : }