Actual source code: fnbasic.c

slepc-3.10.1 2018-10-23
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2018, 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>      /*I "slepcfn.h" I*/
 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: {

 38:   PetscFunctionListDestroy(&FNList);
 39:   FNPackageInitialized = PETSC_FALSE;
 40:   FNRegisterAllCalled  = PETSC_FALSE;
 41:   return(0);
 42: }

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

 49:   Level: developer

 51: .seealso: SlepcInitialize()
 52: @*/
 53: PetscErrorCode FNInitializePackage(void)
 54: {
 55:   char           logList[256];
 56:   PetscBool      opt,pkg;

 60:   if (FNPackageInitialized) return(0);
 61:   FNPackageInitialized = PETSC_TRUE;
 62:   /* Register Classes */
 63:   PetscClassIdRegister("Math Function",&FN_CLASSID);
 64:   /* Register Constructors */
 65:   FNRegisterAll();
 66:   /* Register Events */
 67:   PetscLogEventRegister("FNEvaluate",FN_CLASSID,&FN_Evaluate);
 68:   /* Process info exclusions */
 69:   PetscOptionsGetString(NULL,NULL,"-info_exclude",logList,sizeof(logList),&opt);
 70:   if (opt) {
 71:     PetscStrInList("fn",logList,',',&pkg);
 72:     if (pkg) { PetscInfoDeactivateClass(FN_CLASSID); }
 73:   }
 74:   /* Process summary exclusions */
 75:   PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
 76:   if (opt) {
 77:     PetscStrInList("fn",logList,',',&pkg);
 78:     if (pkg) { PetscLogEventDeactivateClass(FN_CLASSID); }
 79:   }
 80:   /* Register package finalizer */
 81:   PetscRegisterFinalize(FNFinalizePackage);
 82:   return(0);
 83: }

 85: /*@
 86:    FNCreate - Creates an FN context.

 88:    Collective on MPI_Comm

 90:    Input Parameter:
 91: .  comm - MPI communicator

 93:    Output Parameter:
 94: .  newfn - location to put the FN context

 96:    Level: beginner

 98: .seealso: FNDestroy(), FN
 99: @*/
100: PetscErrorCode FNCreate(MPI_Comm comm,FN *newfn)
101: {
102:   FN             fn;

107:   *newfn = 0;
108:   FNInitializePackage();
109:   SlepcHeaderCreate(fn,FN_CLASSID,"FN","Math Function","FN",comm,FNDestroy,FNView);

111:   fn->alpha    = 1.0;
112:   fn->beta     = 1.0;
113:   fn->method   = 0;

115:   fn->nw       = 0;
116:   fn->cw       = 0;
117:   fn->data     = NULL;

119:   *newfn = fn;
120:   return(0);
121: }

123: /*@C
124:    FNSetOptionsPrefix - Sets the prefix used for searching for all
125:    FN options in the database.

127:    Logically Collective on FN

129:    Input Parameters:
130: +  fn - the math function context
131: -  prefix - the prefix string to prepend to all FN option requests

133:    Notes:
134:    A hyphen (-) must NOT be given at the beginning of the prefix name.
135:    The first character of all runtime options is AUTOMATICALLY the
136:    hyphen.

138:    Level: advanced

140: .seealso: FNAppendOptionsPrefix()
141: @*/
142: PetscErrorCode FNSetOptionsPrefix(FN fn,const char *prefix)
143: {

148:   PetscObjectSetOptionsPrefix((PetscObject)fn,prefix);
149:   return(0);
150: }

152: /*@C
153:    FNAppendOptionsPrefix - Appends to the prefix used for searching for all
154:    FN options in the database.

156:    Logically Collective on FN

158:    Input Parameters:
159: +  fn - the math function context
160: -  prefix - the prefix string to prepend to all FN option requests

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

166:    Level: advanced

168: .seealso: FNSetOptionsPrefix()
169: @*/
170: PetscErrorCode FNAppendOptionsPrefix(FN fn,const char *prefix)
171: {

176:   PetscObjectAppendOptionsPrefix((PetscObject)fn,prefix);
177:   return(0);
178: }

180: /*@C
181:    FNGetOptionsPrefix - Gets the prefix used for searching for all
182:    FN options in the database.

184:    Not Collective

186:    Input Parameters:
187: .  fn - the math function context

189:    Output Parameters:
190: .  prefix - pointer to the prefix string used is returned

192:    Note:
193:    On the Fortran side, the user should pass in a string 'prefix' of
194:    sufficient length to hold the prefix.

196:    Level: advanced

198: .seealso: FNSetOptionsPrefix(), FNAppendOptionsPrefix()
199: @*/
200: PetscErrorCode FNGetOptionsPrefix(FN fn,const char *prefix[])
201: {

207:   PetscObjectGetOptionsPrefix((PetscObject)fn,prefix);
208:   return(0);
209: }

211: /*@C
212:    FNSetType - Selects the type for the FN object.

214:    Logically Collective on FN

216:    Input Parameter:
217: +  fn   - the math function context
218: -  type - a known type

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

224:    Level: intermediate

226: .seealso: FNGetType()
227: @*/
228: PetscErrorCode FNSetType(FN fn,FNType type)
229: {
230:   PetscErrorCode ierr,(*r)(FN);
231:   PetscBool      match;


237:   PetscObjectTypeCompare((PetscObject)fn,type,&match);
238:   if (match) return(0);

240:    PetscFunctionListFind(FNList,type,&r);
241:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested FN type %s",type);

243:   if (fn->ops->destroy) { (*fn->ops->destroy)(fn); }
244:   PetscMemzero(fn->ops,sizeof(struct _FNOps));

246:   PetscObjectChangeTypeName((PetscObject)fn,type);
247:   (*r)(fn);
248:   return(0);
249: }

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

254:    Not Collective

256:    Input Parameter:
257: .  fn - the math function context

259:    Output Parameter:
260: .  name - name of the math function

262:    Level: intermediate

264: .seealso: FNSetType()
265: @*/
266: PetscErrorCode FNGetType(FN fn,FNType *type)
267: {
271:   *type = ((PetscObject)fn)->type_name;
272:   return(0);
273: }

275: /*@
276:    FNSetScale - Sets the scaling parameters that define the matematical function.

278:    Logically Collective on FN

280:    Input Parameters:
281: +  fn    - the math function context
282: .  alpha - inner scaling (argument)
283: -  beta  - outer scaling (result)

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

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

295:    Level: intermediate

297: .seealso: FNGetScale(), FNEvaluateFunction()
298: @*/
299: PetscErrorCode FNSetScale(FN fn,PetscScalar alpha,PetscScalar beta)
300: {
305:   fn->alpha = alpha;
306:   fn->beta  = beta;
307:   return(0);
308: }

310: /*@
311:    FNGetScale - Gets the scaling parameters that define the matematical function.

313:    Not Collective

315:    Input Parameter:
316: .  fn    - the math function context

318:    Output Parameters:
319: +  alpha - inner scaling (argument)
320: -  beta  - outer scaling (result)

322:    Level: intermediate

324: .seealso: FNSetScale()
325: @*/
326: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)
327: {
330:   if (alpha) *alpha = fn->alpha;
331:   if (beta)  *beta  = fn->beta;
332:   return(0);
333: }

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

338:    Logically Collective on FN

340:    Input Parameter:
341: +  fn   - the math function context
342: -  meth - an index indentifying the method

344:    Options Database Key:
345: .  -fn_method <meth> - Sets the method

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

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

355:    Level: intermediate

357: .seealso: FNGetMethod(), FNEvaluateFunctionMat()
358: @*/
359: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)
360: {
364:   if (meth<0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
365:   if (meth>FN_MAX_SOLVE) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
366:   fn->method = meth;
367:   return(0);
368: }

370: /*@
371:    FNGetMethod - Gets the method currently used in the FN.

373:    Not Collective

375:    Input Parameter:
376: .  fn - the math function context

378:    Output Parameter:
379: .  meth - identifier of the method

381:    Level: intermediate

383: .seealso: FNSetMethod()
384: @*/
385: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)
386: {
390:   *meth = fn->method;
391:   return(0);
392: }

394: /*@
395:    FNSetParallel - Selects the mode of operation in parallel runs.

397:    Logically Collective on FN

399:    Input Parameter:
400: +  fn    - the math function context
401: -  pmode - the parallel mode

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

406:    Notes:
407:    This is relevant only when the function is evaluated on a matrix, with
408:    either FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().

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

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

420:    Level: advanced

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

433: /*@
434:    FNGetParallel - Gets the mode of operation in parallel runs.

436:    Not Collective

438:    Input Parameter:
439: .  fn - the math function context

441:    Output Parameter:
442: .  pmode - the parallel mode

444:    Level: advanced

446: .seealso: FNSetParallel()
447: @*/
448: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)
449: {
453:   *pmode = fn->pmode;
454:   return(0);
455: }

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

460:    Not collective

462:    Input Parameters:
463: +  fn - the math function context
464: -  x  - the value where the function must be evaluated

466:    Output Parameter:
467: .  y  - the result of f(x)

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

473:    Level: intermediate

475: .seealso: FNEvaluateDerivative(), FNEvaluateFunctionMat(), FNSetScale()
476: @*/
477: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)
478: {
480:   PetscScalar    xf,yf;

486:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
487:   xf = fn->alpha*x;
488:   (*fn->ops->evaluatefunction)(fn,xf,&yf);
489:   *y = fn->beta*yf;
490:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
491:   return(0);
492: }

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

497:    Not Collective

499:    Input Parameters:
500: +  fn - the math function context
501: -  x  - the value where the derivative must be evaluated

503:    Output Parameter:
504: .  y  - the result of f'(x)

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

510:    Level: intermediate

512: .seealso: FNEvaluateFunction()
513: @*/
514: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)
515: {
517:   PetscScalar    xf,yf;

523:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
524:   xf = fn->alpha*x;
525:   (*fn->ops->evaluatederivative)(fn,xf,&yf);
526:   *y = fn->alpha*fn->beta*yf;
527:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
528:   return(0);
529: }

531: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)
532: {
533: #if defined(PETSC_MISSING_LAPACK_SYEV) || defined(SLEPC_MISSING_LAPACK_LACPY)
535:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SYEV/LACPY - Lapack routines are unavailable");
536: #else
538:   PetscInt       i,j;
539:   PetscBLASInt   n,k,ld,lwork,info;
540:   PetscScalar    *Q,*W,*work,a,x,y,one=1.0,zero=0.0;
541:   PetscReal      *eig,dummy;
542: #if defined(PETSC_USE_COMPLEX)
543:   PetscReal      *rwork,rdummy;
544: #endif

547:   PetscBLASIntCast(m,&n);
548:   ld = n;
549:   k = firstonly? 1: n;

551:   /* workspace query and memory allocation */
552:   lwork = -1;
553: #if defined(PETSC_USE_COMPLEX)
554:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,As,&ld,&dummy,&a,&lwork,&rdummy,&info));
555:   PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
556:   PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork);
557: #else
558:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,As,&ld,&dummy,&a,&lwork,&info));
559:   PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
560:   PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work);
561: #endif

563:   /* compute eigendecomposition */
564:   PetscStackCallBLAS("LAPACKlacpy",LAPACKlacpy_("L",&n,&n,As,&ld,Q,&ld));
565: #if defined(PETSC_USE_COMPLEX)
566:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
567: #else
568:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
569: #endif
570:   SlepcCheckLapackInfo("syev",info);

572:   /* W = f(Lambda)*Q' */
573:   for (i=0;i<n;i++) {
574:     x = fn->alpha*eig[i];
575:     (*fn->ops->evaluatefunction)(fn,x,&y);  /* y = f(x) */
576:     for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
577:   }
578:   /* Bs = Q*W */
579:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
580: #if defined(PETSC_USE_COMPLEX)
581:   PetscFree5(eig,Q,W,work,rwork);
582: #else
583:   PetscFree4(eig,Q,W,work);
584: #endif
585:   PetscLogFlops(9.0*n*n*n+2.0*n*n*n);
586:   return(0);
587: #endif
588: }

590: /*
591:    FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
592:    compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
593:    decomposition of A is A=Q*D*Q'
594: */
595: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)
596: {
598:   PetscInt       m;
599:   PetscScalar    *As,*Bs;

602:   MatDenseGetArray(A,&As);
603:   MatDenseGetArray(B,&Bs);
604:   MatGetSize(A,&m,NULL);
605:   FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE);
606:   MatDenseRestoreArray(A,&As);
607:   MatDenseRestoreArray(B,&Bs);
608:   return(0);
609: }

611: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
612: {
614:   PetscBool      set,flg,symm=PETSC_FALSE;
615:   PetscInt       m,n;
616:   PetscMPIInt    size,rank;
617:   PetscScalar    *pF;
618:   Mat            M,F;

621:   /* destination matrix */
622:   F = B?B:A;

624:   /* check symmetry of A */
625:   MatIsHermitianKnown(A,&set,&flg);
626:   symm = set? flg: PETSC_FALSE;

628:   MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
629:   MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
630:   if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
631:     
632:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
633:     if (symm && !fn->method) {  /* prefer diagonalization */
634:       PetscInfo(fn,"Computing matrix function via diagonalization\n");
635:       FNEvaluateFunctionMat_Sym_Default(fn,A,F);
636:     } else {
637:       /* scale argument */
638:       if (fn->alpha!=(PetscScalar)1.0) {
639:         FN_AllocateWorkMat(fn,A,&M);
640:         MatScale(M,fn->alpha);
641:       } else M = A;
642:       if (fn->ops->evaluatefunctionmat[fn->method]) {
643:         (*fn->ops->evaluatefunctionmat[fn->method])(fn,M,F);
644:       } else if (!fn->method) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
645:       else SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
646:       if (fn->alpha!=(PetscScalar)1.0) {
647:         FN_FreeWorkMat(fn,&M);
648:       }
649:       /* scale result */
650:       MatScale(F,fn->beta);
651:     }
652:     PetscFPTrapPop();
653:   }
654:   if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {  /* synchronize */
655:     MatGetSize(A,&m,&n);
656:     MatDenseGetArray(F,&pF);
657:     MPI_Bcast(pF,n*n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
658:     MatDenseRestoreArray(F,&pF);
659:   }
660:   return(0);
661: }

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

667:    Logically Collective on FN

669:    Input Parameters:
670: +  fn - the math function context
671: -  A  - matrix on which the function must be evaluated

673:    Output Parameter:
674: .  B  - (optional) matrix resulting from evaluating f(A)

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

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

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

690:    Level: advanced

692: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
693: @*/
694: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
695: {
697:   PetscBool      match,inplace=PETSC_FALSE;
698:   PetscInt       m,n,n1;

705:   if (B) {
708:   } else inplace = PETSC_TRUE;
709:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&match);
710:   if (!match) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Mat A must be of type seqdense");
711:   MatGetSize(A,&m,&n);
712:   if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %D rows, %D cols)",m,n);
713:   if (!inplace) {
714:     PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&match);
715:     if (!match) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Mat B must be of type seqdense");
716:     n1 = n;
717:     MatGetSize(B,&m,&n);
718:     if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat B is not square (has %D rows, %D cols)",m,n);
719:     if (n1!=n) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrices A and B must have the same dimension");
720:   }

722:   /* evaluate matrix function */
723:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
724:   FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE);
725:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
726:   return(0);
727: }

729: /*
730:    FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
731:    and then copies the first column.
732: */
733: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
734: {
736:   Mat            F;

739:   FN_AllocateWorkMat(fn,A,&F);
740:   if (fn->ops->evaluatefunctionmat[fn->method]) {
741:     (*fn->ops->evaluatefunctionmat[fn->method])(fn,A,F);
742:   } else if (!fn->method) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
743:   else SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
744:   MatGetColumnVector(F,v,0);
745:   FN_FreeWorkMat(fn,&F);
746:   return(0);
747: }

749: /*
750:    FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
751:    compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
752:    decomposition of A is A=Q*D*Q'. Only the first column is computed.
753: */
754: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
755: {
757:   PetscInt       m;
758:   PetscScalar    *As,*vs;

761:   MatDenseGetArray(A,&As);
762:   VecGetArray(v,&vs);
763:   MatGetSize(A,&m,NULL);
764:   FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE);
765:   MatDenseRestoreArray(A,&As);
766:   VecRestoreArray(v,&vs);
767:   return(0);
768: }

770: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
771: {
773:   PetscBool      set,flg,symm=PETSC_FALSE;
774:   PetscInt       m,n;
775:   Mat            M;
776:   PetscMPIInt    size,rank;
777:   PetscScalar    *pv;

780:   /* check symmetry of A */
781:   MatIsHermitianKnown(A,&set,&flg);
782:   symm = set? flg: PETSC_FALSE;

784:   /* evaluate matrix function */
785:   MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
786:   MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
787:   if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
788:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
789:     if (symm && !fn->method) {  /* prefer diagonalization */
790:       PetscInfo(fn,"Computing matrix function via diagonalization\n");
791:       FNEvaluateFunctionMatVec_Sym_Default(fn,A,v);
792:     } else {
793:       /* scale argument */
794:       if (fn->alpha!=(PetscScalar)1.0) {
795:         FN_AllocateWorkMat(fn,A,&M);
796:         MatScale(M,fn->alpha);
797:       } else M = A;
798:       if (fn->ops->evaluatefunctionmatvec[fn->method]) {
799:         (*fn->ops->evaluatefunctionmatvec[fn->method])(fn,M,v);
800:       } else {
801:         FNEvaluateFunctionMatVec_Default(fn,M,v);
802:       }
803:       if (fn->alpha!=(PetscScalar)1.0) {
804:         FN_FreeWorkMat(fn,&M);
805:       }
806:       /* scale result */
807:       VecScale(v,fn->beta);
808:     }
809:     PetscFPTrapPop();
810:   }

812:   /* synchronize */
813:   if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { 
814:     MatGetSize(A,&m,&n);
815:     VecGetArray(v,&pv);
816:     MPI_Bcast(pv,n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
817:     VecRestoreArray(v,&pv);
818:   }
819:   return(0);
820: }

822: /*@
823:    FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
824:    for a given matrix A.

826:    Logically Collective on FN

828:    Input Parameters:
829: +  fn - the math function context
830: -  A  - matrix on which the function must be evaluated

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

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

839:    Level: advanced

841: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
842: @*/
843: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
844: {
846:   PetscBool      match;
847:   PetscInt       m,n;

856:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&match);
857:   if (!match) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Mat A must be of type seqdense");
858:   MatGetSize(A,&m,&n);
859:   if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %D rows, %D cols)",m,n);
860:   VecGetSize(v,&m);
861:   if (m!=n) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrix A and vector v must have the same size");
862:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
863:   FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE);
864:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
865:   return(0);
866: }

868: /*@
869:    FNSetFromOptions - Sets FN options from the options database.

871:    Collective on FN

873:    Input Parameters:
874: .  fn - the math function context

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

879:    Level: beginner
880: @*/
881: PetscErrorCode FNSetFromOptions(FN fn)
882: {
884:   char           type[256];
885:   PetscScalar    array[2];
886:   PetscInt       k,meth;
887:   PetscBool      flg;
888:   FNParallelType pmode;

892:   FNRegisterAll();
893:   PetscObjectOptionsBegin((PetscObject)fn);
894:     PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,256,&flg);
895:     if (flg) {
896:       FNSetType(fn,type);
897:     } else if (!((PetscObject)fn)->type_name) {
898:       FNSetType(fn,FNRATIONAL);
899:     }

901:     k = 2;
902:     array[0] = 0.0; array[1] = 0.0;
903:     PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg);
904:     if (flg) {
905:       if (k<2) array[1] = 1.0;
906:       FNSetScale(fn,array[0],array[1]);
907:     }

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

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

915:     if (fn->ops->setfromoptions) {
916:       (*fn->ops->setfromoptions)(PetscOptionsObject,fn);
917:     }
918:     PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)fn);
919:   PetscOptionsEnd();
920:   return(0);
921: }

923: /*@C
924:    FNView - Prints the FN data structure.

926:    Collective on FN

928:    Input Parameters:
929: +  fn - the math function context
930: -  viewer - optional visualization context

932:    Note:
933:    The available visualization contexts include
934: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
935: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
936:          output where only the first processor opens
937:          the file.  All other processors send their
938:          data to the first processor to print.

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

943:    Level: beginner
944: @*/
945: PetscErrorCode FNView(FN fn,PetscViewer viewer)
946: {
947:   PetscBool      isascii;
948:   PetscInt       tabs;
950:   PetscMPIInt    size;

954:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)fn));
957:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
958:   if (isascii) {
959:     PetscViewerASCIIGetTab(viewer,&tabs);
960:     PetscViewerASCIISetTab(viewer,((PetscObject)fn)->tablevel);
961:     PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer);
962:     MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
963:     if (size>1) {
964:       PetscViewerASCIIPrintf(viewer,"  parallel operation mode: %s\n",FNParallelTypes[fn->pmode]);
965:     }
966:     if (fn->ops->view) {
967:       PetscViewerASCIIPushTab(viewer);
968:       (*fn->ops->view)(fn,viewer);
969:       PetscViewerASCIIPopTab(viewer);
970:     }
971:     PetscViewerASCIISetTab(viewer,tabs);
972:   }
973:   return(0);
974: }

976: /*@
977:    FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
978:    different communicator.

980:    Collective on FN

982:    Input Parameters:
983: +  fn   - the math function context
984: -  comm - MPI communicator

986:    Output Parameter:
987: .  newfn - location to put the new FN context

989:    Note:
990:    In order to use the same MPI communicator as in the original object,
991:    use PetscObjectComm((PetscObject)fn).

993:    Level: developer

995: .seealso: FNCreate()
996: @*/
997: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)
998: {
1000:   FNType         type;
1001:   PetscScalar    alpha,beta;
1002:   PetscInt       meth;
1003:   FNParallelType ptype;

1009:   FNCreate(comm,newfn);
1010:   FNGetType(fn,&type);
1011:   FNSetType(*newfn,type);
1012:   FNGetScale(fn,&alpha,&beta);
1013:   FNSetScale(*newfn,alpha,beta);
1014:   FNGetMethod(fn,&meth);
1015:   FNSetMethod(*newfn,meth);
1016:   FNGetParallel(fn,&ptype);
1017:   FNSetParallel(*newfn,ptype);
1018:   if (fn->ops->duplicate) {
1019:     (*fn->ops->duplicate)(fn,comm,newfn);
1020:   }
1021:   return(0);
1022: }

1024: /*@
1025:    FNDestroy - Destroys FN context that was created with FNCreate().

1027:    Collective on FN

1029:    Input Parameter:
1030: .  fn - the math function context

1032:    Level: beginner

1034: .seealso: FNCreate()
1035: @*/
1036: PetscErrorCode FNDestroy(FN *fn)
1037: {
1039:   PetscInt       i;

1042:   if (!*fn) return(0);
1044:   if (--((PetscObject)(*fn))->refct > 0) { *fn = 0; return(0); }
1045:   if ((*fn)->ops->destroy) { (*(*fn)->ops->destroy)(*fn); }
1046:   for (i=0;i<(*fn)->nw;i++) {
1047:     MatDestroy(&(*fn)->W[i]);
1048:   }
1049:   PetscHeaderDestroy(fn);
1050:   return(0);
1051: }

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

1056:    Not collective

1058:    Input Parameters:
1059: +  name - name of a new user-defined FN
1060: -  function - routine to create context

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

1065:    Level: advanced

1067: .seealso: FNRegisterAll()
1068: @*/
1069: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1070: {

1074:   FNInitializePackage();
1075:   PetscFunctionListAdd(&FNList,name,function);
1076:   return(0);
1077: }