Actual source code: fnbasic.c

slepc-3.18.0 2022-10-01
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 = NULL;
 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_",NULL};

 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:   return 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) return 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:   return 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 = NULL;
 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:   return 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:   return 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:   return 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:   return 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) return 0;

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

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

227:   PetscObjectChangeTypeName((PetscObject)fn,type);
228:   (*r)(fn);
229:   return 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:   return 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:   return 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:   return 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:   return 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:   return 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:   return 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:   return 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:   PetscUseTypeMethod(fn,evaluatefunction,xf,&yf);
462:   *y = fn->beta*yf;
463:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
464:   return 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:   PetscUseTypeMethod(fn,evaluatederivative,xf,&yf);
497:   *y = fn->alpha*fn->beta*yf;
498:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
499:   return 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:   PetscCallBLAS("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:   PetscCallBLAS("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:   PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
532: #else
533:   PetscCallBLAS("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:     PetscUseTypeMethod(fn,evaluatefunction,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:   PetscCallBLAS("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:   return 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:   return 0;
572: }

574: PetscErrorCode FNEvaluateFunctionMat_Basic(FN fn,Mat A,Mat F)
575: {
576:   PetscBool iscuda;

578:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda);
579:   if (iscuda && !fn->ops->evaluatefunctionmatcuda[fn->method]) PetscInfo(fn,"The method %" PetscInt_FMT " is not implemented for CUDA, falling back to CPU version\n",fn->method);
580:   if (iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatcuda[fn->method],A,F);
581:   else if (fn->ops->evaluatefunctionmat[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmat[fn->method],A,F);
582:   else {
585:   }
586:   return 0;
587: }

589: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
590: {
591:   PetscBool      set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
592:   PetscInt       m,n;
593:   PetscMPIInt    size,rank;
594:   PetscScalar    *pF;
595:   Mat            M,F;

597:   /* destination matrix */
598:   F = B?B:A;

600:   /* check symmetry of A */
601:   MatIsHermitianKnown(A,&set,&flg);
602:   symm = set? flg: PETSC_FALSE;

604:   MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
605:   MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
606:   if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
607:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
608:     PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda);
609:     hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
610:     if (!hasspecificmeth && symm && !fn->method) {  /* prefer diagonalization */
611:       PetscInfo(fn,"Computing matrix function via diagonalization\n");
612:       FNEvaluateFunctionMat_Sym_Default(fn,A,F);
613:     } else {
614:       /* scale argument */
615:       if (fn->alpha!=(PetscScalar)1.0) {
616:         FN_AllocateWorkMat(fn,A,&M);
617:         MatScale(M,fn->alpha);
618:       } else M = A;
619:       FNEvaluateFunctionMat_Basic(fn,M,F);
620:       if (fn->alpha!=(PetscScalar)1.0) FN_FreeWorkMat(fn,&M);
621:       /* scale result */
622:       MatScale(F,fn->beta);
623:     }
624:     PetscFPTrapPop();
625:   }
626:   if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {  /* synchronize */
627:     MatGetSize(A,&m,&n);
628:     MatDenseGetArray(F,&pF);
629:     MPI_Bcast(pF,n*n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
630:     MatDenseRestoreArray(F,&pF);
631:   }
632:   return 0;
633: }

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

639:    Logically Collective on fn

641:    Input Parameters:
642: +  fn - the math function context
643: -  A  - matrix on which the function must be evaluated

645:    Output Parameter:
646: .  B  - (optional) matrix resulting from evaluating f(A)

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

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

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

662:    Level: advanced

664: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
665: @*/
666: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
667: {
668:   PetscBool      inplace=PETSC_FALSE;
669:   PetscInt       m,n,n1;
670:   MatType        type;

676:   if (B) {
679:   } else inplace = PETSC_TRUE;
681:   MatGetSize(A,&m,&n);
683:   if (!inplace) {
684:     MatGetType(A,&type);
686:     n1 = n;
687:     MatGetSize(B,&m,&n);
690:   }

692:   /* evaluate matrix function */
693:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
694:   FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE);
695:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
696:   return 0;
697: }

699: /*
700:    FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
701:    and then copies the first column.
702: */
703: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
704: {
705:   Mat            F;

707:   FN_AllocateWorkMat(fn,A,&F);
708:   FNEvaluateFunctionMat_Basic(fn,A,F);
709:   MatGetColumnVector(F,v,0);
710:   FN_FreeWorkMat(fn,&F);
711:   return 0;
712: }

714: /*
715:    FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
716:    compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
717:    decomposition of A is A=Q*D*Q'. Only the first column is computed.
718: */
719: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
720: {
721:   PetscInt          m;
722:   const PetscScalar *As;
723:   PetscScalar       *vs;

725:   MatDenseGetArrayRead(A,&As);
726:   VecGetArray(v,&vs);
727:   MatGetSize(A,&m,NULL);
728:   FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE);
729:   MatDenseRestoreArrayRead(A,&As);
730:   VecRestoreArray(v,&vs);
731:   return 0;
732: }

734: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
735: {
736:   PetscBool      set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
737:   PetscInt       m,n;
738:   Mat            M;
739:   PetscMPIInt    size,rank;
740:   PetscScalar    *pv;

742:   /* check symmetry of A */
743:   MatIsHermitianKnown(A,&set,&flg);
744:   symm = set? flg: PETSC_FALSE;

746:   /* evaluate matrix function */
747:   MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
748:   MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
749:   if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
750:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
751:     PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda);
752:     hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
753:     if (!hasspecificmeth && symm && !fn->method) {  /* prefer diagonalization */
754:       PetscInfo(fn,"Computing matrix function via diagonalization\n");
755:       FNEvaluateFunctionMatVec_Sym_Default(fn,A,v);
756:     } else {
757:       /* scale argument */
758:       if (fn->alpha!=(PetscScalar)1.0) {
759:         FN_AllocateWorkMat(fn,A,&M);
760:         MatScale(M,fn->alpha);
761:       } else M = A;
762:       if (iscuda && fn->ops->evaluatefunctionmatveccuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatveccuda[fn->method],M,v);
763:       else if (fn->ops->evaluatefunctionmatvec[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatvec[fn->method],M,v);
764:       else FNEvaluateFunctionMatVec_Default(fn,M,v);
765:       if (fn->alpha!=(PetscScalar)1.0) FN_FreeWorkMat(fn,&M);
766:       /* scale result */
767:       VecScale(v,fn->beta);
768:     }
769:     PetscFPTrapPop();
770:   }

772:   /* synchronize */
773:   if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
774:     MatGetSize(A,&m,&n);
775:     VecGetArray(v,&pv);
776:     MPI_Bcast(pv,n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
777:     VecRestoreArray(v,&pv);
778:   }
779:   return 0;
780: }

782: /*@
783:    FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
784:    for a given matrix A.

786:    Logically Collective on fn

788:    Input Parameters:
789: +  fn - the math function context
790: -  A  - matrix on which the function must be evaluated

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

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

799:    Level: advanced

801: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
802: @*/
803: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
804: {
805:   PetscInt       m,n;
806:   PetscBool      iscuda;

815:   MatGetSize(A,&m,&n);
817:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda);
819:   VecGetSize(v,&m);
821:   PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
822:   FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE);
823:   PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
824:   return 0;
825: }

827: /*@
828:    FNSetFromOptions - Sets FN options from the options database.

830:    Collective on fn

832:    Input Parameters:
833: .  fn - the math function context

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

838:    Level: beginner

840: .seealso: FNSetOptionsPrefix()
841: @*/
842: PetscErrorCode FNSetFromOptions(FN fn)
843: {
844:   char           type[256];
845:   PetscScalar    array[2];
846:   PetscInt       k,meth;
847:   PetscBool      flg;
848:   FNParallelType pmode;

851:   FNRegisterAll();
852:   PetscObjectOptionsBegin((PetscObject)fn);
853:     PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg);
854:     if (flg) FNSetType(fn,type);
855:     else if (!((PetscObject)fn)->type_name) FNSetType(fn,FNRATIONAL);

857:     k = 2;
858:     array[0] = 0.0; array[1] = 0.0;
859:     PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg);
860:     if (flg) {
861:       if (k<2) array[1] = 1.0;
862:       FNSetScale(fn,array[0],array[1]);
863:     }

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

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

871:     PetscTryTypeMethod(fn,setfromoptions,PetscOptionsObject);
872:     PetscObjectProcessOptionsHandlers((PetscObject)fn,PetscOptionsObject);
873:   PetscOptionsEnd();
874:   return 0;
875: }

877: /*@C
878:    FNView - Prints the FN data structure.

880:    Collective on fn

882:    Input Parameters:
883: +  fn - the math function context
884: -  viewer - optional visualization context

886:    Note:
887:    The available visualization contexts include
888: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
889: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
890:          output where only the first processor opens
891:          the file.  All other processors send their
892:          data to the first processor to print.

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

897:    Level: beginner

899: .seealso: FNCreate()
900: @*/
901: PetscErrorCode FNView(FN fn,PetscViewer viewer)
902: {
903:   PetscBool      isascii;
904:   PetscMPIInt    size;

907:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer);
910:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
911:   if (isascii) {
912:     PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer);
913:     MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
914:     if (size>1) PetscViewerASCIIPrintf(viewer,"  parallel operation mode: %s\n",FNParallelTypes[fn->pmode]);
915:     PetscViewerASCIIPushTab(viewer);
916:     PetscTryTypeMethod(fn,view,viewer);
917:     PetscViewerASCIIPopTab(viewer);
918:   }
919:   return 0;
920: }

922: /*@C
923:    FNViewFromOptions - View from options

925:    Collective on FN

927:    Input Parameters:
928: +  fn   - the math function context
929: .  obj  - optional object
930: -  name - command line option

932:    Level: intermediate

934: .seealso: FNView(), FNCreate()
935: @*/
936: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])
937: {
939:   PetscObjectViewFromOptions((PetscObject)fn,obj,name);
940:   return 0;
941: }

943: /*@
944:    FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
945:    different communicator.

947:    Collective on fn

949:    Input Parameters:
950: +  fn   - the math function context
951: -  comm - MPI communicator

953:    Output Parameter:
954: .  newfn - location to put the new FN context

956:    Note:
957:    In order to use the same MPI communicator as in the original object,
958:    use PetscObjectComm((PetscObject)fn).

960:    Level: developer

962: .seealso: FNCreate()
963: @*/
964: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)
965: {
966:   FNType         type;
967:   PetscScalar    alpha,beta;
968:   PetscInt       meth;
969:   FNParallelType ptype;

974:   FNCreate(comm,newfn);
975:   FNGetType(fn,&type);
976:   FNSetType(*newfn,type);
977:   FNGetScale(fn,&alpha,&beta);
978:   FNSetScale(*newfn,alpha,beta);
979:   FNGetMethod(fn,&meth);
980:   FNSetMethod(*newfn,meth);
981:   FNGetParallel(fn,&ptype);
982:   FNSetParallel(*newfn,ptype);
983:   PetscTryTypeMethod(fn,duplicate,comm,newfn);
984:   return 0;
985: }

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

990:    Collective on fn

992:    Input Parameter:
993: .  fn - the math function context

995:    Level: beginner

997: .seealso: FNCreate()
998: @*/
999: PetscErrorCode FNDestroy(FN *fn)
1000: {
1001:   PetscInt       i;

1003:   if (!*fn) return 0;
1005:   if (--((PetscObject)(*fn))->refct > 0) { *fn = NULL; return 0; }
1006:   PetscTryTypeMethod(*fn,destroy);
1007:   for (i=0;i<(*fn)->nw;i++) MatDestroy(&(*fn)->W[i]);
1008:   PetscHeaderDestroy(fn);
1009:   return 0;
1010: }

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

1015:    Not collective

1017:    Input Parameters:
1018: +  name - name of a new user-defined FN
1019: -  function - routine to create context

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

1024:    Level: advanced

1026: .seealso: FNRegisterAll()
1027: @*/
1028: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1029: {
1030:   FNInitializePackage();
1031:   PetscFunctionListAdd(&FNList,name,function);
1032:   return 0;
1033: }