Actual source code: fncombine.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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:    A function that is obtained by combining two other functions (either by
 12:    addition, multiplication, division or composition)

 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: */

 20: #include <slepc/private/fnimpl.h>

 22: typedef struct {
 23:   FN            f1,f2;    /* functions */
 24:   FNCombineType comb;     /* how the functions are combined */
 25: } FN_COMBINE;

 27: static PetscErrorCode FNEvaluateFunction_Combine(FN fn,PetscScalar x,PetscScalar *y)
 28: {
 29:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
 30:   PetscScalar    a,b;

 32:   PetscFunctionBegin;
 33:   PetscCall(FNEvaluateFunction(ctx->f1,x,&a));
 34:   switch (ctx->comb) {
 35:     case FN_COMBINE_ADD:
 36:       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
 37:       *y = a+b;
 38:       break;
 39:     case FN_COMBINE_MULTIPLY:
 40:       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
 41:       *y = a*b;
 42:       break;
 43:     case FN_COMBINE_DIVIDE:
 44:       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
 45:       PetscCheck(b!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Function not defined in the requested value");
 46:       *y = a/b;
 47:       break;
 48:     case FN_COMBINE_COMPOSE:
 49:       PetscCall(FNEvaluateFunction(ctx->f2,a,y));
 50:       break;
 51:   }
 52:   PetscFunctionReturn(PETSC_SUCCESS);
 53: }

 55: static PetscErrorCode FNEvaluateDerivative_Combine(FN fn,PetscScalar x,PetscScalar *yp)
 56: {
 57:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
 58:   PetscScalar    a,b,ap,bp;

 60:   PetscFunctionBegin;
 61:   switch (ctx->comb) {
 62:     case FN_COMBINE_ADD:
 63:       PetscCall(FNEvaluateDerivative(ctx->f1,x,&ap));
 64:       PetscCall(FNEvaluateDerivative(ctx->f2,x,&bp));
 65:       *yp = ap+bp;
 66:       break;
 67:     case FN_COMBINE_MULTIPLY:
 68:       PetscCall(FNEvaluateDerivative(ctx->f1,x,&ap));
 69:       PetscCall(FNEvaluateDerivative(ctx->f2,x,&bp));
 70:       PetscCall(FNEvaluateFunction(ctx->f1,x,&a));
 71:       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
 72:       *yp = ap*b+a*bp;
 73:       break;
 74:     case FN_COMBINE_DIVIDE:
 75:       PetscCall(FNEvaluateDerivative(ctx->f1,x,&ap));
 76:       PetscCall(FNEvaluateDerivative(ctx->f2,x,&bp));
 77:       PetscCall(FNEvaluateFunction(ctx->f1,x,&a));
 78:       PetscCall(FNEvaluateFunction(ctx->f2,x,&b));
 79:       PetscCheck(b!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
 80:       *yp = (ap*b-a*bp)/(b*b);
 81:       break;
 82:     case FN_COMBINE_COMPOSE:
 83:       PetscCall(FNEvaluateFunction(ctx->f1,x,&a));
 84:       PetscCall(FNEvaluateDerivative(ctx->f1,x,&ap));
 85:       PetscCall(FNEvaluateDerivative(ctx->f2,a,yp));
 86:       *yp *= ap;
 87:       break;
 88:   }
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: static PetscErrorCode FNEvaluateFunctionMat_Combine(FN fn,Mat A,Mat B)
 93: {
 94:   FN_COMBINE   *ctx = (FN_COMBINE*)fn->data;
 95:   Mat          W,Z,F;
 96:   PetscBool    iscuda;

 98:   PetscFunctionBegin;
 99:   PetscCall(FN_AllocateWorkMat(fn,A,&W));
100:   switch (ctx->comb) {
101:     case FN_COMBINE_ADD:
102:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE));
103:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,A,B,PETSC_FALSE));
104:       PetscCall(MatAXPY(B,1.0,W,SAME_NONZERO_PATTERN));
105:       break;
106:     case FN_COMBINE_MULTIPLY:
107:       PetscCall(FN_AllocateWorkMat(fn,A,&Z));
108:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE));
109:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,A,Z,PETSC_FALSE));
110:       PetscCall(MatMatMult(W,Z,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B));
111:       PetscCall(FN_FreeWorkMat(fn,&Z));
112:       break;
113:     case FN_COMBINE_DIVIDE:
114:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,A,W,PETSC_FALSE));
115:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,B,PETSC_FALSE));
116:       PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
117:       PetscCall(MatGetFactor(W,iscuda?MATSOLVERCUDA:MATSOLVERPETSC,MAT_FACTOR_LU,&F));
118:       PetscCall(MatLUFactorSymbolic(F,W,NULL,NULL,NULL));
119:       PetscCall(MatLUFactorNumeric(F,W,NULL));
120:       PetscCall(MatMatSolve(F,B,B));
121:       PetscCall(MatDestroy(&F));
122:       break;
123:     case FN_COMBINE_COMPOSE:
124:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE));
125:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,W,B,PETSC_FALSE));
126:       break;
127:   }
128:   PetscCall(FN_FreeWorkMat(fn,&W));
129:   PetscFunctionReturn(PETSC_SUCCESS);
130: }

132: static PetscErrorCode FNEvaluateFunctionMatVec_Combine(FN fn,Mat A,Vec v)
133: {
134:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
135:   PetscBool      iscuda;
136:   Mat            Z,F;
137:   Vec            w;

139:   PetscFunctionBegin;
140:   switch (ctx->comb) {
141:     case FN_COMBINE_ADD:
142:       PetscCall(VecDuplicate(v,&w));
143:       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f1,A,w,PETSC_FALSE));
144:       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f2,A,v,PETSC_FALSE));
145:       PetscCall(VecAXPY(v,1.0,w));
146:       PetscCall(VecDestroy(&w));
147:       break;
148:     case FN_COMBINE_MULTIPLY:
149:       PetscCall(VecDuplicate(v,&w));
150:       PetscCall(FN_AllocateWorkMat(fn,A,&Z));
151:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,Z,PETSC_FALSE));
152:       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f2,A,w,PETSC_FALSE));
153:       PetscCall(MatMult(Z,w,v));
154:       PetscCall(FN_FreeWorkMat(fn,&Z));
155:       PetscCall(VecDestroy(&w));
156:       break;
157:     case FN_COMBINE_DIVIDE:
158:       PetscCall(VecDuplicate(v,&w));
159:       PetscCall(FN_AllocateWorkMat(fn,A,&Z));
160:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f2,A,Z,PETSC_FALSE));
161:       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f1,A,w,PETSC_FALSE));
162:       PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
163:       PetscCall(MatGetFactor(Z,iscuda?MATSOLVERCUDA:MATSOLVERPETSC,MAT_FACTOR_LU,&F));
164:       PetscCall(MatLUFactorSymbolic(F,Z,NULL,NULL,NULL));
165:       PetscCall(MatLUFactorNumeric(F,Z,NULL));
166:       PetscCall(MatSolve(F,w,v));
167:       PetscCall(MatDestroy(&F));
168:       PetscCall(FN_FreeWorkMat(fn,&Z));
169:       PetscCall(VecDestroy(&w));
170:       break;
171:     case FN_COMBINE_COMPOSE:
172:       PetscCall(FN_AllocateWorkMat(fn,A,&Z));
173:       PetscCall(FNEvaluateFunctionMat_Private(ctx->f1,A,Z,PETSC_FALSE));
174:       PetscCall(FNEvaluateFunctionMatVec_Private(ctx->f2,Z,v,PETSC_FALSE));
175:       PetscCall(FN_FreeWorkMat(fn,&Z));
176:       break;
177:   }
178:   PetscFunctionReturn(PETSC_SUCCESS);
179: }

181: static PetscErrorCode FNView_Combine(FN fn,PetscViewer viewer)
182: {
183:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;
184:   PetscBool      isascii;
185:   char           str[50];

187:   PetscFunctionBegin;
188:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
189:   if (isascii) {
190:     if (fn->alpha!=(PetscScalar)1.0 || fn->beta!=(PetscScalar)1.0) {
191:       PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
192:       PetscCall(PetscViewerASCIIPrintf(viewer,"  scaling factors: inner %s",str));
193:       PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_TRUE));
194:       PetscCall(PetscViewerASCIIPrintf(viewer,"outer %s\n",str));
195:     }
196:     switch (ctx->comb) {
197:       case FN_COMBINE_ADD:
198:         PetscCall(PetscViewerASCIIPrintf(viewer,"  two added functions f1+f2\n"));
199:         break;
200:       case FN_COMBINE_MULTIPLY:
201:         PetscCall(PetscViewerASCIIPrintf(viewer,"  two multiplied functions f1*f2\n"));
202:         break;
203:       case FN_COMBINE_DIVIDE:
204:         PetscCall(PetscViewerASCIIPrintf(viewer,"  a quotient of two functions f1/f2\n"));
205:         break;
206:       case FN_COMBINE_COMPOSE:
207:         PetscCall(PetscViewerASCIIPrintf(viewer,"  two composed functions f2(f1(.))\n"));
208:         break;
209:     }
210:     PetscCall(PetscViewerASCIIPushTab(viewer));
211:     PetscCall(FNView(ctx->f1,viewer));
212:     PetscCall(FNView(ctx->f2,viewer));
213:     PetscCall(PetscViewerASCIIPopTab(viewer));
214:   }
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: }

218: static PetscErrorCode FNCombineSetChildren_Combine(FN fn,FNCombineType comb,FN f1,FN f2)
219: {
220:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;

222:   PetscFunctionBegin;
223:   ctx->comb = comb;
224:   PetscCall(PetscObjectReference((PetscObject)f1));
225:   PetscCall(FNDestroy(&ctx->f1));
226:   ctx->f1 = f1;
227:   PetscCall(PetscObjectReference((PetscObject)f2));
228:   PetscCall(FNDestroy(&ctx->f2));
229:   ctx->f2 = f2;
230:   PetscFunctionReturn(PETSC_SUCCESS);
231: }

233: /*@
234:    FNCombineSetChildren - Sets the two child functions that constitute this
235:    combined function, and the way they must be combined.

237:    Logically Collective

239:    Input Parameters:
240: +  fn   - the math function context
241: .  comb - how to combine the functions (addition, multiplication, division or composition)
242: .  f1   - first function
243: -  f2   - second function

245:    Level: intermediate

247: .seealso: FNCombineGetChildren()
248: @*/
249: PetscErrorCode FNCombineSetChildren(FN fn,FNCombineType comb,FN f1,FN f2)
250: {
251:   PetscFunctionBegin;
256:   PetscTryMethod(fn,"FNCombineSetChildren_C",(FN,FNCombineType,FN,FN),(fn,comb,f1,f2));
257:   PetscFunctionReturn(PETSC_SUCCESS);
258: }

260: static PetscErrorCode FNCombineGetChildren_Combine(FN fn,FNCombineType *comb,FN *f1,FN *f2)
261: {
262:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;

264:   PetscFunctionBegin;
265:   if (comb) *comb = ctx->comb;
266:   if (f1) {
267:     if (!ctx->f1) PetscCall(FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f1));
268:     *f1 = ctx->f1;
269:   }
270:   if (f2) {
271:     if (!ctx->f2) PetscCall(FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f2));
272:     *f2 = ctx->f2;
273:   }
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }

277: /*@
278:    FNCombineGetChildren - Gets the two child functions that constitute this
279:    combined function, and the way they are combined.

281:    Not Collective

283:    Input Parameter:
284: .  fn   - the math function context

286:    Output Parameters:
287: +  comb - how to combine the functions (addition, multiplication, division or composition)
288: .  f1   - first function
289: -  f2   - second function

291:    Level: intermediate

293: .seealso: FNCombineSetChildren()
294: @*/
295: PetscErrorCode FNCombineGetChildren(FN fn,FNCombineType *comb,FN *f1,FN *f2)
296: {
297:   PetscFunctionBegin;
299:   PetscUseMethod(fn,"FNCombineGetChildren_C",(FN,FNCombineType*,FN*,FN*),(fn,comb,f1,f2));
300:   PetscFunctionReturn(PETSC_SUCCESS);
301: }

303: static PetscErrorCode FNDuplicate_Combine(FN fn,MPI_Comm comm,FN *newfn)
304: {
305:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data,*ctx2 = (FN_COMBINE*)(*newfn)->data;

307:   PetscFunctionBegin;
308:   ctx2->comb = ctx->comb;
309:   PetscCall(FNDuplicate(ctx->f1,comm,&ctx2->f1));
310:   PetscCall(FNDuplicate(ctx->f2,comm,&ctx2->f2));
311:   PetscFunctionReturn(PETSC_SUCCESS);
312: }

314: static PetscErrorCode FNDestroy_Combine(FN fn)
315: {
316:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;

318:   PetscFunctionBegin;
319:   PetscCall(FNDestroy(&ctx->f1));
320:   PetscCall(FNDestroy(&ctx->f2));
321:   PetscCall(PetscFree(fn->data));
322:   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",NULL));
323:   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",NULL));
324:   PetscFunctionReturn(PETSC_SUCCESS);
325: }

327: SLEPC_EXTERN PetscErrorCode FNCreate_Combine(FN fn)
328: {
329:   FN_COMBINE     *ctx;

331:   PetscFunctionBegin;
332:   PetscCall(PetscNew(&ctx));
333:   fn->data = (void*)ctx;

335:   fn->ops->evaluatefunction          = FNEvaluateFunction_Combine;
336:   fn->ops->evaluatederivative        = FNEvaluateDerivative_Combine;
337:   fn->ops->evaluatefunctionmat[0]    = FNEvaluateFunctionMat_Combine;
338:   fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Combine;
339: #if defined(PETSC_HAVE_CUDA)
340:   fn->ops->evaluatefunctionmatcuda[0]    = FNEvaluateFunctionMat_Combine;
341:   fn->ops->evaluatefunctionmatveccuda[0] = FNEvaluateFunctionMatVec_Combine;
342: #endif
343:   fn->ops->view                      = FNView_Combine;
344:   fn->ops->duplicate                 = FNDuplicate_Combine;
345:   fn->ops->destroy                   = FNDestroy_Combine;
346:   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",FNCombineSetChildren_Combine));
347:   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",FNCombineGetChildren_Combine));
348:   PetscFunctionReturn(PETSC_SUCCESS);
349: }