Actual source code: fnbasic.c
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: /*@
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: [](sec:fn), `SlepcFinalize()`, `FNInitializePackage()`
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: /*@
43: FNInitializePackage - This function initializes everything in the `FN` package.
44: It is called from `PetscDLLibraryRegister_slepc()` when using dynamic libraries, and
45: on the first call to `FNCreate()` when using shared or static libraries.
47: Note:
48: This function never needs to be called by SLEPc users.
50: Level: developer
52: .seealso: [](sec:fn), `FN`, `SlepcInitialize()`, `FNFinalizePackage()`
53: @*/
54: PetscErrorCode FNInitializePackage(void)
55: {
56: char logList[256];
57: PetscBool opt,pkg;
58: PetscClassId classids[1];
60: PetscFunctionBegin;
61: if (FNPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
62: FNPackageInitialized = PETSC_TRUE;
63: /* Register Classes */
64: PetscCall(PetscClassIdRegister("Math Function",&FN_CLASSID));
65: /* Register Constructors */
66: PetscCall(FNRegisterAll());
67: /* Register Events */
68: PetscCall(PetscLogEventRegister("FNEvaluate",FN_CLASSID,&FN_Evaluate));
69: /* Process Info */
70: classids[0] = FN_CLASSID;
71: PetscCall(PetscInfoProcessClass("fn",1,&classids[0]));
72: /* Process summary exclusions */
73: PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
74: if (opt) {
75: PetscCall(PetscStrInList("fn",logList,',',&pkg));
76: if (pkg) PetscCall(PetscLogEventDeactivateClass(FN_CLASSID));
77: }
78: /* Register package finalizer */
79: PetscCall(PetscRegisterFinalize(FNFinalizePackage));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: /*@
84: FNCreate - Creates an `FN` context.
86: Collective
88: Input Parameter:
89: . comm - MPI communicator
91: Output Parameter:
92: . newfn - location to put the `FN` context
94: Level: beginner
96: .seealso: [](sec:fn), `FNDestroy()`, `FN`
97: @*/
98: PetscErrorCode FNCreate(MPI_Comm comm,FN *newfn)
99: {
100: FN fn;
102: PetscFunctionBegin;
103: PetscAssertPointer(newfn,2);
104: PetscCall(FNInitializePackage());
105: PetscCall(SlepcHeaderCreate(fn,FN_CLASSID,"FN","Math Function","FN",comm,FNDestroy,FNView));
107: fn->alpha = 1.0;
108: fn->beta = 1.0;
109: fn->method = 0;
111: fn->nw = 0;
112: fn->cw = 0;
113: fn->data = NULL;
115: *newfn = fn;
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: /*@
120: FNSetOptionsPrefix - Sets the prefix used for searching for all
121: `FN` options in the database.
123: Logically Collective
125: Input Parameters:
126: + fn - the math function context
127: - prefix - the prefix string to prepend to all `FN` option requests
129: Notes:
130: A hyphen (-) must NOT be given at the beginning of the prefix name.
131: The first character of all runtime options is AUTOMATICALLY the
132: hyphen.
134: Level: advanced
136: .seealso: [](sec:fn), `FNAppendOptionsPrefix()`
137: @*/
138: PetscErrorCode FNSetOptionsPrefix(FN fn,const char prefix[])
139: {
140: PetscFunctionBegin;
142: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fn,prefix));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /*@
147: FNAppendOptionsPrefix - Appends to the prefix used for searching for all
148: `FN` options in the database.
150: Logically Collective
152: Input Parameters:
153: + fn - the math function context
154: - prefix - the prefix string to prepend to all `FN` option requests
156: Notes:
157: A hyphen (-) must NOT be given at the beginning of the prefix name.
158: The first character of all runtime options is AUTOMATICALLY the hyphen.
160: Level: advanced
162: .seealso: [](sec:fn), `FNSetOptionsPrefix()`
163: @*/
164: PetscErrorCode FNAppendOptionsPrefix(FN fn,const char prefix[])
165: {
166: PetscFunctionBegin;
168: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)fn,prefix));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: /*@
173: FNGetOptionsPrefix - Gets the prefix used for searching for all
174: `FN` options in the database.
176: Not Collective
178: Input Parameter:
179: . fn - the math function context
181: Output Parameter:
182: . prefix - pointer to the prefix string used is returned
184: Level: advanced
186: .seealso: [](sec:fn), `FNSetOptionsPrefix()`, `FNAppendOptionsPrefix()`
187: @*/
188: PetscErrorCode FNGetOptionsPrefix(FN fn,const char *prefix[])
189: {
190: PetscFunctionBegin;
192: PetscAssertPointer(prefix,2);
193: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)fn,prefix));
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: /*@
198: FNSetType - Selects the type for the `FN` object.
200: Logically Collective
202: Input Parameters:
203: + fn - the math function context
204: - type - a known type
206: Options Database Key:
207: . -fn_type type - sets the `FN` type
209: Note:
210: The default is `FNRATIONAL`, which includes polynomials as a particular
211: case as well as simple functions such as $f(x)=x$ and $f(x)=constant$.
213: Level: intermediate
215: .seealso: [](sec:fn), `FNGetType()`
216: @*/
217: PetscErrorCode FNSetType(FN fn,FNType type)
218: {
219: PetscErrorCode (*r)(FN);
220: PetscBool match;
222: PetscFunctionBegin;
224: PetscAssertPointer(type,2);
226: PetscCall(PetscObjectTypeCompare((PetscObject)fn,type,&match));
227: if (match) PetscFunctionReturn(PETSC_SUCCESS);
229: PetscCall(PetscFunctionListFind(FNList,type,&r));
230: PetscCheck(r,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested FN type %s",type);
232: PetscTryTypeMethod(fn,destroy);
233: PetscCall(PetscMemzero(fn->ops,sizeof(struct _FNOps)));
235: PetscCall(PetscObjectChangeTypeName((PetscObject)fn,type));
236: PetscCall((*r)(fn));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*@
241: FNGetType - Gets the `FN` type name (as a string) from the `FN` context.
243: Not Collective
245: Input Parameter:
246: . fn - the math function context
248: Output Parameter:
249: . type - name of the math function
251: Note:
252: `type` should not be retained for later use as it will be an invalid pointer
253: if the `FNType` of `fn` is changed.
255: Level: intermediate
257: .seealso: [](sec:fn), `FNSetType()`, `PetscObjectTypeCompare()`, `PetscObjectTypeCompareAny()`
258: @*/
259: PetscErrorCode FNGetType(FN fn,FNType *type)
260: {
261: PetscFunctionBegin;
263: PetscAssertPointer(type,2);
264: *type = ((PetscObject)fn)->type_name;
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: /*@
269: FNSetScale - Sets the scaling parameters that define the matematical function.
271: Logically Collective
273: Input Parameters:
274: + fn - the math function context
275: . alpha - inner scaling (argument)
276: - beta - outer scaling (result)
278: Notes:
279: Given a function $f(x)$ specified by the `FN` type, the scaling parameters can
280: be used to realize the function $\beta f(\alpha x)$. So when these values are given,
281: the procedure for function evaluation will first multiply the argument by $\alpha$,
282: then evaluate the function itself, and finally scale the result by $\beta$.
283: Likewise, these values are also considered when evaluating the derivative.
285: If you want to provide only one of the two scaling factors, set the other
286: one to 1.0.
288: Level: intermediate
290: .seealso: [](sec:fn), `FNGetScale()`, `FNEvaluateFunction()`
291: @*/
292: PetscErrorCode FNSetScale(FN fn,PetscScalar alpha,PetscScalar beta)
293: {
294: PetscFunctionBegin;
298: PetscCheck(PetscAbsScalar(alpha)!=0.0 && PetscAbsScalar(beta)!=0.0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_WRONG,"Scaling factors must be nonzero");
299: fn->alpha = alpha;
300: fn->beta = beta;
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: /*@
305: FNGetScale - Gets the scaling parameters that define the matematical function.
307: Not Collective
309: Input Parameter:
310: . fn - the math function context
312: Output Parameters:
313: + alpha - inner scaling (argument)
314: - beta - outer scaling (result)
316: Level: intermediate
318: .seealso: [](sec:fn), `FNSetScale()`
319: @*/
320: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)
321: {
322: PetscFunctionBegin;
324: if (alpha) *alpha = fn->alpha;
325: if (beta) *beta = fn->beta;
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: /*@
330: FNSetMethod - Selects the method to be used to evaluate functions of matrices.
332: Logically Collective
334: Input Parameters:
335: + fn - the math function context
336: - meth - an index identifying the method
338: Options Database Key:
339: . -fn_method meth - sets the method
341: Notes:
342: In some `FN` types there are more than one algorithm available for computing
343: matrix functions. In that case, this function allows choosing the wanted method.
345: If `meth` is currently set to 0 (the default) and the input argument `A` of
346: `FNEvaluateFunctionMat()` is a symmetric/Hermitian matrix, then the computation
347: is done via the eigendecomposition of `A`, rather than with the general algorithm.
349: Level: intermediate
351: .seealso: [](sec:fn), `FNGetMethod()`, `FNEvaluateFunctionMat()`
352: @*/
353: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)
354: {
355: PetscFunctionBegin;
358: PetscCheck(meth>=0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
359: PetscCheck(meth<=FN_MAX_SOLVE,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
360: fn->method = meth;
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: /*@
365: FNGetMethod - Gets the method currently used in the `FN`.
367: Not Collective
369: Input Parameter:
370: . fn - the math function context
372: Output Parameter:
373: . meth - identifier of the method
375: Level: intermediate
377: .seealso: [](sec:fn), `FNSetMethod()`
378: @*/
379: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)
380: {
381: PetscFunctionBegin;
383: PetscAssertPointer(meth,2);
384: *meth = fn->method;
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /*@
389: FNSetParallel - Selects the mode of operation in parallel runs.
391: Logically Collective
393: Input Parameters:
394: + fn - the math function context
395: - pmode - the parallel mode
397: Options Database Key:
398: . -fn_parallel (redundant|synchronized) - sets the parallel mode
400: Notes:
401: This is relevant only when the function is evaluated on a matrix, with
402: either `FNEvaluateFunctionMat()` or `FNEvaluateFunctionMatVec()`.
404: In the `redundant` parallel mode, all processes will make the computation
405: redundantly, starting from the same data, and producing the same result.
406: This result may be slightly different in the different processes if using a
407: multithreaded BLAS library, which may cause issues in ill-conditioned problems.
409: In the `synchronized` parallel mode, only the first MPI process performs the
410: computation and then the computed matrix is broadcast to the other
411: processes in the communicator. This communication is done automatically at
412: the end of `FNEvaluateFunctionMat()` or `FNEvaluateFunctionMatVec()`.
414: Level: advanced
416: .seealso: [](sec:fn), `FNEvaluateFunctionMat()`, `FNEvaluateFunctionMatVec()`, `FNGetParallel()`
417: @*/
418: PetscErrorCode FNSetParallel(FN fn,FNParallelType pmode)
419: {
420: PetscFunctionBegin;
423: fn->pmode = pmode;
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }
427: /*@
428: FNGetParallel - Gets the mode of operation in parallel runs.
430: Not Collective
432: Input Parameter:
433: . fn - the math function context
435: Output Parameter:
436: . pmode - the parallel mode
438: Level: advanced
440: .seealso: [](sec:fn), `FNSetParallel()`
441: @*/
442: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)
443: {
444: PetscFunctionBegin;
446: PetscAssertPointer(pmode,2);
447: *pmode = fn->pmode;
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: /*@
452: FNEvaluateFunction - Computes the value of the function $f(x)$ for a given $x$.
454: Not Collective
456: Input Parameters:
457: + fn - the math function context
458: - x - the value where the function must be evaluated
460: Output Parameter:
461: . y - the result of $f(x)$
463: Note:
464: Scaling factors are taken into account, so the actual function evaluation
465: will return $\beta f(\alpha x)$.
467: Level: intermediate
469: .seealso: [](sec:fn), `FNEvaluateDerivative()`, `FNEvaluateFunctionMat()`, `FNSetScale()`
470: @*/
471: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)
472: {
473: PetscScalar xf,yf;
475: PetscFunctionBegin;
478: PetscAssertPointer(y,3);
479: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
480: xf = fn->alpha*x;
481: PetscUseTypeMethod(fn,evaluatefunction,xf,&yf);
482: *y = fn->beta*yf;
483: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: /*@
488: FNEvaluateDerivative - Computes the value of the derivative $f'(x)$ for a given $x$.
490: Not Collective
492: Input Parameters:
493: + fn - the math function context
494: - x - the value where the derivative must be evaluated
496: Output Parameter:
497: . y - the result of $f'(x)$
499: Note:
500: Scaling factors are taken into account, so the actual derivative evaluation will
501: return $\alpha\beta f'(\alpha x)$.
503: Level: intermediate
505: .seealso: [](sec:fn), `FNEvaluateFunction()`, `FNSetScale()`
506: @*/
507: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)
508: {
509: PetscScalar xf,yf;
511: PetscFunctionBegin;
514: PetscAssertPointer(y,3);
515: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
516: xf = fn->alpha*x;
517: PetscUseTypeMethod(fn,evaluatederivative,xf,&yf);
518: *y = fn->alpha*fn->beta*yf;
519: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,const PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)
524: {
525: PetscInt i,j;
526: PetscBLASInt n,k,ld,lwork;
527: PetscScalar *Q,*W,*work,adummy,a,x,y,one=1.0,zero=0.0;
528: PetscReal *eig,dummy;
529: #if PetscDefined(USE_COMPLEX)
530: PetscReal *rwork,rdummy;
531: #endif
533: PetscFunctionBegin;
534: PetscCall(PetscBLASIntCast(m,&n));
535: ld = n;
536: k = firstonly? 1: n;
538: /* workspace query and memory allocation */
539: lwork = -1;
540: #if PetscDefined(USE_COMPLEX)
541: PetscCallLAPACKInfo("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&rdummy,&info));
542: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
543: PetscCall(PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork));
544: #else
545: PetscCallLAPACKInfo("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&info));
546: PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
547: PetscCall(PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work));
548: #endif
550: /* compute eigendecomposition */
551: for (j=0;j<n;j++) for (i=j;i<n;i++) Q[i+j*ld] = As[i+j*ld];
552: #if PetscDefined(USE_COMPLEX)
553: PetscCallLAPACKInfo("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
554: #else
555: PetscCallLAPACKInfo("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
556: #endif
558: /* W = f(Lambda)*Q' */
559: for (i=0;i<n;i++) {
560: x = fn->alpha*eig[i];
561: PetscUseTypeMethod(fn,evaluatefunction,x,&y); /* y = f(x) */
562: for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
563: }
564: /* Bs = Q*W */
565: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
566: #if PetscDefined(USE_COMPLEX)
567: PetscCall(PetscFree5(eig,Q,W,work,rwork));
568: #else
569: PetscCall(PetscFree4(eig,Q,W,work));
570: #endif
571: PetscCall(PetscLogFlops(9.0*n*n*n+2.0*n*n*n));
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: /*
576: FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
577: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
578: decomposition of A is A=Q*D*Q'
579: */
580: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)
581: {
582: PetscInt m;
583: const PetscScalar *As;
584: PetscScalar *Bs;
586: PetscFunctionBegin;
587: PetscCall(MatDenseGetArrayRead(A,&As));
588: PetscCall(MatDenseGetArray(B,&Bs));
589: PetscCall(MatGetSize(A,&m,NULL));
590: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE));
591: PetscCall(MatDenseRestoreArrayRead(A,&As));
592: PetscCall(MatDenseRestoreArray(B,&Bs));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: static PetscErrorCode FNEvaluateFunctionMat_Basic(FN fn,Mat A,Mat F)
597: {
598: PetscBool iscuda;
600: PetscFunctionBegin;
601: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
602: 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));
603: if (iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatcuda[fn->method],A,F);
604: else if (fn->ops->evaluatefunctionmat[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmat[fn->method],A,F);
605: else {
606: PetscCheck(fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
607: PetscCheck(!fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
608: }
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
613: {
614: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
615: PetscInt m,n;
616: PetscMPIInt size,rank,n2;
617: PetscScalar *pF;
618: Mat M,F;
620: PetscFunctionBegin;
621: /* destination matrix */
622: F = B?B:A;
624: /* check symmetry of A */
625: PetscCall(MatIsHermitianKnown(A,&set,&flg));
626: symm = set? flg: PETSC_FALSE;
628: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
629: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
630: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
631: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
632: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
633: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
634: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
635: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
636: PetscCall(FNEvaluateFunctionMat_Sym_Default(fn,A,F));
637: } else {
638: /* scale argument */
639: if (fn->alpha!=(PetscScalar)1.0) {
640: PetscCall(FN_AllocateWorkMat(fn,A,&M));
641: PetscCall(MatScale(M,fn->alpha));
642: } else M = A;
643: PetscCall(FNEvaluateFunctionMat_Basic(fn,M,F));
644: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
645: /* scale result */
646: PetscCall(MatScale(F,fn->beta));
647: }
648: PetscCall(PetscFPTrapPop());
649: }
650: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { /* synchronize */
651: PetscCall(MatGetSize(A,&m,&n));
652: PetscCall(MatDenseGetArray(F,&pF));
653: PetscCall(PetscMPIIntCast(n*n,&n2));
654: PetscCallMPI(MPI_Bcast(pF,n2,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
655: PetscCall(MatDenseRestoreArray(F,&pF));
656: }
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: /*@
661: FNEvaluateFunctionMat - Computes the value of the function $f(A)$ for a given
662: matrix $A$, where the result is also a matrix.
664: Logically Collective
666: Input Parameters:
667: + fn - the math function context
668: - A - matrix on which the function must be evaluated
670: Output Parameter:
671: . B - (optional) matrix resulting from evaluating $f(A)$
673: Notes:
674: Matrix `A` must be a square sequential dense `Mat`, with all entries equal on
675: all processes (otherwise each process will compute different results).
676: If matrix `B` is provided, it must also be a square sequential dense `Mat`, and
677: both matrices must have the same dimensions. If `B` is `NULL` (or `B`=`A`) then
678: the function will perform an in-place computation, overwriting `A` with $f(A)$.
680: If `A` is known to be real symmetric or complex Hermitian then it is
681: recommended to set the appropriate flag with `MatSetOption()`, because
682: symmetry can sometimes be exploited by the algorithm.
684: Scaling factors are taken into account, so the actual function evaluation
685: will return $\beta f(\alpha A)$.
687: Level: advanced
689: .seealso: [](sec:fn), `FNEvaluateFunction()`, `FNEvaluateFunctionMatVec()`, `FNSetMethod()`
690: @*/
691: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
692: {
693: PetscBool inplace=PETSC_FALSE;
694: PetscInt m,n,n1;
695: MatType type;
697: PetscFunctionBegin;
702: if (B) {
705: } else inplace = PETSC_TRUE;
706: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA); //SlepcMatCheckSeq(A);
707: PetscCall(MatGetSize(A,&m,&n));
708: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
709: if (!inplace) {
710: PetscCall(MatGetType(A,&type));
711: PetscCheckTypeName(B,type);
712: n1 = n;
713: PetscCall(MatGetSize(B,&m,&n));
714: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat B is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
715: PetscCheck(n1==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrices A and B must have the same dimension");
716: }
718: /* evaluate matrix function */
719: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
720: PetscCall(FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE));
721: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: /*
726: FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
727: and then copies the first column.
728: */
729: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
730: {
731: Mat F;
733: PetscFunctionBegin;
734: PetscCall(FN_AllocateWorkMat(fn,A,&F));
735: PetscCall(FNEvaluateFunctionMat_Basic(fn,A,F));
736: PetscCall(MatGetColumnVector(F,v,0));
737: PetscCall(FN_FreeWorkMat(fn,&F));
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: /*
742: FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
743: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
744: decomposition of A is A=Q*D*Q'. Only the first column is computed.
745: */
746: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
747: {
748: PetscInt m;
749: const PetscScalar *As;
750: PetscScalar *vs;
752: PetscFunctionBegin;
753: PetscCall(MatDenseGetArrayRead(A,&As));
754: PetscCall(VecGetArray(v,&vs));
755: PetscCall(MatGetSize(A,&m,NULL));
756: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE));
757: PetscCall(MatDenseRestoreArrayRead(A,&As));
758: PetscCall(VecRestoreArray(v,&vs));
759: PetscFunctionReturn(PETSC_SUCCESS);
760: }
762: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
763: {
764: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
765: PetscInt m,n;
766: Mat M;
767: PetscMPIInt size,rank,n_;
768: PetscScalar *pv;
770: PetscFunctionBegin;
771: /* check symmetry of A */
772: PetscCall(MatIsHermitianKnown(A,&set,&flg));
773: symm = set? flg: PETSC_FALSE;
775: /* evaluate matrix function */
776: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
777: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
778: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
779: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
780: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
781: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
782: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
783: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
784: PetscCall(FNEvaluateFunctionMatVec_Sym_Default(fn,A,v));
785: } else {
786: /* scale argument */
787: if (fn->alpha!=(PetscScalar)1.0) {
788: PetscCall(FN_AllocateWorkMat(fn,A,&M));
789: PetscCall(MatScale(M,fn->alpha));
790: } else M = A;
791: if (iscuda && fn->ops->evaluatefunctionmatveccuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatveccuda[fn->method],M,v);
792: else if (fn->ops->evaluatefunctionmatvec[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatvec[fn->method],M,v);
793: else PetscCall(FNEvaluateFunctionMatVec_Default(fn,M,v));
794: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
795: /* scale result */
796: PetscCall(VecScale(v,fn->beta));
797: }
798: PetscCall(PetscFPTrapPop());
799: }
801: /* synchronize */
802: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
803: PetscCall(MatGetSize(A,&m,&n));
804: PetscCall(VecGetArray(v,&pv));
805: PetscCall(PetscMPIIntCast(n,&n_));
806: PetscCallMPI(MPI_Bcast(pv,n_,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
807: PetscCall(VecRestoreArray(v,&pv));
808: }
809: PetscFunctionReturn(PETSC_SUCCESS);
810: }
812: /*@
813: FNEvaluateFunctionMatVec - Computes the first column of the matrix $f(A)$
814: for a given matrix $A$.
816: Logically Collective
818: Input Parameters:
819: + fn - the math function context
820: - A - matrix on which the function must be evaluated
822: Output Parameter:
823: . v - vector to hold the first column of $f(A)$
825: Notes:
826: This operation is similar to `FNEvaluateFunctionMat()` but returns only
827: the first column of $f(A)$, hence saving computations in most cases.
829: Level: advanced
831: .seealso: [](sec:fn), `FNEvaluateFunction()`, `FNEvaluateFunctionMat()`, `FNSetMethod()`
832: @*/
833: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
834: {
835: PetscInt m,n;
836: PetscBool iscuda;
838: PetscFunctionBegin;
845: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA); //SlepcMatCheckSeq(A);
846: PetscCall(MatGetSize(A,&m,&n));
847: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
848: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
849: PetscCheckTypeName(v,iscuda?VECSEQCUDA:VECSEQ);
850: PetscCall(VecGetSize(v,&m));
851: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrix A and vector v must have the same size");
852: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
853: PetscCall(FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE));
854: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: /*@
859: FNSetFromOptions - Sets `FN` options from the options database.
861: Collective
863: Input Parameter:
864: . fn - the math function context
866: Note:
867: To see all options, run your program with the `-help` option.
869: Level: beginner
871: .seealso: [](sec:fn), `FNSetOptionsPrefix()`
872: @*/
873: PetscErrorCode FNSetFromOptions(FN fn)
874: {
875: char type[256];
876: PetscScalar array[2];
877: PetscInt k,meth;
878: PetscBool flg;
879: FNParallelType pmode;
881: PetscFunctionBegin;
883: PetscCall(FNRegisterAll());
884: PetscObjectOptionsBegin((PetscObject)fn);
885: PetscCall(PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg));
886: if (flg) PetscCall(FNSetType(fn,type));
887: else if (!((PetscObject)fn)->type_name) PetscCall(FNSetType(fn,FNRATIONAL));
889: k = 2;
890: array[0] = 0.0; array[1] = 0.0;
891: PetscCall(PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg));
892: if (flg) {
893: if (k<2) array[1] = 1.0;
894: PetscCall(FNSetScale(fn,array[0],array[1]));
895: }
897: PetscCall(PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg));
898: if (flg) PetscCall(FNSetMethod(fn,meth));
900: PetscCall(PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg));
901: if (flg) PetscCall(FNSetParallel(fn,pmode));
903: PetscTryTypeMethod(fn,setfromoptions,PetscOptionsObject);
904: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fn,PetscOptionsObject));
905: PetscOptionsEnd();
906: PetscFunctionReturn(PETSC_SUCCESS);
907: }
909: /*@
910: FNView - Prints the `FN` data structure.
912: Collective
914: Input Parameters:
915: + fn - the math function context
916: - viewer - optional visualization context
918: Note:
919: The available visualization contexts include
920: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
921: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
922: first process opens the file; all other processes send their data to the
923: first one to print
925: The user can open an alternative visualization context with `PetscViewerASCIIOpen()`
926: to output to a specified file.
928: Use `FNViewFromOptions()` to allow the user to select many different `PetscViewerType`
929: and formats from the options database.
931: Level: beginner
933: .seealso: [](sec:fn), `FNCreate()`, `FNViewFromOptions()`
934: @*/
935: PetscErrorCode FNView(FN fn,PetscViewer viewer)
936: {
937: PetscBool isascii;
938: PetscMPIInt size;
940: PetscFunctionBegin;
942: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer));
944: PetscCheckSameComm(fn,1,viewer,2);
945: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
946: if (isascii) {
947: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer));
948: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
949: if (size>1) PetscCall(PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",FNParallelTypes[fn->pmode]));
950: PetscCall(PetscViewerASCIIPushTab(viewer));
951: PetscTryTypeMethod(fn,view,viewer);
952: PetscCall(PetscViewerASCIIPopTab(viewer));
953: }
954: PetscFunctionReturn(PETSC_SUCCESS);
955: }
957: /*@
958: FNViewFromOptions - View (print) an `FN` object based on values in the options database.
960: Collective
962: Input Parameters:
963: + fn - the math function context
964: . obj - optional object that provides the options prefix used to query the options database
965: - name - command line option
967: Level: intermediate
969: .seealso: [](sec:fn), `FNView()`, `FNCreate()`, `PetscObjectViewFromOptions()`
970: @*/
971: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])
972: {
973: PetscFunctionBegin;
975: PetscCall(PetscObjectViewFromOptions((PetscObject)fn,obj,name));
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: /*@
980: FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
981: different communicator.
983: Collective
985: Input Parameters:
986: + fn - the math function context
987: - comm - MPI communicator
989: Output Parameter:
990: . newfn - location to put the new `FN` context
992: Note:
993: In order to use the same MPI communicator as in the original object,
994: use `PetscObjectComm`((`PetscObject`)`fn`).
996: Level: developer
998: .seealso: [](sec:fn), `FNCreate()`
999: @*/
1000: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)
1001: {
1002: FNType type;
1003: PetscScalar alpha,beta;
1004: PetscInt meth;
1005: FNParallelType ptype;
1007: PetscFunctionBegin;
1010: PetscAssertPointer(newfn,3);
1011: PetscCall(FNCreate(comm,newfn));
1012: PetscCall(FNGetType(fn,&type));
1013: PetscCall(FNSetType(*newfn,type));
1014: PetscCall(FNGetScale(fn,&alpha,&beta));
1015: PetscCall(FNSetScale(*newfn,alpha,beta));
1016: PetscCall(FNGetMethod(fn,&meth));
1017: PetscCall(FNSetMethod(*newfn,meth));
1018: PetscCall(FNGetParallel(fn,&ptype));
1019: PetscCall(FNSetParallel(*newfn,ptype));
1020: PetscTryTypeMethod(fn,duplicate,comm,newfn);
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: /*@
1025: FNDestroy - Destroys an `FN` context that was created with `FNCreate()`.
1027: Collective
1029: Input Parameter:
1030: . fn - the math function context
1032: Level: beginner
1034: .seealso: [](sec:fn), `FNCreate()`
1035: @*/
1036: PetscErrorCode FNDestroy(FN *fn)
1037: {
1038: PetscInt i;
1040: PetscFunctionBegin;
1041: if (!*fn) PetscFunctionReturn(PETSC_SUCCESS);
1043: if (--((PetscObject)*fn)->refct > 0) { *fn = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
1044: PetscTryTypeMethod(*fn,destroy);
1045: for (i=0;i<(*fn)->nw;i++) PetscCall(MatDestroy(&(*fn)->W[i]));
1046: PetscCall(PetscHeaderDestroy(fn));
1047: PetscFunctionReturn(PETSC_SUCCESS);
1048: }
1050: /*@
1051: FNRegister - Adds a mathematical function to the `FN` package.
1053: Not Collective
1055: Input Parameters:
1056: + name - name of a new user-defined `FN`
1057: - function - routine to create the context
1059: Notes:
1060: `FNRegister()` may be called multiple times to add several user-defined functions.
1062: Level: advanced
1064: .seealso: [](sec:fn), `FNRegisterAll()`
1065: @*/
1066: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1067: {
1068: PetscFunctionBegin;
1069: PetscCall(FNInitializePackage());
1070: PetscCall(PetscFunctionListAdd(&FNList,name,function));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }