Actual source code: fncombine.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }