Actual source code: fncombine.c

slepc-main 2024-11-09
Report Typos and Errors
  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;

186:   PetscFunctionBegin;
187:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
188:   if (isascii) {
189:     switch (ctx->comb) {
190:       case FN_COMBINE_ADD:
191:         PetscCall(PetscViewerASCIIPrintf(viewer,"  two added functions f1+f2\n"));
192:         break;
193:       case FN_COMBINE_MULTIPLY:
194:         PetscCall(PetscViewerASCIIPrintf(viewer,"  two multiplied functions f1*f2\n"));
195:         break;
196:       case FN_COMBINE_DIVIDE:
197:         PetscCall(PetscViewerASCIIPrintf(viewer,"  a quotient of two functions f1/f2\n"));
198:         break;
199:       case FN_COMBINE_COMPOSE:
200:         PetscCall(PetscViewerASCIIPrintf(viewer,"  two composed functions f2(f1(.))\n"));
201:         break;
202:     }
203:     PetscCall(PetscViewerASCIIPushTab(viewer));
204:     PetscCall(FNView(ctx->f1,viewer));
205:     PetscCall(FNView(ctx->f2,viewer));
206:     PetscCall(PetscViewerASCIIPopTab(viewer));
207:   }
208:   PetscFunctionReturn(PETSC_SUCCESS);
209: }

211: static PetscErrorCode FNCombineSetChildren_Combine(FN fn,FNCombineType comb,FN f1,FN f2)
212: {
213:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;

215:   PetscFunctionBegin;
216:   ctx->comb = comb;
217:   PetscCall(PetscObjectReference((PetscObject)f1));
218:   PetscCall(FNDestroy(&ctx->f1));
219:   ctx->f1 = f1;
220:   PetscCall(PetscObjectReference((PetscObject)f2));
221:   PetscCall(FNDestroy(&ctx->f2));
222:   ctx->f2 = f2;
223:   PetscFunctionReturn(PETSC_SUCCESS);
224: }

226: /*@
227:    FNCombineSetChildren - Sets the two child functions that constitute this
228:    combined function, and the way they must be combined.

230:    Logically Collective

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

238:    Level: intermediate

240: .seealso: FNCombineGetChildren()
241: @*/
242: PetscErrorCode FNCombineSetChildren(FN fn,FNCombineType comb,FN f1,FN f2)
243: {
244:   PetscFunctionBegin;
249:   PetscTryMethod(fn,"FNCombineSetChildren_C",(FN,FNCombineType,FN,FN),(fn,comb,f1,f2));
250:   PetscFunctionReturn(PETSC_SUCCESS);
251: }

253: static PetscErrorCode FNCombineGetChildren_Combine(FN fn,FNCombineType *comb,FN *f1,FN *f2)
254: {
255:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;

257:   PetscFunctionBegin;
258:   if (comb) *comb = ctx->comb;
259:   if (f1) {
260:     if (!ctx->f1) PetscCall(FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f1));
261:     *f1 = ctx->f1;
262:   }
263:   if (f2) {
264:     if (!ctx->f2) PetscCall(FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f2));
265:     *f2 = ctx->f2;
266:   }
267:   PetscFunctionReturn(PETSC_SUCCESS);
268: }

270: /*@
271:    FNCombineGetChildren - Gets the two child functions that constitute this
272:    combined function, and the way they are combined.

274:    Not Collective

276:    Input Parameter:
277: .  fn   - the math function context

279:    Output Parameters:
280: +  comb - how to combine the functions (addition, multiplication, division or composition)
281: .  f1   - first function
282: -  f2   - second function

284:    Level: intermediate

286: .seealso: FNCombineSetChildren()
287: @*/
288: PetscErrorCode FNCombineGetChildren(FN fn,FNCombineType *comb,FN *f1,FN *f2)
289: {
290:   PetscFunctionBegin;
292:   PetscUseMethod(fn,"FNCombineGetChildren_C",(FN,FNCombineType*,FN*,FN*),(fn,comb,f1,f2));
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }

296: static PetscErrorCode FNDuplicate_Combine(FN fn,MPI_Comm comm,FN *newfn)
297: {
298:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data,*ctx2 = (FN_COMBINE*)(*newfn)->data;

300:   PetscFunctionBegin;
301:   ctx2->comb = ctx->comb;
302:   PetscCall(FNDuplicate(ctx->f1,comm,&ctx2->f1));
303:   PetscCall(FNDuplicate(ctx->f2,comm,&ctx2->f2));
304:   PetscFunctionReturn(PETSC_SUCCESS);
305: }

307: static PetscErrorCode FNDestroy_Combine(FN fn)
308: {
309:   FN_COMBINE     *ctx = (FN_COMBINE*)fn->data;

311:   PetscFunctionBegin;
312:   PetscCall(FNDestroy(&ctx->f1));
313:   PetscCall(FNDestroy(&ctx->f2));
314:   PetscCall(PetscFree(fn->data));
315:   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",NULL));
316:   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",NULL));
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }

320: SLEPC_EXTERN PetscErrorCode FNCreate_Combine(FN fn)
321: {
322:   FN_COMBINE     *ctx;

324:   PetscFunctionBegin;
325:   PetscCall(PetscNew(&ctx));
326:   fn->data = (void*)ctx;

328:   fn->ops->evaluatefunction          = FNEvaluateFunction_Combine;
329:   fn->ops->evaluatederivative        = FNEvaluateDerivative_Combine;
330:   fn->ops->evaluatefunctionmat[0]    = FNEvaluateFunctionMat_Combine;
331:   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:   fn->ops->view                      = FNView_Combine;
337:   fn->ops->duplicate                 = FNDuplicate_Combine;
338:   fn->ops->destroy                   = FNDestroy_Combine;
339:   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",FNCombineSetChildren_Combine));
340:   PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",FNCombineGetChildren_Combine));
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: }