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: [](sec:fn), `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: [](sec:fn), `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: /*MC
328: FNCOMBINE - FNCOMBINE = "combine" - A mathematical function that is defined
329: by combining two previously defined functions, e.g., by addition. Use
330: `FNCombineSetChildren()` to specify which functions must be combined and how.
332: Level: beginner
334: .seealso: [](sec:fn), `FN`, `FNType`, `FNSetType()`, `FNCombineSetChildren()`
335: M*/
337: SLEPC_EXTERN PetscErrorCode FNCreate_Combine(FN fn)
338: {
339: FN_COMBINE *ctx;
341: PetscFunctionBegin;
342: PetscCall(PetscNew(&ctx));
343: fn->data = (void*)ctx;
345: fn->ops->evaluatefunction = FNEvaluateFunction_Combine;
346: fn->ops->evaluatederivative = FNEvaluateDerivative_Combine;
347: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Combine;
348: fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Combine;
349: #if defined(PETSC_HAVE_CUDA)
350: fn->ops->evaluatefunctionmatcuda[0] = FNEvaluateFunctionMat_Combine;
351: fn->ops->evaluatefunctionmatveccuda[0] = FNEvaluateFunctionMatVec_Combine;
352: #endif
353: fn->ops->view = FNView_Combine;
354: fn->ops->duplicate = FNDuplicate_Combine;
355: fn->ops->destroy = FNDestroy_Combine;
356: PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",FNCombineSetChildren_Combine));
357: PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",FNCombineGetChildren_Combine));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }