Actual source code: fnbasic.c

slepc-3.22.1 2024-10-28
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:   PetscFunctionBegin;
 36:   PetscCall(PetscFunctionListDestroy(&FNList));
 37:   FNPackageInitialized = PETSC_FALSE;
 38:   FNRegisterAllCalled  = PETSC_FALSE;
 39:   PetscFunctionReturn(PETSC_SUCCESS);
 40: }

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

 47:   Level: developer

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

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

 80: /*@
 81:    FNCreate - Creates an FN context.

 83:    Collective

 85:    Input Parameter:
 86: .  comm - MPI communicator

 88:    Output Parameter:
 89: .  newfn - location to put the FN context

 91:    Level: beginner

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

 99:   PetscFunctionBegin;
100:   PetscAssertPointer(newfn,2);
101:   PetscCall(FNInitializePackage());
102:   PetscCall(SlepcHeaderCreate(fn,FN_CLASSID,"FN","Math Function","FN",comm,FNDestroy,FNView));

104:   fn->alpha    = 1.0;
105:   fn->beta     = 1.0;
106:   fn->method   = 0;

108:   fn->nw       = 0;
109:   fn->cw       = 0;
110:   fn->data     = NULL;

112:   *newfn = fn;
113:   PetscFunctionReturn(PETSC_SUCCESS);
114: }

116: /*@
117:    FNSetOptionsPrefix - Sets the prefix used for searching for all
118:    FN options in the database.

120:    Logically Collective

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

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

131:    Level: advanced

133: .seealso: FNAppendOptionsPrefix()
134: @*/
135: PetscErrorCode FNSetOptionsPrefix(FN fn,const char *prefix)
136: {
137:   PetscFunctionBegin;
139:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fn,prefix));
140:   PetscFunctionReturn(PETSC_SUCCESS);
141: }

143: /*@
144:    FNAppendOptionsPrefix - Appends to the prefix used for searching for all
145:    FN options in the database.

147:    Logically Collective

149:    Input Parameters:
150: +  fn - the math function context
151: -  prefix - the prefix string to prepend to all FN option requests

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

157:    Level: advanced

159: .seealso: FNSetOptionsPrefix()
160: @*/
161: PetscErrorCode FNAppendOptionsPrefix(FN fn,const char *prefix)
162: {
163:   PetscFunctionBegin;
165:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)fn,prefix));
166:   PetscFunctionReturn(PETSC_SUCCESS);
167: }

169: /*@
170:    FNGetOptionsPrefix - Gets the prefix used for searching for all
171:    FN options in the database.

173:    Not Collective

175:    Input Parameters:
176: .  fn - the math function context

178:    Output Parameters:
179: .  prefix - pointer to the prefix string used is returned

181:    Note:
182:    On the Fortran side, the user should pass in a string 'prefix' of
183:    sufficient length to hold the prefix.

185:    Level: advanced

187: .seealso: FNSetOptionsPrefix(), FNAppendOptionsPrefix()
188: @*/
189: PetscErrorCode FNGetOptionsPrefix(FN fn,const char *prefix[])
190: {
191:   PetscFunctionBegin;
193:   PetscAssertPointer(prefix,2);
194:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)fn,prefix));
195:   PetscFunctionReturn(PETSC_SUCCESS);
196: }

198: /*@
199:    FNSetType - Selects the type for the FN object.

201:    Logically Collective

203:    Input Parameters:
204: +  fn   - the math function context
205: -  type - a known type

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

211:    Level: intermediate

213: .seealso: FNGetType()
214: @*/
215: PetscErrorCode FNSetType(FN fn,FNType type)
216: {
217:   PetscErrorCode (*r)(FN);
218:   PetscBool      match;

220:   PetscFunctionBegin;
222:   PetscAssertPointer(type,2);

224:   PetscCall(PetscObjectTypeCompare((PetscObject)fn,type,&match));
225:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

227:   PetscCall(PetscFunctionListFind(FNList,type,&r));
228:   PetscCheck(r,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested FN type %s",type);

230:   PetscTryTypeMethod(fn,destroy);
231:   PetscCall(PetscMemzero(fn->ops,sizeof(struct _FNOps)));

233:   PetscCall(PetscObjectChangeTypeName((PetscObject)fn,type));
234:   PetscCall((*r)(fn));
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

238: /*@
239:    FNGetType - Gets the FN type name (as a string) from the FN context.

241:    Not Collective

243:    Input Parameter:
244: .  fn - the math function context

246:    Output Parameter:
247: .  type - name of the math function

249:    Level: intermediate

251: .seealso: FNSetType()
252: @*/
253: PetscErrorCode FNGetType(FN fn,FNType *type)
254: {
255:   PetscFunctionBegin;
257:   PetscAssertPointer(type,2);
258:   *type = ((PetscObject)fn)->type_name;
259:   PetscFunctionReturn(PETSC_SUCCESS);
260: }

262: /*@
263:    FNSetScale - Sets the scaling parameters that define the matematical function.

265:    Logically Collective

267:    Input Parameters:
268: +  fn    - the math function context
269: .  alpha - inner scaling (argument)
270: -  beta  - outer scaling (result)

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

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

282:    Level: intermediate

284: .seealso: FNGetScale(), FNEvaluateFunction()
285: @*/
286: PetscErrorCode FNSetScale(FN fn,PetscScalar alpha,PetscScalar beta)
287: {
288:   PetscFunctionBegin;
292:   PetscCheck(PetscAbsScalar(alpha)!=0.0 && PetscAbsScalar(beta)!=0.0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_WRONG,"Scaling factors must be nonzero");
293:   fn->alpha = alpha;
294:   fn->beta  = beta;
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: /*@
299:    FNGetScale - Gets the scaling parameters that define the matematical function.

301:    Not Collective

303:    Input Parameter:
304: .  fn    - the math function context

306:    Output Parameters:
307: +  alpha - inner scaling (argument)
308: -  beta  - outer scaling (result)

310:    Level: intermediate

312: .seealso: FNSetScale()
313: @*/
314: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)
315: {
316:   PetscFunctionBegin;
318:   if (alpha) *alpha = fn->alpha;
319:   if (beta)  *beta  = fn->beta;
320:   PetscFunctionReturn(PETSC_SUCCESS);
321: }

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

326:    Logically Collective

328:    Input Parameters:
329: +  fn   - the math function context
330: -  meth - an index identifying the method

332:    Options Database Key:
333: .  -fn_method <meth> - Sets the method

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

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

343:    Level: intermediate

345: .seealso: FNGetMethod(), FNEvaluateFunctionMat()
346: @*/
347: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)
348: {
349:   PetscFunctionBegin;
352:   PetscCheck(meth>=0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
353:   PetscCheck(meth<=FN_MAX_SOLVE,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
354:   fn->method = meth;
355:   PetscFunctionReturn(PETSC_SUCCESS);
356: }

358: /*@
359:    FNGetMethod - Gets the method currently used in the FN.

361:    Not Collective

363:    Input Parameter:
364: .  fn - the math function context

366:    Output Parameter:
367: .  meth - identifier of the method

369:    Level: intermediate

371: .seealso: FNSetMethod()
372: @*/
373: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)
374: {
375:   PetscFunctionBegin;
377:   PetscAssertPointer(meth,2);
378:   *meth = fn->method;
379:   PetscFunctionReturn(PETSC_SUCCESS);
380: }

382: /*@
383:    FNSetParallel - Selects the mode of operation in parallel runs.

385:    Logically Collective

387:    Input Parameters:
388: +  fn    - the math function context
389: -  pmode - the parallel mode

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

394:    Notes:
395:    This is relevant only when the function is evaluated on a matrix, with
396:    either FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().

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

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

408:    Level: advanced

410: .seealso: FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec(), FNGetParallel()
411: @*/
412: PetscErrorCode FNSetParallel(FN fn,FNParallelType pmode)
413: {
414:   PetscFunctionBegin;
417:   fn->pmode = pmode;
418:   PetscFunctionReturn(PETSC_SUCCESS);
419: }

421: /*@
422:    FNGetParallel - Gets the mode of operation in parallel runs.

424:    Not Collective

426:    Input Parameter:
427: .  fn - the math function context

429:    Output Parameter:
430: .  pmode - the parallel mode

432:    Level: advanced

434: .seealso: FNSetParallel()
435: @*/
436: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)
437: {
438:   PetscFunctionBegin;
440:   PetscAssertPointer(pmode,2);
441:   *pmode = fn->pmode;
442:   PetscFunctionReturn(PETSC_SUCCESS);
443: }

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

448:    Not Collective

450:    Input Parameters:
451: +  fn - the math function context
452: -  x  - the value where the function must be evaluated

454:    Output Parameter:
455: .  y  - the result of f(x)

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

461:    Level: intermediate

463: .seealso: FNEvaluateDerivative(), FNEvaluateFunctionMat(), FNSetScale()
464: @*/
465: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)
466: {
467:   PetscScalar    xf,yf;

469:   PetscFunctionBegin;
472:   PetscAssertPointer(y,3);
473:   PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
474:   xf = fn->alpha*x;
475:   PetscUseTypeMethod(fn,evaluatefunction,xf,&yf);
476:   *y = fn->beta*yf;
477:   PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

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

484:    Not Collective

486:    Input Parameters:
487: +  fn - the math function context
488: -  x  - the value where the derivative must be evaluated

490:    Output Parameter:
491: .  y  - the result of f'(x)

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

497:    Level: intermediate

499: .seealso: FNEvaluateFunction()
500: @*/
501: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)
502: {
503:   PetscScalar    xf,yf;

505:   PetscFunctionBegin;
508:   PetscAssertPointer(y,3);
509:   PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
510:   xf = fn->alpha*x;
511:   PetscUseTypeMethod(fn,evaluatederivative,xf,&yf);
512:   *y = fn->alpha*fn->beta*yf;
513:   PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
514:   PetscFunctionReturn(PETSC_SUCCESS);
515: }

517: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,const PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)
518: {
519:   PetscInt       i,j;
520:   PetscBLASInt   n,k,ld,lwork,info;
521:   PetscScalar    *Q,*W,*work,adummy,a,x,y,one=1.0,zero=0.0;
522:   PetscReal      *eig,dummy;
523: #if defined(PETSC_USE_COMPLEX)
524:   PetscReal      *rwork,rdummy;
525: #endif

527:   PetscFunctionBegin;
528:   PetscCall(PetscBLASIntCast(m,&n));
529:   ld = n;
530:   k = firstonly? 1: n;

532:   /* workspace query and memory allocation */
533:   lwork = -1;
534: #if defined(PETSC_USE_COMPLEX)
535:   PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&rdummy,&info));
536:   PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
537:   PetscCall(PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork));
538: #else
539:   PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&info));
540:   PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
541:   PetscCall(PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work));
542: #endif

544:   /* compute eigendecomposition */
545:   for (j=0;j<n;j++) for (i=j;i<n;i++) Q[i+j*ld] = As[i+j*ld];
546: #if defined(PETSC_USE_COMPLEX)
547:   PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
548: #else
549:   PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
550: #endif
551:   SlepcCheckLapackInfo("syev",info);

553:   /* W = f(Lambda)*Q' */
554:   for (i=0;i<n;i++) {
555:     x = fn->alpha*eig[i];
556:     PetscUseTypeMethod(fn,evaluatefunction,x,&y);  /* y = f(x) */
557:     for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
558:   }
559:   /* Bs = Q*W */
560:   PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
561: #if defined(PETSC_USE_COMPLEX)
562:   PetscCall(PetscFree5(eig,Q,W,work,rwork));
563: #else
564:   PetscCall(PetscFree4(eig,Q,W,work));
565: #endif
566:   PetscCall(PetscLogFlops(9.0*n*n*n+2.0*n*n*n));
567:   PetscFunctionReturn(PETSC_SUCCESS);
568: }

570: /*
571:    FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
572:    compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
573:    decomposition of A is A=Q*D*Q'
574: */
575: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)
576: {
577:   PetscInt          m;
578:   const PetscScalar *As;
579:   PetscScalar       *Bs;

581:   PetscFunctionBegin;
582:   PetscCall(MatDenseGetArrayRead(A,&As));
583:   PetscCall(MatDenseGetArray(B,&Bs));
584:   PetscCall(MatGetSize(A,&m,NULL));
585:   PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE));
586:   PetscCall(MatDenseRestoreArrayRead(A,&As));
587:   PetscCall(MatDenseRestoreArray(B,&Bs));
588:   PetscFunctionReturn(PETSC_SUCCESS);
589: }

591: static PetscErrorCode FNEvaluateFunctionMat_Basic(FN fn,Mat A,Mat F)
592: {
593:   PetscBool iscuda;

595:   PetscFunctionBegin;
596:   PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
597:   if (iscuda && !fn->ops->evaluatefunctionmatcuda[fn->method]) PetscCall(PetscInfo(fn,"The method %" PetscInt_FMT " is not implemented for CUDA, falling back to CPU version\n",fn->method));
598:   if (iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatcuda[fn->method],A,F);
599:   else if (fn->ops->evaluatefunctionmat[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmat[fn->method],A,F);
600:   else {
601:     PetscCheck(fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
602:     PetscCheck(!fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
603:   }
604:   PetscFunctionReturn(PETSC_SUCCESS);
605: }

607: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
608: {
609:   PetscBool      set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
610:   PetscInt       m,n;
611:   PetscMPIInt    size,rank,n2;
612:   PetscScalar    *pF;
613:   Mat            M,F;

615:   PetscFunctionBegin;
616:   /* destination matrix */
617:   F = B?B:A;

619:   /* check symmetry of A */
620:   PetscCall(MatIsHermitianKnown(A,&set,&flg));
621:   symm = set? flg: PETSC_FALSE;

623:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
624:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
625:   if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
626:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
627:     PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
628:     hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
629:     if (!hasspecificmeth && symm && !fn->method) {  /* prefer diagonalization */
630:       PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
631:       PetscCall(FNEvaluateFunctionMat_Sym_Default(fn,A,F));
632:     } else {
633:       /* scale argument */
634:       if (fn->alpha!=(PetscScalar)1.0) {
635:         PetscCall(FN_AllocateWorkMat(fn,A,&M));
636:         PetscCall(MatScale(M,fn->alpha));
637:       } else M = A;
638:       PetscCall(FNEvaluateFunctionMat_Basic(fn,M,F));
639:       if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
640:       /* scale result */
641:       PetscCall(MatScale(F,fn->beta));
642:     }
643:     PetscCall(PetscFPTrapPop());
644:   }
645:   if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {  /* synchronize */
646:     PetscCall(MatGetSize(A,&m,&n));
647:     PetscCall(MatDenseGetArray(F,&pF));
648:     PetscCall(PetscMPIIntCast(n*n,&n2));
649:     PetscCallMPI(MPI_Bcast(pF,n2,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
650:     PetscCall(MatDenseRestoreArray(F,&pF));
651:   }
652:   PetscFunctionReturn(PETSC_SUCCESS);
653: }

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

659:    Logically Collective

661:    Input Parameters:
662: +  fn - the math function context
663: -  A  - matrix on which the function must be evaluated

665:    Output Parameter:
666: .  B  - (optional) matrix resulting from evaluating f(A)

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

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

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

682:    Level: advanced

684: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
685: @*/
686: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
687: {
688:   PetscBool      inplace=PETSC_FALSE;
689:   PetscInt       m,n,n1;
690:   MatType        type;

692:   PetscFunctionBegin;
697:   if (B) {
700:   } else inplace = PETSC_TRUE;
701:   PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA);   //SlepcMatCheckSeq(A);
702:   PetscCall(MatGetSize(A,&m,&n));
703:   PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
704:   if (!inplace) {
705:     PetscCall(MatGetType(A,&type));
706:     PetscCheckTypeName(B,type);
707:     n1 = n;
708:     PetscCall(MatGetSize(B,&m,&n));
709:     PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat B is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
710:     PetscCheck(n1==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrices A and B must have the same dimension");
711:   }

713:   /* evaluate matrix function */
714:   PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
715:   PetscCall(FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE));
716:   PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
717:   PetscFunctionReturn(PETSC_SUCCESS);
718: }

720: /*
721:    FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
722:    and then copies the first column.
723: */
724: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
725: {
726:   Mat            F;

728:   PetscFunctionBegin;
729:   PetscCall(FN_AllocateWorkMat(fn,A,&F));
730:   PetscCall(FNEvaluateFunctionMat_Basic(fn,A,F));
731:   PetscCall(MatGetColumnVector(F,v,0));
732:   PetscCall(FN_FreeWorkMat(fn,&F));
733:   PetscFunctionReturn(PETSC_SUCCESS);
734: }

736: /*
737:    FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
738:    compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
739:    decomposition of A is A=Q*D*Q'. Only the first column is computed.
740: */
741: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
742: {
743:   PetscInt          m;
744:   const PetscScalar *As;
745:   PetscScalar       *vs;

747:   PetscFunctionBegin;
748:   PetscCall(MatDenseGetArrayRead(A,&As));
749:   PetscCall(VecGetArray(v,&vs));
750:   PetscCall(MatGetSize(A,&m,NULL));
751:   PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE));
752:   PetscCall(MatDenseRestoreArrayRead(A,&As));
753:   PetscCall(VecRestoreArray(v,&vs));
754:   PetscFunctionReturn(PETSC_SUCCESS);
755: }

757: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
758: {
759:   PetscBool      set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
760:   PetscInt       m,n;
761:   Mat            M;
762:   PetscMPIInt    size,rank,n_;
763:   PetscScalar    *pv;

765:   PetscFunctionBegin;
766:   /* check symmetry of A */
767:   PetscCall(MatIsHermitianKnown(A,&set,&flg));
768:   symm = set? flg: PETSC_FALSE;

770:   /* evaluate matrix function */
771:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
772:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
773:   if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
774:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
775:     PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
776:     hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
777:     if (!hasspecificmeth && symm && !fn->method) {  /* prefer diagonalization */
778:       PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
779:       PetscCall(FNEvaluateFunctionMatVec_Sym_Default(fn,A,v));
780:     } else {
781:       /* scale argument */
782:       if (fn->alpha!=(PetscScalar)1.0) {
783:         PetscCall(FN_AllocateWorkMat(fn,A,&M));
784:         PetscCall(MatScale(M,fn->alpha));
785:       } else M = A;
786:       if (iscuda && fn->ops->evaluatefunctionmatveccuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatveccuda[fn->method],M,v);
787:       else if (fn->ops->evaluatefunctionmatvec[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatvec[fn->method],M,v);
788:       else PetscCall(FNEvaluateFunctionMatVec_Default(fn,M,v));
789:       if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
790:       /* scale result */
791:       PetscCall(VecScale(v,fn->beta));
792:     }
793:     PetscCall(PetscFPTrapPop());
794:   }

796:   /* synchronize */
797:   if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
798:     PetscCall(MatGetSize(A,&m,&n));
799:     PetscCall(VecGetArray(v,&pv));
800:     PetscCall(PetscMPIIntCast(n,&n_));
801:     PetscCallMPI(MPI_Bcast(pv,n_,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
802:     PetscCall(VecRestoreArray(v,&pv));
803:   }
804:   PetscFunctionReturn(PETSC_SUCCESS);
805: }

807: /*@
808:    FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
809:    for a given matrix A.

811:    Logically Collective

813:    Input Parameters:
814: +  fn - the math function context
815: -  A  - matrix on which the function must be evaluated

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

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

824:    Level: advanced

826: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
827: @*/
828: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
829: {
830:   PetscInt       m,n;
831:   PetscBool      iscuda;

833:   PetscFunctionBegin;
840:   PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA);   //SlepcMatCheckSeq(A);
841:   PetscCall(MatGetSize(A,&m,&n));
842:   PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
843:   PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
844:   PetscCheckTypeName(v,iscuda?VECSEQCUDA:VECSEQ);
845:   PetscCall(VecGetSize(v,&m));
846:   PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrix A and vector v must have the same size");
847:   PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
848:   PetscCall(FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE));
849:   PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
850:   PetscFunctionReturn(PETSC_SUCCESS);
851: }

853: /*@
854:    FNSetFromOptions - Sets FN options from the options database.

856:    Collective

858:    Input Parameters:
859: .  fn - the math function context

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

864:    Level: beginner

866: .seealso: FNSetOptionsPrefix()
867: @*/
868: PetscErrorCode FNSetFromOptions(FN fn)
869: {
870:   char           type[256];
871:   PetscScalar    array[2];
872:   PetscInt       k,meth;
873:   PetscBool      flg;
874:   FNParallelType pmode;

876:   PetscFunctionBegin;
878:   PetscCall(FNRegisterAll());
879:   PetscObjectOptionsBegin((PetscObject)fn);
880:     PetscCall(PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg));
881:     if (flg) PetscCall(FNSetType(fn,type));
882:     else if (!((PetscObject)fn)->type_name) PetscCall(FNSetType(fn,FNRATIONAL));

884:     k = 2;
885:     array[0] = 0.0; array[1] = 0.0;
886:     PetscCall(PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg));
887:     if (flg) {
888:       if (k<2) array[1] = 1.0;
889:       PetscCall(FNSetScale(fn,array[0],array[1]));
890:     }

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

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

898:     PetscTryTypeMethod(fn,setfromoptions,PetscOptionsObject);
899:     PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fn,PetscOptionsObject));
900:   PetscOptionsEnd();
901:   PetscFunctionReturn(PETSC_SUCCESS);
902: }

904: /*@
905:    FNView - Prints the FN data structure.

907:    Collective

909:    Input Parameters:
910: +  fn - the math function context
911: -  viewer - optional visualization context

913:    Note:
914:    The available visualization contexts include
915: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
916: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
917:          output where only the first processor opens
918:          the file.  All other processors send their
919:          data to the first processor to print.

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

924:    Level: beginner

926: .seealso: FNCreate()
927: @*/
928: PetscErrorCode FNView(FN fn,PetscViewer viewer)
929: {
930:   PetscBool      isascii;
931:   PetscMPIInt    size;

933:   PetscFunctionBegin;
935:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer));
937:   PetscCheckSameComm(fn,1,viewer,2);
938:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
939:   if (isascii) {
940:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer));
941:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
942:     if (size>1) PetscCall(PetscViewerASCIIPrintf(viewer,"  parallel operation mode: %s\n",FNParallelTypes[fn->pmode]));
943:     PetscCall(PetscViewerASCIIPushTab(viewer));
944:     PetscTryTypeMethod(fn,view,viewer);
945:     PetscCall(PetscViewerASCIIPopTab(viewer));
946:   }
947:   PetscFunctionReturn(PETSC_SUCCESS);
948: }

950: /*@
951:    FNViewFromOptions - View from options

953:    Collective

955:    Input Parameters:
956: +  fn   - the math function context
957: .  obj  - optional object
958: -  name - command line option

960:    Level: intermediate

962: .seealso: FNView(), FNCreate()
963: @*/
964: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])
965: {
966:   PetscFunctionBegin;
968:   PetscCall(PetscObjectViewFromOptions((PetscObject)fn,obj,name));
969:   PetscFunctionReturn(PETSC_SUCCESS);
970: }

972: /*@
973:    FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
974:    different communicator.

976:    Collective

978:    Input Parameters:
979: +  fn   - the math function context
980: -  comm - MPI communicator

982:    Output Parameter:
983: .  newfn - location to put the new FN context

985:    Note:
986:    In order to use the same MPI communicator as in the original object,
987:    use PetscObjectComm((PetscObject)fn).

989:    Level: developer

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

1000:   PetscFunctionBegin;
1003:   PetscAssertPointer(newfn,3);
1004:   PetscCall(FNCreate(comm,newfn));
1005:   PetscCall(FNGetType(fn,&type));
1006:   PetscCall(FNSetType(*newfn,type));
1007:   PetscCall(FNGetScale(fn,&alpha,&beta));
1008:   PetscCall(FNSetScale(*newfn,alpha,beta));
1009:   PetscCall(FNGetMethod(fn,&meth));
1010:   PetscCall(FNSetMethod(*newfn,meth));
1011:   PetscCall(FNGetParallel(fn,&ptype));
1012:   PetscCall(FNSetParallel(*newfn,ptype));
1013:   PetscTryTypeMethod(fn,duplicate,comm,newfn);
1014:   PetscFunctionReturn(PETSC_SUCCESS);
1015: }

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

1020:    Collective

1022:    Input Parameter:
1023: .  fn - the math function context

1025:    Level: beginner

1027: .seealso: FNCreate()
1028: @*/
1029: PetscErrorCode FNDestroy(FN *fn)
1030: {
1031:   PetscInt       i;

1033:   PetscFunctionBegin;
1034:   if (!*fn) PetscFunctionReturn(PETSC_SUCCESS);
1036:   if (--((PetscObject)*fn)->refct > 0) { *fn = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
1037:   PetscTryTypeMethod(*fn,destroy);
1038:   for (i=0;i<(*fn)->nw;i++) PetscCall(MatDestroy(&(*fn)->W[i]));
1039:   PetscCall(PetscHeaderDestroy(fn));
1040:   PetscFunctionReturn(PETSC_SUCCESS);
1041: }

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

1046:    Not Collective

1048:    Input Parameters:
1049: +  name - name of a new user-defined FN
1050: -  function - routine to create context

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

1055:    Level: advanced

1057: .seealso: FNRegisterAll()
1058: @*/
1059: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1060: {
1061:   PetscFunctionBegin;
1062:   PetscCall(FNInitializePackage());
1063:   PetscCall(PetscFunctionListAdd(&FNList,name,function));
1064:   PetscFunctionReturn(PETSC_SUCCESS);
1065: }