Actual source code: fnbasic.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:    Basic FN routines
 12: */

 14: #include <slepc/private/fnimpl.h>
 15: #include <slepcblaslapack.h>

 17: PetscFunctionList FNList = 0;
 18: PetscBool         FNRegisterAllCalled = PETSC_FALSE;
 19: PetscClassId      FN_CLASSID = 0;
 20: PetscLogEvent     FN_Evaluate = 0;
 21: static PetscBool  FNPackageInitialized = PETSC_FALSE;

 23: const char *FNParallelTypes[] = {"REDUNDANT","SYNCHRONIZED","FNParallelType","FN_PARALLEL_",0};

 25: /*@C
 26:    FNFinalizePackage - This function destroys everything in the Slepc interface
 27:    to the FN package. It is called from SlepcFinalize().

 29:    Level: developer

 31: .seealso: SlepcFinalize()
 32: @*/
 33: PetscErrorCode FNFinalizePackage(void)
 34: {
 35:   PetscFunctionListDestroy(&FNList);
 36:   FNPackageInitialized = PETSC_FALSE;
 37:   FNRegisterAllCalled  = PETSC_FALSE;
 38:   PetscFunctionReturn(0);
 39: }

 41: /*@C
 42:   FNInitializePackage - This function initializes everything in the FN package.
 43:   It is called from PetscDLLibraryRegister() when using dynamic libraries, and
 44:   on the first call to FNCreate() when using static libraries.

 46:   Level: developer

 48: .seealso: SlepcInitialize()
 49: @*/
 50: PetscErrorCode FNInitializePackage(void)
 51: {
 52:   char           logList[256];
 53:   PetscBool      opt,pkg;
 54:   PetscClassId   classids[1];

 56:   if (FNPackageInitialized) PetscFunctionReturn(0);
 57:   FNPackageInitialized = PETSC_TRUE;
 58:   /* Register Classes */
 59:   PetscClassIdRegister("Math Function",&FN_CLASSID);
 60:   /* Register Constructors */
 61:   FNRegisterAll();
 62:   /* Register Events */
 63:   PetscLogEventRegister("FNEvaluate",FN_CLASSID,&FN_Evaluate);
 64:   /* Process Info */
 65:   classids[0] = FN_CLASSID;
 66:   PetscInfoProcessClass("fn",1,&classids[0]);
 67:   /* Process summary exclusions */
 68:   PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
 69:   if (opt) {
 70:     PetscStrInList("fn",logList,',',&pkg);
 71:     if (pkg) PetscLogEventDeactivateClass(FN_CLASSID);
 72:   }
 73:   /* Register package finalizer */
 74:   PetscRegisterFinalize(FNFinalizePackage);
 75:   PetscFunctionReturn(0);
 76: }

 78: /*@
 79:    FNCreate - Creates an FN context.

 81:    Collective

 83:    Input Parameter:
 84: .  comm - MPI communicator

 86:    Output Parameter:
 87: .  newfn - location to put the FN context

 89:    Level: beginner

 91: .seealso: FNDestroy(), FN
 92: @*/
 93: PetscErrorCode FNCreate(MPI_Comm comm,FN *newfn)
 94: {
 95:   FN             fn;

 98:   *newfn = 0;
 99:   FNInitializePackage();
100:   SlepcHeaderCreate(fn,FN_CLASSID,"FN","Math Function","FN",comm,FNDestroy,FNView);

102:   fn->alpha    = 1.0;
103:   fn->beta     = 1.0;
104:   fn->method   = 0;

106:   fn->nw       = 0;
107:   fn->cw       = 0;
108:   fn->data     = NULL;

110:   *newfn = fn;
111:   PetscFunctionReturn(0);
112: }

114: /*@C
115:    FNSetOptionsPrefix - Sets the prefix used for searching for all
116:    FN options in the database.

118:    Logically Collective on fn

120:    Input Parameters:
121: +  fn - the math function context
122: -  prefix - the prefix string to prepend to all FN option requests

124:    Notes:
125:    A hyphen (-) must NOT be given at the beginning of the prefix name.
126:    The first character of all runtime options is AUTOMATICALLY the
127:    hyphen.

129:    Level: advanced

131: .seealso: FNAppendOptionsPrefix()
132: @*/
133: PetscErrorCode FNSetOptionsPrefix(FN fn,const char *prefix)
134: {
136:   PetscObjectSetOptionsPrefix((PetscObject)fn,prefix);
137:   PetscFunctionReturn(0);
138: }

140: /*@C
141:    FNAppendOptionsPrefix - Appends to the prefix used for searching for all
142:    FN options in the database.

144:    Logically Collective on fn

146:    Input Parameters:
147: +  fn - the math function context
148: -  prefix - the prefix string to prepend to all FN option requests

150:    Notes:
151:    A hyphen (-) must NOT be given at the beginning of the prefix name.
152:    The first character of all runtime options is AUTOMATICALLY the hyphen.

154:    Level: advanced

156: .seealso: FNSetOptionsPrefix()
157: @*/
158: PetscErrorCode FNAppendOptionsPrefix(FN fn,const char *prefix)
159: {
161:   PetscObjectAppendOptionsPrefix((PetscObject)fn,prefix);
162:   PetscFunctionReturn(0);
163: }

165: /*@C
166:    FNGetOptionsPrefix - Gets the prefix used for searching for all
167:    FN options in the database.

169:    Not Collective

171:    Input Parameters:
172: .  fn - the math function context

174:    Output Parameters:
175: .  prefix - pointer to the prefix string used is returned

177:    Note:
178:    On the Fortran side, the user should pass in a string 'prefix' of
179:    sufficient length to hold the prefix.

181:    Level: advanced

183: .seealso: FNSetOptionsPrefix(), FNAppendOptionsPrefix()
184: @*/
185: PetscErrorCode FNGetOptionsPrefix(FN fn,const char *prefix[])
186: {
189:   PetscObjectGetOptionsPrefix((PetscObject)fn,prefix);
190:   PetscFunctionReturn(0);
191: }

193: /*@C
194:    FNSetType - Selects the type for the FN object.

196:    Logically Collective on fn

198:    Input Parameters:
199: +  fn   - the math function context
200: -  type - a known type

202:    Notes:
203:    The default is FNRATIONAL, which includes polynomials as a particular
204:    case as well as simple functions such as f(x)=x and f(x)=constant.

206:    Level: intermediate

208: .seealso: FNGetType()
209: @*/
210: PetscErrorCode FNSetType(FN fn,FNType type)
211: {
212:   PetscErrorCode (*r)(FN);
213:   PetscBool      match;


218:   PetscObjectTypeCompare((PetscObject)fn,type,&match);
219:   if (match) PetscFunctionReturn(0);

221:   PetscFunctionListFind(FNList,type,&r);

224:   if (fn->ops->destroy) (*fn->ops->destroy)(fn);
225:   PetscMemzero(fn->ops,sizeof(struct _FNOps));

227:   PetscObjectChangeTypeName((PetscObject)fn,type);
228:   (*r)(fn);
229:   PetscFunctionReturn(0);
230: }

232: /*@C
233:    FNGetType - Gets the FN type name (as a string) from the FN context.

235:    Not Collective

237:    Input Parameter:
238: .  fn - the math function context

240:    Output Parameter:
241: .  type - name of the math function

243:    Level: intermediate

245: .seealso: FNSetType()
246: @*/
247: PetscErrorCode FNGetType(FN fn,FNType *type)
248: {
251:   *type = ((PetscObject)fn)->type_name;
252:   PetscFunctionReturn(0);
253: }

255: /*@
256:    FNSetScale - Sets the scaling parameters that define the matematical function.

258:    Logically Collective on fn

260:    Input Parameters:
261: +  fn    - the math function context
262: .  alpha - inner scaling (argument)
263: -  beta  - outer scaling (result)

265:    Notes:
266:    Given a function f(x) specified by the FN type, the scaling parameters can
267:    be used to realize the function beta*f(alpha*x). So when these values are given,
268:    the procedure for function evaluation will first multiply the argument by alpha,
269:    then evaluate the function itself, and finally scale the result by beta.
270:    Likewise, these values are also considered when evaluating the derivative.

272:    If you want to provide only one of the two scaling factors, set the other
273:    one to 1.0.

275:    Level: intermediate

277: .seealso: FNGetScale(), FNEvaluateFunction()
278: @*/
279: PetscErrorCode FNSetScale(FN fn,PetscScalar alpha,PetscScalar beta)
280: {
285:   fn->alpha = alpha;
286:   fn->beta  = beta;
287:   PetscFunctionReturn(0);
288: }

290: /*@
291:    FNGetScale - Gets the scaling parameters that define the matematical function.

293:    Not Collective

295:    Input Parameter:
296: .  fn    - the math function context

298:    Output Parameters:
299: +  alpha - inner scaling (argument)
300: -  beta  - outer scaling (result)

302:    Level: intermediate

304: .seealso: FNSetScale()
305: @*/
306: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)
307: {
309:   if (alpha) *alpha = fn->alpha;
310:   if (beta)  *beta  = fn->beta;
311:   PetscFunctionReturn(0);
312: }

314: /*@
315:    FNSetMethod - Selects the method to be used to evaluate functions of matrices.

317:    Logically Collective on fn

319:    Input Parameters:
320: +  fn   - the math function context
321: -  meth - an index identifying the method

323:    Options Database Key:
324: .  -fn_method <meth> - Sets the method

326:    Notes:
327:    In some FN types there are more than one algorithms available for computing
328:    matrix functions. In that case, this function allows choosing the wanted method.

330:    If meth is currently set to 0 (the default) and the input argument A of
331:    FNEvaluateFunctionMat() is a symmetric/Hermitian matrix, then the computation
332:    is done via the eigendecomposition of A, rather than with the general algorithm.

334:    Level: intermediate

336: .seealso: FNGetMethod(), FNEvaluateFunctionMat()
337: @*/
338: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)
339: {
344:   fn->method = meth;
345:   PetscFunctionReturn(0);
346: }

348: /*@
349:    FNGetMethod - Gets the method currently used in the FN.

351:    Not Collective

353:    Input Parameter:
354: .  fn - the math function context

356:    Output Parameter:
357: .  meth - identifier of the method

359:    Level: intermediate

361: .seealso: FNSetMethod()
362: @*/
363: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)
364: {
367:   *meth = fn->method;
368:   PetscFunctionReturn(0);
369: }

371: /*@
372:    FNSetParallel - Selects the mode of operation in parallel runs.

374:    Logically Collective on fn

376:    Input Parameters:
377: +  fn    - the math function context
378: -  pmode - the parallel mode

380:    Options Database Key:
381: .  -fn_parallel <mode> - Sets the parallel mode, either 'redundant' or 'synchronized'

383:    Notes:
384:    This is relevant only when the function is evaluated on a matrix, with
385:    either FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().

387:    In the 'redundant' parallel mode, all processes will make the computation
388:    redundantly, starting from the same data, and producing the same result.
389:    This result may be slightly different in the different processes if using a
390:    multithreaded BLAS library, which may cause issues in ill-conditioned problems.

392:    In the 'synchronized' parallel mode, only the first MPI process performs the
393:    computation and then the computed matrix is broadcast to the other
394:    processes in the communicator. This communication is done automatically at
395:    the end of FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().

397:    Level: advanced

399: .seealso: FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec(), FNGetParallel()
400: @*/
401: PetscErrorCode FNSetParallel(FN fn,FNParallelType pmode)
402: {
405:   fn->pmode = pmode;
406:   PetscFunctionReturn(0);
407: }

409: /*@
410:    FNGetParallel - Gets the mode of operation in parallel runs.

412:    Not Collective

414:    Input Parameter:
415: .  fn - the math function context

417:    Output Parameter:
418: .  pmode - the parallel mode

420:    Level: advanced

422: .seealso: FNSetParallel()
423: @*/
424: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)
425: {
428:   *pmode = fn->pmode;
429:   PetscFunctionReturn(0);
430: }

432: /*@
433:    FNEvaluateFunction - Computes the value of the function f(x) for a given x.

435:    Not collective

437:    Input Parameters:
438: +  fn - the math function context
439: -  x  - the value where the function must be evaluated

441:    Output Parameter:
442: .  y  - the result of f(x)

444:    Note:
445:    Scaling factors are taken into account, so the actual function evaluation
446:    will return beta*f(alpha*x).

448:    Level: intermediate

450: .seealso: FNEvaluateDerivative(), FNEvaluateFunctionMat(), FNSetScale()
451: @*/
452: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)
453: {
454:   PetscScalar    xf,yf;

459:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
460:   xf = fn->alpha*x;
461:   (*fn->ops->evaluatefunction)(fn,xf,&yf);
462:   *y = fn->beta*yf;
463:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
464:   PetscFunctionReturn(0);
465: }

467: /*@
468:    FNEvaluateDerivative - Computes the value of the derivative f'(x) for a given x.

470:    Not Collective

472:    Input Parameters:
473: +  fn - the math function context
474: -  x  - the value where the derivative must be evaluated

476:    Output Parameter:
477: .  y  - the result of f'(x)

479:    Note:
480:    Scaling factors are taken into account, so the actual derivative evaluation will
481:    return alpha*beta*f'(alpha*x).

483:    Level: intermediate

485: .seealso: FNEvaluateFunction()
486: @*/
487: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)
488: {
489:   PetscScalar    xf,yf;

494:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
495:   xf = fn->alpha*x;
496:   (*fn->ops->evaluatederivative)(fn,xf,&yf);
497:   *y = fn->alpha*fn->beta*yf;
498:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
499:   PetscFunctionReturn(0);
500: }

502: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,const PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)
503: {
504:   PetscInt       i,j;
505:   PetscBLASInt   n,k,ld,lwork,info;
506:   PetscScalar    *Q,*W,*work,adummy,a,x,y,one=1.0,zero=0.0;
507:   PetscReal      *eig,dummy;
508: #if defined(PETSC_USE_COMPLEX)
509:   PetscReal      *rwork,rdummy;
510: #endif

512:   PetscBLASIntCast(m,&n);
513:   ld = n;
514:   k = firstonly? 1: n;

516:   /* workspace query and memory allocation */
517:   lwork = -1;
518: #if defined(PETSC_USE_COMPLEX)
519:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&rdummy,&info));
520:   PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
521:   PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork);
522: #else
523:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&info));
524:   PetscBLASIntCast((PetscInt)a,&lwork);
525:   PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work);
526: #endif

528:   /* compute eigendecomposition */
529:   for (j=0;j<n;j++) for (i=j;i<n;i++) Q[i+j*ld] = As[i+j*ld];
530: #if defined(PETSC_USE_COMPLEX)
531:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
532: #else
533:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
534: #endif
535:   SlepcCheckLapackInfo("syev",info);

537:   /* W = f(Lambda)*Q' */
538:   for (i=0;i<n;i++) {
539:     x = fn->alpha*eig[i];
540:     (*fn->ops->evaluatefunction)(fn,x,&y);  /* y = f(x) */
541:     for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
542:   }
543:   /* Bs = Q*W */
544:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
545: #if defined(PETSC_USE_COMPLEX)
546:   PetscFree5(eig,Q,W,work,rwork);
547: #else
548:   PetscFree4(eig,Q,W,work);
549: #endif
550:   PetscLogFlops(9.0*n*n*n+2.0*n*n*n);
551:   PetscFunctionReturn(0);
552: }

554: /*
555:    FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
556:    compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
557:    decomposition of A is A=Q*D*Q'
558: */
559: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)
560: {
561:   PetscInt          m;
562:   const PetscScalar *As;
563:   PetscScalar       *Bs;

565:   MatDenseGetArrayRead(A,&As);
566:   MatDenseGetArray(B,&Bs);
567:   MatGetSize(A,&m,NULL);
568:   FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE);
569:   MatDenseRestoreArrayRead(A,&As);
570:   MatDenseRestoreArray(B,&Bs);
571:   PetscFunctionReturn(0);
572: }

574: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
575: {
576:   PetscBool      set,flg,symm=PETSC_FALSE;
577:   PetscInt       m,n;
578:   PetscMPIInt    size,rank;
579:   PetscScalar    *pF;
580:   Mat            M,F;

582:   /* destination matrix */
583:   F = B?B:A;

585:   /* check symmetry of A */
586:   MatIsHermitianKnown(A,&set,&flg);
587:   symm = set? flg: PETSC_FALSE;

589:   MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
590:   MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
591:   if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {

593:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
594:     if (symm && !fn->method) {  /* prefer diagonalization */
595:       PetscInfo(fn,"Computing matrix function via diagonalization\n");
596:       FNEvaluateFunctionMat_Sym_Default(fn,A,F);
597:     } else {
598:       /* scale argument */
599:       if (fn->alpha!=(PetscScalar)1.0) {
600:         FN_AllocateWorkMat(fn,A,&M);
601:         MatScale(M,fn->alpha);
602:       } else M = A;
603:       if (fn->ops->evaluatefunctionmat[fn->method]) (*fn->ops->evaluatefunctionmat[fn->method])(fn,M,F);
604:       else {
607:       }
608:       if (fn->alpha!=(PetscScalar)1.0) FN_FreeWorkMat(fn,&M);
609:       /* scale result */
610:       MatScale(F,fn->beta);
611:     }
612:     PetscFPTrapPop();
613:   }
614:   if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {  /* synchronize */
615:     MatGetSize(A,&m,&n);
616:     MatDenseGetArray(F,&pF);
617:     MPI_Bcast(pF,n*n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
618:     MatDenseRestoreArray(F,&pF);
619:   }
620:   PetscFunctionReturn(0);
621: }

623: /*@
624:    FNEvaluateFunctionMat - Computes the value of the function f(A) for a given
625:    matrix A, where the result is also a matrix.

627:    Logically Collective on fn

629:    Input Parameters:
630: +  fn - the math function context
631: -  A  - matrix on which the function must be evaluated

633:    Output Parameter:
634: .  B  - (optional) matrix resulting from evaluating f(A)

636:    Notes:
637:    Matrix A must be a square sequential dense Mat, with all entries equal on
638:    all processes (otherwise each process will compute different results).
639:    If matrix B is provided, it must also be a square sequential dense Mat, and
640:    both matrices must have the same dimensions. If B is NULL (or B=A) then the
641:    function will perform an in-place computation, overwriting A with f(A).

643:    If A is known to be real symmetric or complex Hermitian then it is
644:    recommended to set the appropriate flag with MatSetOption(), because
645:    symmetry can sometimes be exploited by the algorithm.

647:    Scaling factors are taken into account, so the actual function evaluation
648:    will return beta*f(alpha*A).

650:    Level: advanced

652: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
653: @*/
654: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
655: {
656:   PetscBool      inplace=PETSC_FALSE;
657:   PetscInt       m,n,n1;

663:   if (B) {
666:   } else inplace = PETSC_TRUE;
668:   MatGetSize(A,&m,&n);
670:   if (!inplace) {
672:     n1 = n;
673:     MatGetSize(B,&m,&n);
676:   }

678:   /* evaluate matrix function */
679:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
680:   FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE);
681:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
682:   PetscFunctionReturn(0);
683: }

685: /*
686:    FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
687:    and then copies the first column.
688: */
689: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
690: {
691:   Mat            F;

693:   FN_AllocateWorkMat(fn,A,&F);
694:   if (fn->ops->evaluatefunctionmat[fn->method]) (*fn->ops->evaluatefunctionmat[fn->method])(fn,A,F);
695:   else {
698:   }
699:   MatGetColumnVector(F,v,0);
700:   FN_FreeWorkMat(fn,&F);
701:   PetscFunctionReturn(0);
702: }

704: /*
705:    FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
706:    compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
707:    decomposition of A is A=Q*D*Q'. Only the first column is computed.
708: */
709: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
710: {
711:   PetscInt          m;
712:   const PetscScalar *As;
713:   PetscScalar       *vs;

715:   MatDenseGetArrayRead(A,&As);
716:   VecGetArray(v,&vs);
717:   MatGetSize(A,&m,NULL);
718:   FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE);
719:   MatDenseRestoreArrayRead(A,&As);
720:   VecRestoreArray(v,&vs);
721:   PetscFunctionReturn(0);
722: }

724: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
725: {
726:   PetscBool      set,flg,symm=PETSC_FALSE;
727:   PetscInt       m,n;
728:   Mat            M;
729:   PetscMPIInt    size,rank;
730:   PetscScalar    *pv;

732:   /* check symmetry of A */
733:   MatIsHermitianKnown(A,&set,&flg);
734:   symm = set? flg: PETSC_FALSE;

736:   /* evaluate matrix function */
737:   MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
738:   MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
739:   if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
740:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
741:     if (symm && !fn->method) {  /* prefer diagonalization */
742:       PetscInfo(fn,"Computing matrix function via diagonalization\n");
743:       FNEvaluateFunctionMatVec_Sym_Default(fn,A,v);
744:     } else {
745:       /* scale argument */
746:       if (fn->alpha!=(PetscScalar)1.0) {
747:         FN_AllocateWorkMat(fn,A,&M);
748:         MatScale(M,fn->alpha);
749:       } else M = A;
750:       if (fn->ops->evaluatefunctionmatvec[fn->method]) (*fn->ops->evaluatefunctionmatvec[fn->method])(fn,M,v);
751:       else FNEvaluateFunctionMatVec_Default(fn,M,v);
752:       if (fn->alpha!=(PetscScalar)1.0) FN_FreeWorkMat(fn,&M);
753:       /* scale result */
754:       VecScale(v,fn->beta);
755:     }
756:     PetscFPTrapPop();
757:   }

759:   /* synchronize */
760:   if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
761:     MatGetSize(A,&m,&n);
762:     VecGetArray(v,&pv);
763:     MPI_Bcast(pv,n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
764:     VecRestoreArray(v,&pv);
765:   }
766:   PetscFunctionReturn(0);
767: }

769: /*@
770:    FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
771:    for a given matrix A.

773:    Logically Collective on fn

775:    Input Parameters:
776: +  fn - the math function context
777: -  A  - matrix on which the function must be evaluated

779:    Output Parameter:
780: .  v  - vector to hold the first column of f(A)

782:    Notes:
783:    This operation is similar to FNEvaluateFunctionMat() but returns only
784:    the first column of f(A), hence saving computations in most cases.

786:    Level: advanced

788: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
789: @*/
790: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
791: {
792:   PetscInt       m,n;

801:   MatGetSize(A,&m,&n);
803:   VecGetSize(v,&m);
805:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
806:   FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE);
807:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
808:   PetscFunctionReturn(0);
809: }

811: /*@
812:    FNSetFromOptions - Sets FN options from the options database.

814:    Collective on fn

816:    Input Parameters:
817: .  fn - the math function context

819:    Notes:
820:    To see all options, run your program with the -help option.

822:    Level: beginner

824: .seealso: FNSetOptionsPrefix()
825: @*/
826: PetscErrorCode FNSetFromOptions(FN fn)
827: {
829:   char           type[256];
830:   PetscScalar    array[2];
831:   PetscInt       k,meth;
832:   PetscBool      flg;
833:   FNParallelType pmode;

836:   FNRegisterAll();
837:   ierr = PetscObjectOptionsBegin((PetscObject)fn);
838:     PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg);
839:     if (flg) FNSetType(fn,type);
840:     else if (!((PetscObject)fn)->type_name) FNSetType(fn,FNRATIONAL);

842:     k = 2;
843:     array[0] = 0.0; array[1] = 0.0;
844:     PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg);
845:     if (flg) {
846:       if (k<2) array[1] = 1.0;
847:       FNSetScale(fn,array[0],array[1]);
848:     }

850:     PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg);
851:     if (flg) FNSetMethod(fn,meth);

853:     PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg);
854:     if (flg) FNSetParallel(fn,pmode);

856:     if (fn->ops->setfromoptions) (*fn->ops->setfromoptions)(PetscOptionsObject,fn);
857:     PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)fn);
858:   ierr = PetscOptionsEnd();
859:   PetscFunctionReturn(0);
860: }

862: /*@C
863:    FNView - Prints the FN data structure.

865:    Collective on fn

867:    Input Parameters:
868: +  fn - the math function context
869: -  viewer - optional visualization context

871:    Note:
872:    The available visualization contexts include
873: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
874: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
875:          output where only the first processor opens
876:          the file.  All other processors send their
877:          data to the first processor to print.

879:    The user can open an alternative visualization context with
880:    PetscViewerASCIIOpen() - output to a specified file.

882:    Level: beginner

884: .seealso: FNCreate()
885: @*/
886: PetscErrorCode FNView(FN fn,PetscViewer viewer)
887: {
888:   PetscBool      isascii;
889:   PetscMPIInt    size;

892:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer);
895:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
896:   if (isascii) {
897:     PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer);
898:     MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
899:     if (size>1) PetscViewerASCIIPrintf(viewer,"  parallel operation mode: %s\n",FNParallelTypes[fn->pmode]);
900:     if (fn->ops->view) {
901:       PetscViewerASCIIPushTab(viewer);
902:       (*fn->ops->view)(fn,viewer);
903:       PetscViewerASCIIPopTab(viewer);
904:     }
905:   }
906:   PetscFunctionReturn(0);
907: }

909: /*@C
910:    FNViewFromOptions - View from options

912:    Collective on FN

914:    Input Parameters:
915: +  fn   - the math function context
916: .  obj  - optional object
917: -  name - command line option

919:    Level: intermediate

921: .seealso: FNView(), FNCreate()
922: @*/
923: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])
924: {
926:   PetscObjectViewFromOptions((PetscObject)fn,obj,name);
927:   PetscFunctionReturn(0);
928: }

930: /*@
931:    FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
932:    different communicator.

934:    Collective on fn

936:    Input Parameters:
937: +  fn   - the math function context
938: -  comm - MPI communicator

940:    Output Parameter:
941: .  newfn - location to put the new FN context

943:    Note:
944:    In order to use the same MPI communicator as in the original object,
945:    use PetscObjectComm((PetscObject)fn).

947:    Level: developer

949: .seealso: FNCreate()
950: @*/
951: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)
952: {
953:   FNType         type;
954:   PetscScalar    alpha,beta;
955:   PetscInt       meth;
956:   FNParallelType ptype;

961:   FNCreate(comm,newfn);
962:   FNGetType(fn,&type);
963:   FNSetType(*newfn,type);
964:   FNGetScale(fn,&alpha,&beta);
965:   FNSetScale(*newfn,alpha,beta);
966:   FNGetMethod(fn,&meth);
967:   FNSetMethod(*newfn,meth);
968:   FNGetParallel(fn,&ptype);
969:   FNSetParallel(*newfn,ptype);
970:   if (fn->ops->duplicate) (*fn->ops->duplicate)(fn,comm,newfn);
971:   PetscFunctionReturn(0);
972: }

974: /*@C
975:    FNDestroy - Destroys FN context that was created with FNCreate().

977:    Collective on fn

979:    Input Parameter:
980: .  fn - the math function context

982:    Level: beginner

984: .seealso: FNCreate()
985: @*/
986: PetscErrorCode FNDestroy(FN *fn)
987: {
988:   PetscInt       i;

990:   if (!*fn) PetscFunctionReturn(0);
992:   if (--((PetscObject)(*fn))->refct > 0) { *fn = 0; PetscFunctionReturn(0); }
993:   if ((*fn)->ops->destroy) (*(*fn)->ops->destroy)(*fn);
994:   for (i=0;i<(*fn)->nw;i++) MatDestroy(&(*fn)->W[i]);
995:   PetscHeaderDestroy(fn);
996:   PetscFunctionReturn(0);
997: }

999: /*@C
1000:    FNRegister - Adds a mathematical function to the FN package.

1002:    Not collective

1004:    Input Parameters:
1005: +  name - name of a new user-defined FN
1006: -  function - routine to create context

1008:    Notes:
1009:    FNRegister() may be called multiple times to add several user-defined functions.

1011:    Level: advanced

1013: .seealso: FNRegisterAll()
1014: @*/
1015: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1016: {
1017:   FNInitializePackage();
1018:   PetscFunctionListAdd(&FNList,name,function);
1019:   PetscFunctionReturn(0);
1020: }