Actual source code: fncombine.c
slepc-3.22.1 2024-10-28
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: }