Actual source code: fncombine.c
slepc-3.17.2 2022-08-09
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>
21: #include <slepcblaslapack.h>
23: typedef struct {
24: FN f1,f2; /* functions */
25: FNCombineType comb; /* how the functions are combined */
26: } FN_COMBINE;
28: PetscErrorCode FNEvaluateFunction_Combine(FN fn,PetscScalar x,PetscScalar *y)
29: {
30: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
31: PetscScalar a,b;
33: FNEvaluateFunction(ctx->f1,x,&a);
34: switch (ctx->comb) {
35: case FN_COMBINE_ADD:
36: FNEvaluateFunction(ctx->f2,x,&b);
37: *y = a+b;
38: break;
39: case FN_COMBINE_MULTIPLY:
40: FNEvaluateFunction(ctx->f2,x,&b);
41: *y = a*b;
42: break;
43: case FN_COMBINE_DIVIDE:
44: FNEvaluateFunction(ctx->f2,x,&b);
46: *y = a/b;
47: break;
48: case FN_COMBINE_COMPOSE:
49: FNEvaluateFunction(ctx->f2,a,y);
50: break;
51: }
52: PetscFunctionReturn(0);
53: }
55: 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: switch (ctx->comb) {
61: case FN_COMBINE_ADD:
62: FNEvaluateDerivative(ctx->f1,x,&ap);
63: FNEvaluateDerivative(ctx->f2,x,&bp);
64: *yp = ap+bp;
65: break;
66: case FN_COMBINE_MULTIPLY:
67: FNEvaluateDerivative(ctx->f1,x,&ap);
68: FNEvaluateDerivative(ctx->f2,x,&bp);
69: FNEvaluateFunction(ctx->f1,x,&a);
70: FNEvaluateFunction(ctx->f2,x,&b);
71: *yp = ap*b+a*bp;
72: break;
73: case FN_COMBINE_DIVIDE:
74: FNEvaluateDerivative(ctx->f1,x,&ap);
75: FNEvaluateDerivative(ctx->f2,x,&bp);
76: FNEvaluateFunction(ctx->f1,x,&a);
77: FNEvaluateFunction(ctx->f2,x,&b);
79: *yp = (ap*b-a*bp)/(b*b);
80: break;
81: case FN_COMBINE_COMPOSE:
82: FNEvaluateFunction(ctx->f1,x,&a);
83: FNEvaluateDerivative(ctx->f1,x,&ap);
84: FNEvaluateDerivative(ctx->f2,a,yp);
85: *yp *= ap;
86: break;
87: }
88: PetscFunctionReturn(0);
89: }
91: PetscErrorCode FNEvaluateFunctionMat_Combine(FN fn,Mat A,Mat B)
92: {
93: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
94: PetscScalar *Ba,*Wa,one=1.0,zero=0.0;
95: const PetscScalar *Za;
96: PetscBLASInt n,ld,ld2,inc=1,*ipiv,info;
97: PetscInt m;
98: Mat W,Z;
100: FN_AllocateWorkMat(fn,A,&W);
101: MatGetSize(A,&m,NULL);
102: PetscBLASIntCast(m,&n);
103: ld = n;
104: ld2 = ld*ld;
106: switch (ctx->comb) {
107: case FN_COMBINE_ADD:
108: FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE);
109: FNEvaluateFunctionMat_Private(ctx->f2,A,B,PETSC_FALSE);
110: MatDenseGetArray(B,&Ba);
111: MatDenseGetArray(W,&Wa);
112: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&ld2,&one,Wa,&inc,Ba,&inc));
113: PetscLogFlops(1.0*n*n);
114: MatDenseRestoreArray(B,&Ba);
115: MatDenseRestoreArray(W,&Wa);
116: break;
117: case FN_COMBINE_MULTIPLY:
118: FN_AllocateWorkMat(fn,A,&Z);
119: FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE);
120: FNEvaluateFunctionMat_Private(ctx->f2,A,Z,PETSC_FALSE);
121: MatDenseGetArray(B,&Ba);
122: MatDenseGetArray(W,&Wa);
123: MatDenseGetArrayRead(Z,&Za);
124: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,Wa,&ld,Za,&ld,&zero,Ba,&ld));
125: PetscLogFlops(2.0*n*n*n);
126: MatDenseRestoreArray(B,&Ba);
127: MatDenseRestoreArray(W,&Wa);
128: MatDenseRestoreArrayRead(Z,&Za);
129: FN_FreeWorkMat(fn,&Z);
130: break;
131: case FN_COMBINE_DIVIDE:
132: FNEvaluateFunctionMat_Private(ctx->f2,A,W,PETSC_FALSE);
133: FNEvaluateFunctionMat_Private(ctx->f1,A,B,PETSC_FALSE);
134: PetscMalloc1(ld,&ipiv);
135: MatDenseGetArray(B,&Ba);
136: MatDenseGetArray(W,&Wa);
137: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&n,Wa,&ld,ipiv,Ba,&ld,&info));
138: SlepcCheckLapackInfo("gesv",info);
139: PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n*n);
140: MatDenseRestoreArray(B,&Ba);
141: MatDenseRestoreArray(W,&Wa);
142: PetscFree(ipiv);
143: break;
144: case FN_COMBINE_COMPOSE:
145: FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE);
146: FNEvaluateFunctionMat_Private(ctx->f2,W,B,PETSC_FALSE);
147: break;
148: }
150: FN_FreeWorkMat(fn,&W);
151: PetscFunctionReturn(0);
152: }
154: PetscErrorCode FNEvaluateFunctionMatVec_Combine(FN fn,Mat A,Vec v)
155: {
156: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
157: PetscScalar *va,*Za;
158: PetscBLASInt n,ld,*ipiv,info,one=1;
159: PetscInt m;
160: Mat Z;
161: Vec w;
163: MatGetSize(A,&m,NULL);
164: PetscBLASIntCast(m,&n);
165: ld = n;
167: switch (ctx->comb) {
168: case FN_COMBINE_ADD:
169: VecDuplicate(v,&w);
170: FNEvaluateFunctionMatVec(ctx->f1,A,w);
171: FNEvaluateFunctionMatVec(ctx->f2,A,v);
172: VecAXPY(v,1.0,w);
173: VecDestroy(&w);
174: break;
175: case FN_COMBINE_MULTIPLY:
176: VecDuplicate(v,&w);
177: FN_AllocateWorkMat(fn,A,&Z);
178: FNEvaluateFunctionMat_Private(ctx->f1,A,Z,PETSC_FALSE);
179: FNEvaluateFunctionMatVec_Private(ctx->f2,A,w,PETSC_FALSE);
180: MatMult(Z,w,v);
181: FN_FreeWorkMat(fn,&Z);
182: VecDestroy(&w);
183: break;
184: case FN_COMBINE_DIVIDE:
185: VecDuplicate(v,&w);
186: FN_AllocateWorkMat(fn,A,&Z);
187: FNEvaluateFunctionMat_Private(ctx->f2,A,Z,PETSC_FALSE);
188: FNEvaluateFunctionMatVec_Private(ctx->f1,A,v,PETSC_FALSE);
189: PetscMalloc1(ld,&ipiv);
190: MatDenseGetArray(Z,&Za);
191: VecGetArray(v,&va);
192: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&one,Za,&ld,ipiv,va,&ld,&info));
193: SlepcCheckLapackInfo("gesv",info);
194: PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n);
195: VecRestoreArray(v,&va);
196: MatDenseRestoreArray(Z,&Za);
197: PetscFree(ipiv);
198: FN_FreeWorkMat(fn,&Z);
199: VecDestroy(&w);
200: break;
201: case FN_COMBINE_COMPOSE:
202: FN_AllocateWorkMat(fn,A,&Z);
203: FNEvaluateFunctionMat_Private(ctx->f1,A,Z,PETSC_FALSE);
204: FNEvaluateFunctionMatVec_Private(ctx->f2,Z,v,PETSC_FALSE);
205: FN_FreeWorkMat(fn,&Z);
206: break;
207: }
208: PetscFunctionReturn(0);
209: }
211: PetscErrorCode FNView_Combine(FN fn,PetscViewer viewer)
212: {
213: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
214: PetscBool isascii;
216: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
217: if (isascii) {
218: switch (ctx->comb) {
219: case FN_COMBINE_ADD:
220: PetscViewerASCIIPrintf(viewer," Two added functions f1+f2\n");
221: break;
222: case FN_COMBINE_MULTIPLY:
223: PetscViewerASCIIPrintf(viewer," Two multiplied functions f1*f2\n");
224: break;
225: case FN_COMBINE_DIVIDE:
226: PetscViewerASCIIPrintf(viewer," A quotient of two functions f1/f2\n");
227: break;
228: case FN_COMBINE_COMPOSE:
229: PetscViewerASCIIPrintf(viewer," Two composed functions f2(f1(.))\n");
230: break;
231: }
232: PetscViewerASCIIPushTab(viewer);
233: FNView(ctx->f1,viewer);
234: FNView(ctx->f2,viewer);
235: PetscViewerASCIIPopTab(viewer);
236: }
237: PetscFunctionReturn(0);
238: }
240: static PetscErrorCode FNCombineSetChildren_Combine(FN fn,FNCombineType comb,FN f1,FN f2)
241: {
242: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
244: ctx->comb = comb;
245: PetscObjectReference((PetscObject)f1);
246: FNDestroy(&ctx->f1);
247: ctx->f1 = f1;
248: PetscLogObjectParent((PetscObject)fn,(PetscObject)ctx->f1);
249: PetscObjectReference((PetscObject)f2);
250: FNDestroy(&ctx->f2);
251: ctx->f2 = f2;
252: PetscLogObjectParent((PetscObject)fn,(PetscObject)ctx->f2);
253: PetscFunctionReturn(0);
254: }
256: /*@
257: FNCombineSetChildren - Sets the two child functions that constitute this
258: combined function, and the way they must be combined.
260: Logically Collective on fn
262: Input Parameters:
263: + fn - the math function context
264: . comb - how to combine the functions (addition, multiplication, division or composition)
265: . f1 - first function
266: - f2 - second function
268: Level: intermediate
270: .seealso: FNCombineGetChildren()
271: @*/
272: PetscErrorCode FNCombineSetChildren(FN fn,FNCombineType comb,FN f1,FN f2)
273: {
278: PetscTryMethod(fn,"FNCombineSetChildren_C",(FN,FNCombineType,FN,FN),(fn,comb,f1,f2));
279: PetscFunctionReturn(0);
280: }
282: static PetscErrorCode FNCombineGetChildren_Combine(FN fn,FNCombineType *comb,FN *f1,FN *f2)
283: {
284: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
286: if (comb) *comb = ctx->comb;
287: if (f1) {
288: if (!ctx->f1) {
289: FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f1);
290: PetscLogObjectParent((PetscObject)fn,(PetscObject)ctx->f1);
291: }
292: *f1 = ctx->f1;
293: }
294: if (f2) {
295: if (!ctx->f2) {
296: FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f2);
297: PetscLogObjectParent((PetscObject)fn,(PetscObject)ctx->f2);
298: }
299: *f2 = ctx->f2;
300: }
301: PetscFunctionReturn(0);
302: }
304: /*@
305: FNCombineGetChildren - Gets the two child functions that constitute this
306: combined function, and the way they are combined.
308: Not Collective
310: Input Parameter:
311: . fn - the math function context
313: Output Parameters:
314: + comb - how to combine the functions (addition, multiplication, division or composition)
315: . f1 - first function
316: - f2 - second function
318: Level: intermediate
320: .seealso: FNCombineSetChildren()
321: @*/
322: PetscErrorCode FNCombineGetChildren(FN fn,FNCombineType *comb,FN *f1,FN *f2)
323: {
325: PetscUseMethod(fn,"FNCombineGetChildren_C",(FN,FNCombineType*,FN*,FN*),(fn,comb,f1,f2));
326: PetscFunctionReturn(0);
327: }
329: PetscErrorCode FNDuplicate_Combine(FN fn,MPI_Comm comm,FN *newfn)
330: {
331: FN_COMBINE *ctx = (FN_COMBINE*)fn->data,*ctx2 = (FN_COMBINE*)(*newfn)->data;
333: ctx2->comb = ctx->comb;
334: FNDuplicate(ctx->f1,comm,&ctx2->f1);
335: FNDuplicate(ctx->f2,comm,&ctx2->f2);
336: PetscFunctionReturn(0);
337: }
339: PetscErrorCode FNDestroy_Combine(FN fn)
340: {
341: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
343: FNDestroy(&ctx->f1);
344: FNDestroy(&ctx->f2);
345: PetscFree(fn->data);
346: PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",NULL);
347: PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",NULL);
348: PetscFunctionReturn(0);
349: }
351: SLEPC_EXTERN PetscErrorCode FNCreate_Combine(FN fn)
352: {
353: FN_COMBINE *ctx;
355: PetscNewLog(fn,&ctx);
356: fn->data = (void*)ctx;
358: fn->ops->evaluatefunction = FNEvaluateFunction_Combine;
359: fn->ops->evaluatederivative = FNEvaluateDerivative_Combine;
360: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Combine;
361: fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Combine;
362: fn->ops->view = FNView_Combine;
363: fn->ops->duplicate = FNDuplicate_Combine;
364: fn->ops->destroy = FNDestroy_Combine;
365: PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",FNCombineSetChildren_Combine);
366: PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",FNCombineGetChildren_Combine);
367: PetscFunctionReturn(0);
368: }