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: /*@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: [](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: /*@C
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,info;
527: PetscScalar *Q,*W,*work,adummy,a,x,y,one=1.0,zero=0.0;
528: PetscReal *eig,dummy;
529: #if defined(PETSC_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 defined(PETSC_USE_COMPLEX)
541: PetscCallBLAS("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: PetscCallBLAS("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 defined(PETSC_USE_COMPLEX)
553: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
554: #else
555: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
556: #endif
557: SlepcCheckLapackInfo("syev",info);
559: /* W = f(Lambda)*Q' */
560: for (i=0;i<n;i++) {
561: x = fn->alpha*eig[i];
562: PetscUseTypeMethod(fn,evaluatefunction,x,&y); /* y = f(x) */
563: for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
564: }
565: /* Bs = Q*W */
566: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
567: #if defined(PETSC_USE_COMPLEX)
568: PetscCall(PetscFree5(eig,Q,W,work,rwork));
569: #else
570: PetscCall(PetscFree4(eig,Q,W,work));
571: #endif
572: PetscCall(PetscLogFlops(9.0*n*n*n+2.0*n*n*n));
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /*
577: FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
578: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
579: decomposition of A is A=Q*D*Q'
580: */
581: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)
582: {
583: PetscInt m;
584: const PetscScalar *As;
585: PetscScalar *Bs;
587: PetscFunctionBegin;
588: PetscCall(MatDenseGetArrayRead(A,&As));
589: PetscCall(MatDenseGetArray(B,&Bs));
590: PetscCall(MatGetSize(A,&m,NULL));
591: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE));
592: PetscCall(MatDenseRestoreArrayRead(A,&As));
593: PetscCall(MatDenseRestoreArray(B,&Bs));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: static PetscErrorCode FNEvaluateFunctionMat_Basic(FN fn,Mat A,Mat F)
598: {
599: PetscBool iscuda;
601: PetscFunctionBegin;
602: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
603: 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));
604: if (iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatcuda[fn->method],A,F);
605: else if (fn->ops->evaluatefunctionmat[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmat[fn->method],A,F);
606: else {
607: PetscCheck(fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
608: PetscCheck(!fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
609: }
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
614: {
615: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
616: PetscInt m,n;
617: PetscMPIInt size,rank,n2;
618: PetscScalar *pF;
619: Mat M,F;
621: PetscFunctionBegin;
622: /* destination matrix */
623: F = B?B:A;
625: /* check symmetry of A */
626: PetscCall(MatIsHermitianKnown(A,&set,&flg));
627: symm = set? flg: PETSC_FALSE;
629: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
630: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
631: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
632: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
633: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
634: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
635: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
636: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
637: PetscCall(FNEvaluateFunctionMat_Sym_Default(fn,A,F));
638: } else {
639: /* scale argument */
640: if (fn->alpha!=(PetscScalar)1.0) {
641: PetscCall(FN_AllocateWorkMat(fn,A,&M));
642: PetscCall(MatScale(M,fn->alpha));
643: } else M = A;
644: PetscCall(FNEvaluateFunctionMat_Basic(fn,M,F));
645: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
646: /* scale result */
647: PetscCall(MatScale(F,fn->beta));
648: }
649: PetscCall(PetscFPTrapPop());
650: }
651: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { /* synchronize */
652: PetscCall(MatGetSize(A,&m,&n));
653: PetscCall(MatDenseGetArray(F,&pF));
654: PetscCall(PetscMPIIntCast(n*n,&n2));
655: PetscCallMPI(MPI_Bcast(pF,n2,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
656: PetscCall(MatDenseRestoreArray(F,&pF));
657: }
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: /*@
662: FNEvaluateFunctionMat - Computes the value of the function $f(A)$ for a given
663: matrix $A$, where the result is also a matrix.
665: Logically Collective
667: Input Parameters:
668: + fn - the math function context
669: - A - matrix on which the function must be evaluated
671: Output Parameter:
672: . B - (optional) matrix resulting from evaluating $f(A)$
674: Notes:
675: Matrix `A` must be a square sequential dense `Mat`, with all entries equal on
676: all processes (otherwise each process will compute different results).
677: If matrix `B` is provided, it must also be a square sequential dense `Mat`, and
678: both matrices must have the same dimensions. If `B` is `NULL` (or `B`=`A`) then
679: the function will perform an in-place computation, overwriting `A` with $f(A)$.
681: If `A` is known to be real symmetric or complex Hermitian then it is
682: recommended to set the appropriate flag with `MatSetOption()`, because
683: symmetry can sometimes be exploited by the algorithm.
685: Scaling factors are taken into account, so the actual function evaluation
686: will return $\beta f(\alpha A)$.
688: Level: advanced
690: .seealso: [](sec:fn), `FNEvaluateFunction()`, `FNEvaluateFunctionMatVec()`, `FNSetMethod()`
691: @*/
692: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
693: {
694: PetscBool inplace=PETSC_FALSE;
695: PetscInt m,n,n1;
696: MatType type;
698: PetscFunctionBegin;
703: if (B) {
706: } else inplace = PETSC_TRUE;
707: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA); //SlepcMatCheckSeq(A);
708: PetscCall(MatGetSize(A,&m,&n));
709: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
710: if (!inplace) {
711: PetscCall(MatGetType(A,&type));
712: PetscCheckTypeName(B,type);
713: n1 = n;
714: PetscCall(MatGetSize(B,&m,&n));
715: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat B is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
716: PetscCheck(n1==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrices A and B must have the same dimension");
717: }
719: /* evaluate matrix function */
720: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
721: PetscCall(FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE));
722: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: /*
727: FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
728: and then copies the first column.
729: */
730: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
731: {
732: Mat F;
734: PetscFunctionBegin;
735: PetscCall(FN_AllocateWorkMat(fn,A,&F));
736: PetscCall(FNEvaluateFunctionMat_Basic(fn,A,F));
737: PetscCall(MatGetColumnVector(F,v,0));
738: PetscCall(FN_FreeWorkMat(fn,&F));
739: PetscFunctionReturn(PETSC_SUCCESS);
740: }
742: /*
743: FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
744: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
745: decomposition of A is A=Q*D*Q'. Only the first column is computed.
746: */
747: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
748: {
749: PetscInt m;
750: const PetscScalar *As;
751: PetscScalar *vs;
753: PetscFunctionBegin;
754: PetscCall(MatDenseGetArrayRead(A,&As));
755: PetscCall(VecGetArray(v,&vs));
756: PetscCall(MatGetSize(A,&m,NULL));
757: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE));
758: PetscCall(MatDenseRestoreArrayRead(A,&As));
759: PetscCall(VecRestoreArray(v,&vs));
760: PetscFunctionReturn(PETSC_SUCCESS);
761: }
763: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
764: {
765: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
766: PetscInt m,n;
767: Mat M;
768: PetscMPIInt size,rank,n_;
769: PetscScalar *pv;
771: PetscFunctionBegin;
772: /* check symmetry of A */
773: PetscCall(MatIsHermitianKnown(A,&set,&flg));
774: symm = set? flg: PETSC_FALSE;
776: /* evaluate matrix function */
777: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
778: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
779: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
780: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
781: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
782: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
783: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
784: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
785: PetscCall(FNEvaluateFunctionMatVec_Sym_Default(fn,A,v));
786: } else {
787: /* scale argument */
788: if (fn->alpha!=(PetscScalar)1.0) {
789: PetscCall(FN_AllocateWorkMat(fn,A,&M));
790: PetscCall(MatScale(M,fn->alpha));
791: } else M = A;
792: if (iscuda && fn->ops->evaluatefunctionmatveccuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatveccuda[fn->method],M,v);
793: else if (fn->ops->evaluatefunctionmatvec[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatvec[fn->method],M,v);
794: else PetscCall(FNEvaluateFunctionMatVec_Default(fn,M,v));
795: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
796: /* scale result */
797: PetscCall(VecScale(v,fn->beta));
798: }
799: PetscCall(PetscFPTrapPop());
800: }
802: /* synchronize */
803: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
804: PetscCall(MatGetSize(A,&m,&n));
805: PetscCall(VecGetArray(v,&pv));
806: PetscCall(PetscMPIIntCast(n,&n_));
807: PetscCallMPI(MPI_Bcast(pv,n_,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
808: PetscCall(VecRestoreArray(v,&pv));
809: }
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
813: /*@
814: FNEvaluateFunctionMatVec - Computes the first column of the matrix $f(A)$
815: for a given matrix $A$.
817: Logically Collective
819: Input Parameters:
820: + fn - the math function context
821: - A - matrix on which the function must be evaluated
823: Output Parameter:
824: . v - vector to hold the first column of $f(A)$
826: Notes:
827: This operation is similar to `FNEvaluateFunctionMat()` but returns only
828: the first column of $f(A)$, hence saving computations in most cases.
830: Level: advanced
832: .seealso: [](sec:fn), `FNEvaluateFunction()`, `FNEvaluateFunctionMat()`, `FNSetMethod()`
833: @*/
834: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
835: {
836: PetscInt m,n;
837: PetscBool iscuda;
839: PetscFunctionBegin;
846: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA); //SlepcMatCheckSeq(A);
847: PetscCall(MatGetSize(A,&m,&n));
848: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
849: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
850: PetscCheckTypeName(v,iscuda?VECSEQCUDA:VECSEQ);
851: PetscCall(VecGetSize(v,&m));
852: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrix A and vector v must have the same size");
853: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
854: PetscCall(FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE));
855: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
859: /*@
860: FNSetFromOptions - Sets `FN` options from the options database.
862: Collective
864: Input Parameter:
865: . fn - the math function context
867: Note:
868: To see all options, run your program with the `-help` option.
870: Level: beginner
872: .seealso: [](sec:fn), `FNSetOptionsPrefix()`
873: @*/
874: PetscErrorCode FNSetFromOptions(FN fn)
875: {
876: char type[256];
877: PetscScalar array[2];
878: PetscInt k,meth;
879: PetscBool flg;
880: FNParallelType pmode;
882: PetscFunctionBegin;
884: PetscCall(FNRegisterAll());
885: PetscObjectOptionsBegin((PetscObject)fn);
886: PetscCall(PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg));
887: if (flg) PetscCall(FNSetType(fn,type));
888: else if (!((PetscObject)fn)->type_name) PetscCall(FNSetType(fn,FNRATIONAL));
890: k = 2;
891: array[0] = 0.0; array[1] = 0.0;
892: PetscCall(PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg));
893: if (flg) {
894: if (k<2) array[1] = 1.0;
895: PetscCall(FNSetScale(fn,array[0],array[1]));
896: }
898: PetscCall(PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg));
899: if (flg) PetscCall(FNSetMethod(fn,meth));
901: PetscCall(PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg));
902: if (flg) PetscCall(FNSetParallel(fn,pmode));
904: PetscTryTypeMethod(fn,setfromoptions,PetscOptionsObject);
905: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fn,PetscOptionsObject));
906: PetscOptionsEnd();
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: /*@
911: FNView - Prints the `FN` data structure.
913: Collective
915: Input Parameters:
916: + fn - the math function context
917: - viewer - optional visualization context
919: Note:
920: The available visualization contexts include
921: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
922: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
923: first process opens the file; all other processes send their data to the
924: first one to print
926: The user can open an alternative visualization context with `PetscViewerASCIIOpen()`
927: to output to a specified file.
929: Use `FNViewFromOptions()` to allow the user to select many different `PetscViewerType`
930: and formats from the options database.
932: Level: beginner
934: .seealso: [](sec:fn), `FNCreate()`, `FNViewFromOptions()`
935: @*/
936: PetscErrorCode FNView(FN fn,PetscViewer viewer)
937: {
938: PetscBool isascii;
939: PetscMPIInt size;
941: PetscFunctionBegin;
943: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer));
945: PetscCheckSameComm(fn,1,viewer,2);
946: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
947: if (isascii) {
948: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer));
949: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
950: if (size>1) PetscCall(PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",FNParallelTypes[fn->pmode]));
951: PetscCall(PetscViewerASCIIPushTab(viewer));
952: PetscTryTypeMethod(fn,view,viewer);
953: PetscCall(PetscViewerASCIIPopTab(viewer));
954: }
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: /*@
959: FNViewFromOptions - View (print) an `FN` object based on values in the options database.
961: Collective
963: Input Parameters:
964: + fn - the math function context
965: . obj - optional object that provides the options prefix used to query the options database
966: - name - command line option
968: Level: intermediate
970: .seealso: [](sec:fn), `FNView()`, `FNCreate()`, `PetscObjectViewFromOptions()`
971: @*/
972: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])
973: {
974: PetscFunctionBegin;
976: PetscCall(PetscObjectViewFromOptions((PetscObject)fn,obj,name));
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: /*@
981: FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
982: different communicator.
984: Collective
986: Input Parameters:
987: + fn - the math function context
988: - comm - MPI communicator
990: Output Parameter:
991: . newfn - location to put the new `FN` context
993: Note:
994: In order to use the same MPI communicator as in the original object,
995: use `PetscObjectComm`((`PetscObject`)`fn`).
997: Level: developer
999: .seealso: [](sec:fn), `FNCreate()`
1000: @*/
1001: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)
1002: {
1003: FNType type;
1004: PetscScalar alpha,beta;
1005: PetscInt meth;
1006: FNParallelType ptype;
1008: PetscFunctionBegin;
1011: PetscAssertPointer(newfn,3);
1012: PetscCall(FNCreate(comm,newfn));
1013: PetscCall(FNGetType(fn,&type));
1014: PetscCall(FNSetType(*newfn,type));
1015: PetscCall(FNGetScale(fn,&alpha,&beta));
1016: PetscCall(FNSetScale(*newfn,alpha,beta));
1017: PetscCall(FNGetMethod(fn,&meth));
1018: PetscCall(FNSetMethod(*newfn,meth));
1019: PetscCall(FNGetParallel(fn,&ptype));
1020: PetscCall(FNSetParallel(*newfn,ptype));
1021: PetscTryTypeMethod(fn,duplicate,comm,newfn);
1022: PetscFunctionReturn(PETSC_SUCCESS);
1023: }
1025: /*@
1026: FNDestroy - Destroys an `FN` context that was created with `FNCreate()`.
1028: Collective
1030: Input Parameter:
1031: . fn - the math function context
1033: Level: beginner
1035: .seealso: [](sec:fn), `FNCreate()`
1036: @*/
1037: PetscErrorCode FNDestroy(FN *fn)
1038: {
1039: PetscInt i;
1041: PetscFunctionBegin;
1042: if (!*fn) PetscFunctionReturn(PETSC_SUCCESS);
1044: if (--((PetscObject)*fn)->refct > 0) { *fn = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
1045: PetscTryTypeMethod(*fn,destroy);
1046: for (i=0;i<(*fn)->nw;i++) PetscCall(MatDestroy(&(*fn)->W[i]));
1047: PetscCall(PetscHeaderDestroy(fn));
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: /*@C
1052: FNRegister - Adds a mathematical function to the `FN` package.
1054: Not Collective
1056: Input Parameters:
1057: + name - name of a new user-defined `FN`
1058: - function - routine to create the context
1060: Notes:
1061: `FNRegister()` may be called multiple times to add several user-defined functions.
1063: Level: advanced
1065: .seealso: [](sec:fn), `FNRegisterAll()`
1066: @*/
1067: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1068: {
1069: PetscFunctionBegin;
1070: PetscCall(FNInitializePackage());
1071: PetscCall(PetscFunctionListAdd(&FNList,name,function));
1072: PetscFunctionReturn(PETSC_SUCCESS);
1073: }