Actual source code: fncombine.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcfn.h" I*/
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: {
31: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
32: PetscScalar a,b;
35: FNEvaluateFunction(ctx->f1,x,&a);
36: switch (ctx->comb) {
37: case FN_COMBINE_ADD:
38: FNEvaluateFunction(ctx->f2,x,&b);
39: *y = a+b;
40: break;
41: case FN_COMBINE_MULTIPLY:
42: FNEvaluateFunction(ctx->f2,x,&b);
43: *y = a*b;
44: break;
45: case FN_COMBINE_DIVIDE:
46: FNEvaluateFunction(ctx->f2,x,&b);
47: if (b==0.0) SETERRQ(PETSC_COMM_SELF,1,"Function not defined in the requested value");
48: *y = a/b;
49: break;
50: case FN_COMBINE_COMPOSE:
51: FNEvaluateFunction(ctx->f2,a,y);
52: break;
53: }
54: return(0);
55: }
57: PetscErrorCode FNEvaluateDerivative_Combine(FN fn,PetscScalar x,PetscScalar *yp)
58: {
60: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
61: PetscScalar a,b,ap,bp;
64: switch (ctx->comb) {
65: case FN_COMBINE_ADD:
66: FNEvaluateDerivative(ctx->f1,x,&ap);
67: FNEvaluateDerivative(ctx->f2,x,&bp);
68: *yp = ap+bp;
69: break;
70: case FN_COMBINE_MULTIPLY:
71: FNEvaluateDerivative(ctx->f1,x,&ap);
72: FNEvaluateDerivative(ctx->f2,x,&bp);
73: FNEvaluateFunction(ctx->f1,x,&a);
74: FNEvaluateFunction(ctx->f2,x,&b);
75: *yp = ap*b+a*bp;
76: break;
77: case FN_COMBINE_DIVIDE:
78: FNEvaluateDerivative(ctx->f1,x,&ap);
79: FNEvaluateDerivative(ctx->f2,x,&bp);
80: FNEvaluateFunction(ctx->f1,x,&a);
81: FNEvaluateFunction(ctx->f2,x,&b);
82: if (b==0.0) SETERRQ(PETSC_COMM_SELF,1,"Derivative not defined in the requested value");
83: *yp = (ap*b-a*bp)/(b*b);
84: break;
85: case FN_COMBINE_COMPOSE:
86: FNEvaluateFunction(ctx->f1,x,&a);
87: FNEvaluateDerivative(ctx->f1,x,&ap);
88: FNEvaluateDerivative(ctx->f2,a,yp);
89: *yp *= ap;
90: break;
91: }
92: return(0);
93: }
95: PetscErrorCode FNEvaluateFunctionMat_Combine(FN fn,Mat A,Mat B)
96: {
97: #if defined(PETSC_MISSING_LAPACK_GESV)
99: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESV - Lapack routines are unavailable");
100: #else
102: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
103: PetscScalar *Aa,*Ba,*Wa,*Za,one=1.0,zero=0.0;
104: PetscBLASInt n,ld,ld2,inc=1,*ipiv,info;
105: PetscInt m;
106: Mat W,Z;
109: FN_AllocateWorkMat(fn,A,&W);
110: MatDenseGetArray(A,&Aa);
111: MatDenseGetArray(B,&Ba);
112: MatDenseGetArray(W,&Wa);
113: MatGetSize(A,&m,NULL);
114: PetscBLASIntCast(m,&n);
115: ld = n;
116: ld2 = ld*ld;
118: switch (ctx->comb) {
119: case FN_COMBINE_ADD:
120: FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE);
121: FNEvaluateFunctionMat_Private(ctx->f2,A,B,PETSC_FALSE);
122: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&ld2,&one,Wa,&inc,Ba,&inc));
123: PetscLogFlops(1.0*n*n);
124: break;
125: case FN_COMBINE_MULTIPLY:
126: FN_AllocateWorkMat(fn,A,&Z);
127: MatDenseGetArray(Z,&Za);
128: FNEvaluateFunctionMat_Private(ctx->f1,A,W,PETSC_FALSE);
129: FNEvaluateFunctionMat_Private(ctx->f2,A,Z,PETSC_FALSE);
130: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,Wa,&ld,Za,&ld,&zero,Ba,&ld));
131: PetscLogFlops(2.0*n*n*n);
132: MatDenseRestoreArray(Z,&Za);
133: FN_FreeWorkMat(fn,&Z);
134: break;
135: case FN_COMBINE_DIVIDE:
136: FNEvaluateFunctionMat_Private(ctx->f2,A,W,PETSC_FALSE);
137: FNEvaluateFunctionMat_Private(ctx->f1,A,B,PETSC_FALSE);
138: PetscMalloc1(ld,&ipiv);
139: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&n,Wa,&ld,ipiv,Ba,&ld,&info));
140: SlepcCheckLapackInfo("gesv",info);
141: PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n*n);
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: MatDenseRestoreArray(A,&Aa);
151: MatDenseRestoreArray(B,&Ba);
152: MatDenseRestoreArray(W,&Wa);
153: FN_FreeWorkMat(fn,&W);
154: return(0);
155: #endif
156: }
158: PetscErrorCode FNEvaluateFunctionMatVec_Combine(FN fn,Mat A,Vec v)
159: {
160: #if defined(PETSC_MISSING_LAPACK_GESV)
162: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESV - Lapack routines are unavailable");
163: #else
165: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
166: PetscScalar *va,*Za;
167: PetscBLASInt n,ld,*ipiv,info,one=1;
168: PetscInt m;
169: Mat Z;
170: Vec w;
173: MatGetSize(A,&m,NULL);
174: PetscBLASIntCast(m,&n);
175: ld = n;
177: switch (ctx->comb) {
178: case FN_COMBINE_ADD:
179: VecDuplicate(v,&w);
180: FNEvaluateFunctionMatVec(ctx->f1,A,w);
181: FNEvaluateFunctionMatVec(ctx->f2,A,v);
182: VecAXPY(v,1.0,w);
183: VecDestroy(&w);
184: break;
185: case FN_COMBINE_MULTIPLY:
186: VecDuplicate(v,&w);
187: FN_AllocateWorkMat(fn,A,&Z);
188: FNEvaluateFunctionMat_Private(ctx->f1,A,Z,PETSC_FALSE);
189: FNEvaluateFunctionMatVec_Private(ctx->f2,A,w,PETSC_FALSE);
190: MatMult(Z,w,v);
191: FN_FreeWorkMat(fn,&Z);
192: VecDestroy(&w);
193: break;
194: case FN_COMBINE_DIVIDE:
195: VecDuplicate(v,&w);
196: FN_AllocateWorkMat(fn,A,&Z);
197: FNEvaluateFunctionMat_Private(ctx->f2,A,Z,PETSC_FALSE);
198: FNEvaluateFunctionMatVec_Private(ctx->f1,A,v,PETSC_FALSE);
199: PetscMalloc1(ld,&ipiv);
200: MatDenseGetArray(Z,&Za);
201: VecGetArray(v,&va);
202: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&one,Za,&ld,ipiv,va,&ld,&info));
203: SlepcCheckLapackInfo("gesv",info);
204: PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n);
205: VecRestoreArray(v,&va);
206: MatDenseRestoreArray(Z,&Za);
207: PetscFree(ipiv);
208: FN_FreeWorkMat(fn,&Z);
209: VecDestroy(&w);
210: break;
211: case FN_COMBINE_COMPOSE:
212: FN_AllocateWorkMat(fn,A,&Z);
213: FNEvaluateFunctionMat_Private(ctx->f1,A,Z,PETSC_FALSE);
214: FNEvaluateFunctionMatVec_Private(ctx->f2,Z,v,PETSC_FALSE);
215: FN_FreeWorkMat(fn,&Z);
216: break;
217: }
218: return(0);
219: #endif
220: }
222: PetscErrorCode FNView_Combine(FN fn,PetscViewer viewer)
223: {
225: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
226: PetscBool isascii;
229: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
230: if (isascii) {
231: switch (ctx->comb) {
232: case FN_COMBINE_ADD:
233: PetscViewerASCIIPrintf(viewer," Two added functions f1+f2\n");
234: break;
235: case FN_COMBINE_MULTIPLY:
236: PetscViewerASCIIPrintf(viewer," Two multiplied functions f1*f2\n");
237: break;
238: case FN_COMBINE_DIVIDE:
239: PetscViewerASCIIPrintf(viewer," A quotient of two functions f1/f2\n");
240: break;
241: case FN_COMBINE_COMPOSE:
242: PetscViewerASCIIPrintf(viewer," Two composed functions f2(f1(.))\n");
243: break;
244: }
245: PetscViewerASCIIPushTab(viewer);
246: FNView(ctx->f1,viewer);
247: FNView(ctx->f2,viewer);
248: PetscViewerASCIIPopTab(viewer);
249: }
250: return(0);
251: }
253: static PetscErrorCode FNCombineSetChildren_Combine(FN fn,FNCombineType comb,FN f1,FN f2)
254: {
256: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
259: ctx->comb = comb;
260: PetscObjectReference((PetscObject)f1);
261: FNDestroy(&ctx->f1);
262: ctx->f1 = f1;
263: PetscLogObjectParent((PetscObject)fn,(PetscObject)ctx->f1);
264: PetscObjectReference((PetscObject)f2);
265: FNDestroy(&ctx->f2);
266: ctx->f2 = f2;
267: PetscLogObjectParent((PetscObject)fn,(PetscObject)ctx->f2);
268: return(0);
269: }
271: /*@
272: FNCombineSetChildren - Sets the two child functions that constitute this
273: combined function, and the way they must be combined.
275: Logically Collective on FN
277: Input Parameters:
278: + fn - the math function context
279: . comb - how to combine the functions (addition, multiplication, division or composition)
280: . f1 - first function
281: - f2 - second function
283: Level: intermediate
285: .seealso: FNCombineGetChildren()
286: @*/
287: PetscErrorCode FNCombineSetChildren(FN fn,FNCombineType comb,FN f1,FN f2)
288: {
296: PetscTryMethod(fn,"FNCombineSetChildren_C",(FN,FNCombineType,FN,FN),(fn,comb,f1,f2));
297: return(0);
298: }
300: static PetscErrorCode FNCombineGetChildren_Combine(FN fn,FNCombineType *comb,FN *f1,FN *f2)
301: {
303: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
306: if (comb) *comb = ctx->comb;
307: if (f1) {
308: if (!ctx->f1) {
309: FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f1);
310: PetscLogObjectParent((PetscObject)fn,(PetscObject)ctx->f1);
311: }
312: *f1 = ctx->f1;
313: }
314: if (f2) {
315: if (!ctx->f2) {
316: FNCreate(PetscObjectComm((PetscObject)fn),&ctx->f2);
317: PetscLogObjectParent((PetscObject)fn,(PetscObject)ctx->f2);
318: }
319: *f2 = ctx->f2;
320: }
321: return(0);
322: }
324: /*@
325: FNCombineGetChildren - Gets the two child functions that constitute this
326: combined function, and the way they are combined.
328: Not Collective
330: Input Parameter:
331: . fn - the math function context
333: Output Parameters:
334: + comb - how to combine the functions (addition, multiplication, division or composition)
335: . f1 - first function
336: - f2 - second function
338: Level: intermediate
340: .seealso: FNCombineSetChildren()
341: @*/
342: PetscErrorCode FNCombineGetChildren(FN fn,FNCombineType *comb,FN *f1,FN *f2)
343: {
348: PetscUseMethod(fn,"FNCombineGetChildren_C",(FN,FNCombineType*,FN*,FN*),(fn,comb,f1,f2));
349: return(0);
350: }
352: PetscErrorCode FNDuplicate_Combine(FN fn,MPI_Comm comm,FN *newfn)
353: {
355: FN_COMBINE *ctx = (FN_COMBINE*)fn->data,*ctx2 = (FN_COMBINE*)(*newfn)->data;
358: ctx2->comb = ctx->comb;
359: FNDuplicate(ctx->f1,comm,&ctx2->f1);
360: FNDuplicate(ctx->f2,comm,&ctx2->f2);
361: return(0);
362: }
364: PetscErrorCode FNDestroy_Combine(FN fn)
365: {
367: FN_COMBINE *ctx = (FN_COMBINE*)fn->data;
370: FNDestroy(&ctx->f1);
371: FNDestroy(&ctx->f2);
372: PetscFree(fn->data);
373: PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",NULL);
374: PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",NULL);
375: return(0);
376: }
378: SLEPC_EXTERN PetscErrorCode FNCreate_Combine(FN fn)
379: {
381: FN_COMBINE *ctx;
384: PetscNewLog(fn,&ctx);
385: fn->data = (void*)ctx;
387: fn->ops->evaluatefunction = FNEvaluateFunction_Combine;
388: fn->ops->evaluatederivative = FNEvaluateDerivative_Combine;
389: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Combine;
390: fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Combine;
391: fn->ops->view = FNView_Combine;
392: fn->ops->duplicate = FNDuplicate_Combine;
393: fn->ops->destroy = FNDestroy_Combine;
394: PetscObjectComposeFunction((PetscObject)fn,"FNCombineSetChildren_C",FNCombineSetChildren_Combine);
395: PetscObjectComposeFunction((PetscObject)fn,"FNCombineGetChildren_C",FNCombineGetChildren_Combine);
396: return(0);
397: }