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             :    A function that is obtained by combining two other functions (either by
      12             :    addition, multiplication, division or composition)
      13             : 
      14             :       addition:          f(x) = f1(x)+f2(x)
      15             :       multiplication:    f(x) = f1(x)*f2(x)
      16             :       division:          f(x) = f1(x)/f2(x)      f(A) = f2(A)\f1(A)
      17             :       composition:       f(x) = f2(f1(x))
      18             : */
      19             : 
      20             : #include <slepc/private/fnimpl.h>      /*I "slepcfn.h" I*/
      21             : 
      22             : typedef struct {
      23             :   FN            f1,f2;    /* functions */
      24             :   FNCombineType comb;     /* how the functions are combined */
      25             : } FN_COMBINE;
      26             : 
      27         180 : static PetscErrorCode FNEvaluateFunction_Combine(FN fn,PetscScalar x,PetscScalar *y)
      28             : {
      29         180 :   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
      30         180 :   PetscScalar    a,b;
      31             : 
      32         180 :   PetscFunctionBegin;
      33         180 :   PetscCall(FNEvaluateFunction(ctx->f1,x,&a));
      34         180 :   switch (ctx->comb) {
      35          52 :     case FN_COMBINE_ADD:
      36          52 :       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
      37          52 :       *y = a+b;
      38          52 :       break;
      39          42 :     case FN_COMBINE_MULTIPLY:
      40          42 :       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
      41          42 :       *y = a*b;
      42          42 :       break;
      43          42 :     case FN_COMBINE_DIVIDE:
      44          42 :       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
      45          42 :       PetscCheck(b!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Function not defined in the requested value");
      46          42 :       *y = a/b;
      47          42 :       break;
      48          44 :     case FN_COMBINE_COMPOSE:
      49          44 :       PetscCall(FNEvaluateFunction(ctx->f2,a,y));
      50             :       break;
      51             :   }
      52         180 :   PetscFunctionReturn(PETSC_SUCCESS);
      53             : }
      54             : 
      55           8 : static PetscErrorCode FNEvaluateDerivative_Combine(FN fn,PetscScalar x,PetscScalar *yp)
      56             : {
      57           8 :   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
      58           8 :   PetscScalar    a,b,ap,bp;
      59             : 
      60           8 :   PetscFunctionBegin;
      61           8 :   switch (ctx->comb) {
      62           2 :     case FN_COMBINE_ADD:
      63           2 :       PetscCall(FNEvaluateDerivative(ctx->f1,x,&ap));
      64           2 :       PetscCall(FNEvaluateDerivative(ctx->f2,x,&bp));
      65           2 :       *yp = ap+bp;
      66           2 :       break;
      67           2 :     case FN_COMBINE_MULTIPLY:
      68           2 :       PetscCall(FNEvaluateDerivative(ctx->f1,x,&ap));
      69           2 :       PetscCall(FNEvaluateDerivative(ctx->f2,x,&bp));
      70           2 :       PetscCall(FNEvaluateFunction(ctx->f1,x,&a));
      71           2 :       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
      72           2 :       *yp = ap*b+a*bp;
      73           2 :       break;
      74           2 :     case FN_COMBINE_DIVIDE:
      75           2 :       PetscCall(FNEvaluateDerivative(ctx->f1,x,&ap));
      76           2 :       PetscCall(FNEvaluateDerivative(ctx->f2,x,&bp));
      77           2 :       PetscCall(FNEvaluateFunction(ctx->f1,x,&a));
      78           2 :       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
      79           2 :       PetscCheck(b!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
      80           2 :       *yp = (ap*b-a*bp)/(b*b);
      81           2 :       break;
      82           2 :     case FN_COMBINE_COMPOSE:
      83           2 :       PetscCall(FNEvaluateFunction(ctx->f1,x,&a));
      84           2 :       PetscCall(FNEvaluateDerivative(ctx->f1,x,&ap));
      85           2 :       PetscCall(FNEvaluateDerivative(ctx->f2,a,yp));
      86           2 :       *yp *= ap;
      87           2 :       break;
      88             :   }
      89           8 :   PetscFunctionReturn(PETSC_SUCCESS);
      90             : }
      91             : 
      92           9 : static PetscErrorCode FNEvaluateFunctionMat_Combine(FN fn,Mat A,Mat B)
      93             : {
      94           9 :   FN_COMBINE   *ctx = (FN_COMBINE*)fn->data;
      95           9 :   Mat          W,Z,F;
      96           9 :   PetscBool    iscuda;
      97             : 
      98           9 :   PetscFunctionBegin;
      99           9 :   PetscCall(FN_AllocateWorkMat(fn,A,&W));
     100           9 :   switch (ctx->comb) {
     101           3 :     case FN_COMBINE_ADD:
     102           3 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE));
     103           3 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,A,B,PETSC_FALSE));
     104           3 :       PetscCall(MatAXPY(B,1.0,W,SAME_NONZERO_PATTERN));
     105             :       break;
     106           2 :     case FN_COMBINE_MULTIPLY:
     107           2 :       PetscCall(FN_AllocateWorkMat(fn,A,&Z));
     108           2 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE));
     109           2 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,A,Z,PETSC_FALSE));
     110           2 :       PetscCall(MatMatMult(W,Z,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B));
     111           2 :       PetscCall(FN_FreeWorkMat(fn,&Z));
     112             :       break;
     113           2 :     case FN_COMBINE_DIVIDE:
     114           2 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,A,W,PETSC_FALSE));
     115           2 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,B,PETSC_FALSE));
     116           2 :       PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
     117           4 :       PetscCall(MatGetFactor(W,iscuda?MATSOLVERCUDA:MATSOLVERPETSC,MAT_FACTOR_LU,&F));
     118           2 :       PetscCall(MatLUFactorSymbolic(F,W,NULL,NULL,NULL));
     119           2 :       PetscCall(MatLUFactorNumeric(F,W,NULL));
     120           2 :       PetscCall(MatMatSolve(F,B,B));
     121           2 :       PetscCall(MatDestroy(&F));
     122             :       break;
     123           2 :     case FN_COMBINE_COMPOSE:
     124           2 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE));
     125           2 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,W,B,PETSC_FALSE));
     126             :       break;
     127             :   }
     128           9 :   PetscCall(FN_FreeWorkMat(fn,&W));
     129           9 :   PetscFunctionReturn(PETSC_SUCCESS);
     130             : }
     131             : 
     132          48 : static PetscErrorCode FNEvaluateFunctionMatVec_Combine(FN fn,Mat A,Vec v)
     133             : {
     134          48 :   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
     135          48 :   PetscBool      iscuda;
     136          48 :   Mat            Z,F;
     137          48 :   Vec            w;
     138             : 
     139          48 :   PetscFunctionBegin;
     140          48 :   switch (ctx->comb) {
     141          22 :     case FN_COMBINE_ADD:
     142          22 :       PetscCall(VecDuplicate(v,&w));
     143          22 :       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f1,A,w,PETSC_FALSE));
     144          22 :       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f2,A,v,PETSC_FALSE));
     145          22 :       PetscCall(VecAXPY(v,1.0,w));
     146          22 :       PetscCall(VecDestroy(&w));
     147             :       break;
     148           2 :     case FN_COMBINE_MULTIPLY:
     149           2 :       PetscCall(VecDuplicate(v,&w));
     150           2 :       PetscCall(FN_AllocateWorkMat(fn,A,&Z));
     151           2 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,Z,PETSC_FALSE));
     152           2 :       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f2,A,w,PETSC_FALSE));
     153           2 :       PetscCall(MatMult(Z,w,v));
     154           2 :       PetscCall(FN_FreeWorkMat(fn,&Z));
     155           2 :       PetscCall(VecDestroy(&w));
     156             :       break;
     157          22 :     case FN_COMBINE_DIVIDE:
     158          22 :       PetscCall(VecDuplicate(v,&w));
     159          22 :       PetscCall(FN_AllocateWorkMat(fn,A,&Z));
     160          22 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,A,Z,PETSC_FALSE));
     161          22 :       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f1,A,w,PETSC_FALSE));
     162          22 :       PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
     163          44 :       PetscCall(MatGetFactor(Z,iscuda?MATSOLVERCUDA:MATSOLVERPETSC,MAT_FACTOR_LU,&F));
     164          22 :       PetscCall(MatLUFactorSymbolic(F,Z,NULL,NULL,NULL));
     165          22 :       PetscCall(MatLUFactorNumeric(F,Z,NULL));
     166          22 :       PetscCall(MatSolve(F,w,v));
     167          22 :       PetscCall(MatDestroy(&F));
     168          22 :       PetscCall(FN_FreeWorkMat(fn,&Z));
     169          22 :       PetscCall(VecDestroy(&w));
     170             :       break;
     171           2 :     case FN_COMBINE_COMPOSE:
     172           2 :       PetscCall(FN_AllocateWorkMat(fn,A,&Z));
     173           2 :       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,Z,PETSC_FALSE));
     174           2 :       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f2,Z,v,PETSC_FALSE));
     175           2 :       PetscCall(FN_FreeWorkMat(fn,&Z));
     176             :       break;
     177             :   }
     178          48 :   PetscFunctionReturn(PETSC_SUCCESS);
     179             : }
     180             : 
     181           6 : static PetscErrorCode FNView_Combine(FN fn,PetscViewer viewer)
     182             : {
     183           6 :   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
     184           6 :   PetscBool      isascii;
     185             : 
     186           6 :   PetscFunctionBegin;
     187           6 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     188           6 :   if (isascii) {
     189           6 :     switch (ctx->comb) {
     190           2 :       case FN_COMBINE_ADD:
     191           2 :         PetscCall(PetscViewerASCIIPrintf(viewer,"  two added functions f1+f2\n"));
     192             :         break;
     193           2 :       case FN_COMBINE_MULTIPLY:
     194           2 :         PetscCall(PetscViewerASCIIPrintf(viewer,"  two multiplied functions f1*f2\n"));
     195             :         break;
     196           0 :       case FN_COMBINE_DIVIDE:
     197           0 :         PetscCall(PetscViewerASCIIPrintf(viewer,"  a quotient of two functions f1/f2\n"));
     198             :         break;
     199           2 :       case FN_COMBINE_COMPOSE:
     200           2 :         PetscCall(PetscViewerASCIIPrintf(viewer,"  two composed functions f2(f1(.))\n"));
     201             :         break;
     202             :     }
     203           6 :     PetscCall(PetscViewerASCIIPushTab(viewer));
     204           6 :     PetscCall(FNView(ctx->f1,viewer));
     205           6 :     PetscCall(FNView(ctx->f2,viewer));
     206           6 :     PetscCall(PetscViewerASCIIPopTab(viewer));
     207             :   }
     208           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     209             : }
     210             : 
     211          11 : static PetscErrorCode FNCombineSetChildren_Combine(FN fn,FNCombineType comb,FN f1,FN f2)
     212             : {
     213          11 :   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
     214             : 
     215          11 :   PetscFunctionBegin;
     216          11 :   ctx->comb = comb;
     217          11 :   PetscCall(PetscObjectReference((PetscObject)f1));
     218          11 :   PetscCall(FNDestroy(&ctx->f1));
     219          11 :   ctx->f1 = f1;
     220          11 :   PetscCall(PetscObjectReference((PetscObject)f2));
     221          11 :   PetscCall(FNDestroy(&ctx->f2));
     222          11 :   ctx->f2 = f2;
     223          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     224             : }
     225             : 
     226             : /*@
     227             :    FNCombineSetChildren - Sets the two child functions that constitute this
     228             :    combined function, and the way they must be combined.
     229             : 
     230             :    Logically Collective
     231             : 
     232             :    Input Parameters:
     233             : +  fn   - the math function context
     234             : .  comb - how to combine the functions (addition, multiplication, division or composition)
     235             : .  f1   - first function
     236             : -  f2   - second function
     237             : 
     238             :    Level: intermediate
     239             : 
     240             : .seealso: FNCombineGetChildren()
     241             : @*/
     242          11 : PetscErrorCode FNCombineSetChildren(FN fn,FNCombineType comb,FN f1,FN f2)
     243             : {
     244          11 :   PetscFunctionBegin;
     245          11 :   PetscValidHeaderSpecific(fn,FN_CLASSID,1);
     246          33 :   PetscValidLogicalCollectiveEnum(fn,comb,2);
     247          11 :   PetscValidHeaderSpecific(f1,FN_CLASSID,3);
     248          11 :   PetscValidHeaderSpecific(f2,FN_CLASSID,4);
     249          11 :   PetscTryMethod(fn,"FNCombineSetChildren_C",(FN,FNCombineType,FN,FN),(fn,comb,f1,f2));
     250          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     251             : }
     252             : 
     253           3 : static PetscErrorCode FNCombineGetChildren_Combine(FN fn,FNCombineType *comb,FN *f1,FN *f2)
     254             : {
     255           3 :   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
     256             : 
     257           3 :   PetscFunctionBegin;
     258           3 :   if (comb) *comb = ctx->comb;
     259           3 :   if (f1) {
     260           3 :     if (!ctx->f1) PetscCall(FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f1));
     261           3 :     *f1 = ctx->f1;
     262             :   }
     263           3 :   if (f2) {
     264           3 :     if (!ctx->f2) PetscCall(FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f2));
     265           3 :     *f2 = ctx->f2;
     266             :   }
     267           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     268             : }
     269             : 
     270             : /*@
     271             :    FNCombineGetChildren - Gets the two child functions that constitute this
     272             :    combined function, and the way they are combined.
     273             : 
     274             :    Not Collective
     275             : 
     276             :    Input Parameter:
     277             : .  fn   - the math function context
     278             : 
     279             :    Output Parameters:
     280             : +  comb - how to combine the functions (addition, multiplication, division or composition)
     281             : .  f1   - first function
     282             : -  f2   - second function
     283             : 
     284             :    Level: intermediate
     285             : 
     286             : .seealso: FNCombineSetChildren()
     287             : @*/
     288           3 : PetscErrorCode FNCombineGetChildren(FN fn,FNCombineType *comb,FN *f1,FN *f2)
     289             : {
     290           3 :   PetscFunctionBegin;
     291           3 :   PetscValidHeaderSpecific(fn,FN_CLASSID,1);
     292           3 :   PetscUseMethod(fn,"FNCombineGetChildren_C",(FN,FNCombineType*,FN*,FN*),(fn,comb,f1,f2));
     293           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     294             : }
     295             : 
     296           4 : static PetscErrorCode FNDuplicate_Combine(FN fn,MPI_Comm comm,FN *newfn)
     297             : {
     298           4 :   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data,*ctx2 = (FN_COMBINE*)(*newfn)->data;
     299             : 
     300           4 :   PetscFunctionBegin;
     301           4 :   ctx2->comb = ctx->comb;
     302           4 :   PetscCall(FNDuplicate(ctx->f1,comm,&ctx2->f1));
     303           4 :   PetscCall(FNDuplicate(ctx->f2,comm,&ctx2->f2));
     304           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     305             : }
     306             : 
     307          15 : static PetscErrorCode FNDestroy_Combine(FN fn)
     308             : {
     309          15 :   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
     310             : 
     311          15 :   PetscFunctionBegin;
     312          15 :   PetscCall(FNDestroy(&ctx->f1));
     313          15 :   PetscCall(FNDestroy(&ctx->f2));
     314          15 :   PetscCall(PetscFree(fn->data));
     315          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",NULL));
     316          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",NULL));
     317          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     318             : }
     319             : 
     320          15 : SLEPC_EXTERN PetscErrorCode FNCreate_Combine(FN fn)
     321             : {
     322          15 :   FN_COMBINE     *ctx;
     323             : 
     324          15 :   PetscFunctionBegin;
     325          15 :   PetscCall(PetscNew(&ctx));
     326          15 :   fn->data = (void*)ctx;
     327             : 
     328          15 :   fn->ops->evaluatefunction          = FNEvaluateFunction_Combine;
     329          15 :   fn->ops->evaluatederivative        = FNEvaluateDerivative_Combine;
     330          15 :   fn->ops->evaluatefunctionmat[0]    = FNEvaluateFunctionMat_Combine;
     331          15 :   fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Combine;
     332             : #if defined(PETSC_HAVE_CUDA)
     333             :   fn->ops->evaluatefunctionmatcuda[0]    = FNEvaluateFunctionMat_Combine;
     334             :   fn->ops->evaluatefunctionmatveccuda[0] = FNEvaluateFunctionMatVec_Combine;
     335             : #endif
     336          15 :   fn->ops->view                      = FNView_Combine;
     337          15 :   fn->ops->duplicate                 = FNDuplicate_Combine;
     338          15 :   fn->ops->destroy                   = FNDestroy_Combine;
     339          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",FNCombineSetChildren_Combine));
     340          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",FNCombineGetChildren_Combine));
     341          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     342             : }