Actual source code: fncombine.c

slepc-3.17.2 2022-08-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>
 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: }