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 `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: Level: intermediate
253: .seealso: [](sec:fn), `FNSetType()`
254: @*/
255: PetscErrorCode FNGetType(FN fn,FNType *type)
256: {
257: PetscFunctionBegin;
259: PetscAssertPointer(type,2);
260: *type = ((PetscObject)fn)->type_name;
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: /*@
265: FNSetScale - Sets the scaling parameters that define the matematical function.
267: Logically Collective
269: Input Parameters:
270: + fn - the math function context
271: . alpha - inner scaling (argument)
272: - beta - outer scaling (result)
274: Notes:
275: Given a function $f(x)$ specified by the `FN` type, the scaling parameters can
276: be used to realize the function $\beta f(\alpha x)$. So when these values are given,
277: the procedure for function evaluation will first multiply the argument by $\alpha$,
278: then evaluate the function itself, and finally scale the result by $\beta$.
279: Likewise, these values are also considered when evaluating the derivative.
281: If you want to provide only one of the two scaling factors, set the other
282: one to 1.0.
284: Level: intermediate
286: .seealso: [](sec:fn), `FNGetScale()`, `FNEvaluateFunction()`
287: @*/
288: PetscErrorCode FNSetScale(FN fn,PetscScalar alpha,PetscScalar beta)
289: {
290: PetscFunctionBegin;
294: PetscCheck(PetscAbsScalar(alpha)!=0.0 && PetscAbsScalar(beta)!=0.0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_WRONG,"Scaling factors must be nonzero");
295: fn->alpha = alpha;
296: fn->beta = beta;
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /*@
301: FNGetScale - Gets the scaling parameters that define the matematical function.
303: Not Collective
305: Input Parameter:
306: . fn - the math function context
308: Output Parameters:
309: + alpha - inner scaling (argument)
310: - beta - outer scaling (result)
312: Level: intermediate
314: .seealso: [](sec:fn), `FNSetScale()`
315: @*/
316: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)
317: {
318: PetscFunctionBegin;
320: if (alpha) *alpha = fn->alpha;
321: if (beta) *beta = fn->beta;
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*@
326: FNSetMethod - Selects the method to be used to evaluate functions of matrices.
328: Logically Collective
330: Input Parameters:
331: + fn - the math function context
332: - meth - an index identifying the method
334: Options Database Key:
335: . -fn_method <meth> - Sets the method
337: Notes:
338: In some `FN` types there are more than one algorithm available for computing
339: matrix functions. In that case, this function allows choosing the wanted method.
341: If `meth` is currently set to 0 (the default) and the input argument `A` of
342: `FNEvaluateFunctionMat()` is a symmetric/Hermitian matrix, then the computation
343: is done via the eigendecomposition of `A`, rather than with the general algorithm.
345: Level: intermediate
347: .seealso: [](sec:fn), `FNGetMethod()`, `FNEvaluateFunctionMat()`
348: @*/
349: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)
350: {
351: PetscFunctionBegin;
354: PetscCheck(meth>=0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
355: PetscCheck(meth<=FN_MAX_SOLVE,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
356: fn->method = meth;
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: /*@
361: FNGetMethod - Gets the method currently used in the `FN`.
363: Not Collective
365: Input Parameter:
366: . fn - the math function context
368: Output Parameter:
369: . meth - identifier of the method
371: Level: intermediate
373: .seealso: [](sec:fn), `FNSetMethod()`
374: @*/
375: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)
376: {
377: PetscFunctionBegin;
379: PetscAssertPointer(meth,2);
380: *meth = fn->method;
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*@
385: FNSetParallel - Selects the mode of operation in parallel runs.
387: Logically Collective
389: Input Parameters:
390: + fn - the math function context
391: - pmode - the parallel mode
393: Options Database Key:
394: . -fn_parallel <mode> - Sets the parallel mode, either `redundant` or `synchronized`
396: Notes:
397: This is relevant only when the function is evaluated on a matrix, with
398: either `FNEvaluateFunctionMat()` or `FNEvaluateFunctionMatVec()`.
400: In the `redundant` parallel mode, all processes will make the computation
401: redundantly, starting from the same data, and producing the same result.
402: This result may be slightly different in the different processes if using a
403: multithreaded BLAS library, which may cause issues in ill-conditioned problems.
405: In the `synchronized` parallel mode, only the first MPI process performs the
406: computation and then the computed matrix is broadcast to the other
407: processes in the communicator. This communication is done automatically at
408: the end of `FNEvaluateFunctionMat()` or `FNEvaluateFunctionMatVec()`.
410: Level: advanced
412: .seealso: [](sec:fn), `FNEvaluateFunctionMat()` or `FNEvaluateFunctionMatVec()`, `FNGetParallel()`
413: @*/
414: PetscErrorCode FNSetParallel(FN fn,FNParallelType pmode)
415: {
416: PetscFunctionBegin;
419: fn->pmode = pmode;
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: /*@
424: FNGetParallel - Gets the mode of operation in parallel runs.
426: Not Collective
428: Input Parameter:
429: . fn - the math function context
431: Output Parameter:
432: . pmode - the parallel mode
434: Level: advanced
436: .seealso: [](sec:fn), `FNSetParallel()`
437: @*/
438: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)
439: {
440: PetscFunctionBegin;
442: PetscAssertPointer(pmode,2);
443: *pmode = fn->pmode;
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: /*@
448: FNEvaluateFunction - Computes the value of the function $f(x)$ for a given $x$.
450: Not Collective
452: Input Parameters:
453: + fn - the math function context
454: - x - the value where the function must be evaluated
456: Output Parameter:
457: . y - the result of $f(x)$
459: Note:
460: Scaling factors are taken into account, so the actual function evaluation
461: will return $\beta f(\alpha x)$.
463: Level: intermediate
465: .seealso: [](sec:fn), `FNEvaluateDerivative()`, `FNEvaluateFunctionMat()`, `FNSetScale()`
466: @*/
467: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)
468: {
469: PetscScalar xf,yf;
471: PetscFunctionBegin;
474: PetscAssertPointer(y,3);
475: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
476: xf = fn->alpha*x;
477: PetscUseTypeMethod(fn,evaluatefunction,xf,&yf);
478: *y = fn->beta*yf;
479: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: /*@
484: FNEvaluateDerivative - Computes the value of the derivative $f'(x)$ for a given $x$.
486: Not Collective
488: Input Parameters:
489: + fn - the math function context
490: - x - the value where the derivative must be evaluated
492: Output Parameter:
493: . y - the result of $f'(x)$
495: Note:
496: Scaling factors are taken into account, so the actual derivative evaluation will
497: return $\alpha\beta f'(\alpha x)$.
499: Level: intermediate
501: .seealso: [](sec:fn), `FNEvaluateFunction()`, `FNSetScale()`
502: @*/
503: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)
504: {
505: PetscScalar xf,yf;
507: PetscFunctionBegin;
510: PetscAssertPointer(y,3);
511: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
512: xf = fn->alpha*x;
513: PetscUseTypeMethod(fn,evaluatederivative,xf,&yf);
514: *y = fn->alpha*fn->beta*yf;
515: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,const PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)
520: {
521: PetscInt i,j;
522: PetscBLASInt n,k,ld,lwork,info;
523: PetscScalar *Q,*W,*work,adummy,a,x,y,one=1.0,zero=0.0;
524: PetscReal *eig,dummy;
525: #if defined(PETSC_USE_COMPLEX)
526: PetscReal *rwork,rdummy;
527: #endif
529: PetscFunctionBegin;
530: PetscCall(PetscBLASIntCast(m,&n));
531: ld = n;
532: k = firstonly? 1: n;
534: /* workspace query and memory allocation */
535: lwork = -1;
536: #if defined(PETSC_USE_COMPLEX)
537: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&rdummy,&info));
538: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
539: PetscCall(PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork));
540: #else
541: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&info));
542: PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
543: PetscCall(PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work));
544: #endif
546: /* compute eigendecomposition */
547: for (j=0;j<n;j++) for (i=j;i<n;i++) Q[i+j*ld] = As[i+j*ld];
548: #if defined(PETSC_USE_COMPLEX)
549: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
550: #else
551: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
552: #endif
553: SlepcCheckLapackInfo("syev",info);
555: /* W = f(Lambda)*Q' */
556: for (i=0;i<n;i++) {
557: x = fn->alpha*eig[i];
558: PetscUseTypeMethod(fn,evaluatefunction,x,&y); /* y = f(x) */
559: for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
560: }
561: /* Bs = Q*W */
562: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
563: #if defined(PETSC_USE_COMPLEX)
564: PetscCall(PetscFree5(eig,Q,W,work,rwork));
565: #else
566: PetscCall(PetscFree4(eig,Q,W,work));
567: #endif
568: PetscCall(PetscLogFlops(9.0*n*n*n+2.0*n*n*n));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: /*
573: FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
574: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
575: decomposition of A is A=Q*D*Q'
576: */
577: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)
578: {
579: PetscInt m;
580: const PetscScalar *As;
581: PetscScalar *Bs;
583: PetscFunctionBegin;
584: PetscCall(MatDenseGetArrayRead(A,&As));
585: PetscCall(MatDenseGetArray(B,&Bs));
586: PetscCall(MatGetSize(A,&m,NULL));
587: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE));
588: PetscCall(MatDenseRestoreArrayRead(A,&As));
589: PetscCall(MatDenseRestoreArray(B,&Bs));
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: static PetscErrorCode FNEvaluateFunctionMat_Basic(FN fn,Mat A,Mat F)
594: {
595: PetscBool iscuda;
597: PetscFunctionBegin;
598: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
599: 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));
600: if (iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatcuda[fn->method],A,F);
601: else if (fn->ops->evaluatefunctionmat[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmat[fn->method],A,F);
602: else {
603: PetscCheck(fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
604: PetscCheck(!fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
605: }
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
610: {
611: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
612: PetscInt m,n;
613: PetscMPIInt size,rank,n2;
614: PetscScalar *pF;
615: Mat M,F;
617: PetscFunctionBegin;
618: /* destination matrix */
619: F = B?B:A;
621: /* check symmetry of A */
622: PetscCall(MatIsHermitianKnown(A,&set,&flg));
623: symm = set? flg: PETSC_FALSE;
625: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
626: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
627: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
628: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
629: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
630: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
631: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
632: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
633: PetscCall(FNEvaluateFunctionMat_Sym_Default(fn,A,F));
634: } else {
635: /* scale argument */
636: if (fn->alpha!=(PetscScalar)1.0) {
637: PetscCall(FN_AllocateWorkMat(fn,A,&M));
638: PetscCall(MatScale(M,fn->alpha));
639: } else M = A;
640: PetscCall(FNEvaluateFunctionMat_Basic(fn,M,F));
641: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
642: /* scale result */
643: PetscCall(MatScale(F,fn->beta));
644: }
645: PetscCall(PetscFPTrapPop());
646: }
647: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { /* synchronize */
648: PetscCall(MatGetSize(A,&m,&n));
649: PetscCall(MatDenseGetArray(F,&pF));
650: PetscCall(PetscMPIIntCast(n*n,&n2));
651: PetscCallMPI(MPI_Bcast(pF,n2,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
652: PetscCall(MatDenseRestoreArray(F,&pF));
653: }
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: /*@
658: FNEvaluateFunctionMat - Computes the value of the function $f(A)$ for a given
659: matrix $A$, where the result is also a matrix.
661: Logically Collective
663: Input Parameters:
664: + fn - the math function context
665: - A - matrix on which the function must be evaluated
667: Output Parameter:
668: . B - (optional) matrix resulting from evaluating $f(A)$
670: Notes:
671: Matrix `A` must be a square sequential dense `Mat`, with all entries equal on
672: all processes (otherwise each process will compute different results).
673: If matrix `B` is provided, it must also be a square sequential dense `Mat`, and
674: both matrices must have the same dimensions. If `B` is `NULL` (or `B`=`A`) then
675: the function will perform an in-place computation, overwriting `A` with $f(A)$.
677: If `A` is known to be real symmetric or complex Hermitian then it is
678: recommended to set the appropriate flag with `MatSetOption()`, because
679: symmetry can sometimes be exploited by the algorithm.
681: Scaling factors are taken into account, so the actual function evaluation
682: will return $\beta f(\alpha A)$.
684: Level: advanced
686: .seealso: [](sec:fn), `FNEvaluateFunction()`, `FNEvaluateFunctionMatVec()`, `FNSetMethod()`
687: @*/
688: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
689: {
690: PetscBool inplace=PETSC_FALSE;
691: PetscInt m,n,n1;
692: MatType type;
694: PetscFunctionBegin;
699: if (B) {
702: } else inplace = PETSC_TRUE;
703: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA); //SlepcMatCheckSeq(A);
704: PetscCall(MatGetSize(A,&m,&n));
705: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
706: if (!inplace) {
707: PetscCall(MatGetType(A,&type));
708: PetscCheckTypeName(B,type);
709: n1 = n;
710: PetscCall(MatGetSize(B,&m,&n));
711: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat B is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
712: PetscCheck(n1==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrices A and B must have the same dimension");
713: }
715: /* evaluate matrix function */
716: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
717: PetscCall(FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE));
718: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }
722: /*
723: FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
724: and then copies the first column.
725: */
726: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
727: {
728: Mat F;
730: PetscFunctionBegin;
731: PetscCall(FN_AllocateWorkMat(fn,A,&F));
732: PetscCall(FNEvaluateFunctionMat_Basic(fn,A,F));
733: PetscCall(MatGetColumnVector(F,v,0));
734: PetscCall(FN_FreeWorkMat(fn,&F));
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: /*
739: FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
740: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
741: decomposition of A is A=Q*D*Q'. Only the first column is computed.
742: */
743: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
744: {
745: PetscInt m;
746: const PetscScalar *As;
747: PetscScalar *vs;
749: PetscFunctionBegin;
750: PetscCall(MatDenseGetArrayRead(A,&As));
751: PetscCall(VecGetArray(v,&vs));
752: PetscCall(MatGetSize(A,&m,NULL));
753: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE));
754: PetscCall(MatDenseRestoreArrayRead(A,&As));
755: PetscCall(VecRestoreArray(v,&vs));
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
760: {
761: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
762: PetscInt m,n;
763: Mat M;
764: PetscMPIInt size,rank,n_;
765: PetscScalar *pv;
767: PetscFunctionBegin;
768: /* check symmetry of A */
769: PetscCall(MatIsHermitianKnown(A,&set,&flg));
770: symm = set? flg: PETSC_FALSE;
772: /* evaluate matrix function */
773: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
774: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
775: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
776: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
777: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
778: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
779: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
780: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
781: PetscCall(FNEvaluateFunctionMatVec_Sym_Default(fn,A,v));
782: } else {
783: /* scale argument */
784: if (fn->alpha!=(PetscScalar)1.0) {
785: PetscCall(FN_AllocateWorkMat(fn,A,&M));
786: PetscCall(MatScale(M,fn->alpha));
787: } else M = A;
788: if (iscuda && fn->ops->evaluatefunctionmatveccuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatveccuda[fn->method],M,v);
789: else if (fn->ops->evaluatefunctionmatvec[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatvec[fn->method],M,v);
790: else PetscCall(FNEvaluateFunctionMatVec_Default(fn,M,v));
791: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
792: /* scale result */
793: PetscCall(VecScale(v,fn->beta));
794: }
795: PetscCall(PetscFPTrapPop());
796: }
798: /* synchronize */
799: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
800: PetscCall(MatGetSize(A,&m,&n));
801: PetscCall(VecGetArray(v,&pv));
802: PetscCall(PetscMPIIntCast(n,&n_));
803: PetscCallMPI(MPI_Bcast(pv,n_,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
804: PetscCall(VecRestoreArray(v,&pv));
805: }
806: PetscFunctionReturn(PETSC_SUCCESS);
807: }
809: /*@
810: FNEvaluateFunctionMatVec - Computes the first column of the matrix $f(A)$
811: for a given matrix $A$.
813: Logically Collective
815: Input Parameters:
816: + fn - the math function context
817: - A - matrix on which the function must be evaluated
819: Output Parameter:
820: . v - vector to hold the first column of $f(A)$
822: Notes:
823: This operation is similar to `FNEvaluateFunctionMat()` but returns only
824: the first column of $f(A)$, hence saving computations in most cases.
826: Level: advanced
828: .seealso: [](sec:fn), `FNEvaluateFunction()`, `FNEvaluateFunctionMat()`, `FNSetMethod()`
829: @*/
830: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
831: {
832: PetscInt m,n;
833: PetscBool iscuda;
835: PetscFunctionBegin;
842: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA); //SlepcMatCheckSeq(A);
843: PetscCall(MatGetSize(A,&m,&n));
844: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
845: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
846: PetscCheckTypeName(v,iscuda?VECSEQCUDA:VECSEQ);
847: PetscCall(VecGetSize(v,&m));
848: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrix A and vector v must have the same size");
849: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
850: PetscCall(FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE));
851: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: /*@
856: FNSetFromOptions - Sets `FN` options from the options database.
858: Collective
860: Input Parameter:
861: . fn - the math function context
863: Note:
864: To see all options, run your program with the `-help` option.
866: Level: beginner
868: .seealso: [](sec:fn), `FNSetOptionsPrefix()`
869: @*/
870: PetscErrorCode FNSetFromOptions(FN fn)
871: {
872: char type[256];
873: PetscScalar array[2];
874: PetscInt k,meth;
875: PetscBool flg;
876: FNParallelType pmode;
878: PetscFunctionBegin;
880: PetscCall(FNRegisterAll());
881: PetscObjectOptionsBegin((PetscObject)fn);
882: PetscCall(PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg));
883: if (flg) PetscCall(FNSetType(fn,type));
884: else if (!((PetscObject)fn)->type_name) PetscCall(FNSetType(fn,FNRATIONAL));
886: k = 2;
887: array[0] = 0.0; array[1] = 0.0;
888: PetscCall(PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg));
889: if (flg) {
890: if (k<2) array[1] = 1.0;
891: PetscCall(FNSetScale(fn,array[0],array[1]));
892: }
894: PetscCall(PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg));
895: if (flg) PetscCall(FNSetMethod(fn,meth));
897: PetscCall(PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg));
898: if (flg) PetscCall(FNSetParallel(fn,pmode));
900: PetscTryTypeMethod(fn,setfromoptions,PetscOptionsObject);
901: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fn,PetscOptionsObject));
902: PetscOptionsEnd();
903: PetscFunctionReturn(PETSC_SUCCESS);
904: }
906: /*@
907: FNView - Prints the `FN` data structure.
909: Collective
911: Input Parameters:
912: + fn - the math function context
913: - viewer - optional visualization context
915: Note:
916: The available visualization contexts include
917: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
918: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
919: first process opens the file; all other processes send their data to the
920: first one to print
922: The user can open an alternative visualization context with `PetscViewerASCIIOpen()`
923: to output to a specified file.
925: Use `FNViewFromOptions()` to allow the user to select many different `PetscViewerType`
926: and formats from the options database.
928: Level: beginner
930: .seealso: [](sec:fn), `FNCreate()`, `FNViewFromOptions()`
931: @*/
932: PetscErrorCode FNView(FN fn,PetscViewer viewer)
933: {
934: PetscBool isascii;
935: PetscMPIInt size;
937: PetscFunctionBegin;
939: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer));
941: PetscCheckSameComm(fn,1,viewer,2);
942: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
943: if (isascii) {
944: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer));
945: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
946: if (size>1) PetscCall(PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",FNParallelTypes[fn->pmode]));
947: PetscCall(PetscViewerASCIIPushTab(viewer));
948: PetscTryTypeMethod(fn,view,viewer);
949: PetscCall(PetscViewerASCIIPopTab(viewer));
950: }
951: PetscFunctionReturn(PETSC_SUCCESS);
952: }
954: /*@
955: FNViewFromOptions - View (print) an `FN` object based on values in the options database.
957: Collective
959: Input Parameters:
960: + fn - the math function context
961: . obj - optional object that provides the options prefix used to query the options database
962: - name - command line option
964: Level: intermediate
966: .seealso: [](sec:fn), `FNView()`, `FNCreate()`, `PetscObjectViewFromOptions()`
967: @*/
968: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])
969: {
970: PetscFunctionBegin;
972: PetscCall(PetscObjectViewFromOptions((PetscObject)fn,obj,name));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: /*@
977: FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
978: different communicator.
980: Collective
982: Input Parameters:
983: + fn - the math function context
984: - comm - MPI communicator
986: Output Parameter:
987: . newfn - location to put the new `FN` context
989: Note:
990: In order to use the same MPI communicator as in the original object,
991: use `PetscObjectComm`((`PetscObject`)`fn`).
993: Level: developer
995: .seealso: [](sec:fn), `FNCreate()`
996: @*/
997: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)
998: {
999: FNType type;
1000: PetscScalar alpha,beta;
1001: PetscInt meth;
1002: FNParallelType ptype;
1004: PetscFunctionBegin;
1007: PetscAssertPointer(newfn,3);
1008: PetscCall(FNCreate(comm,newfn));
1009: PetscCall(FNGetType(fn,&type));
1010: PetscCall(FNSetType(*newfn,type));
1011: PetscCall(FNGetScale(fn,&alpha,&beta));
1012: PetscCall(FNSetScale(*newfn,alpha,beta));
1013: PetscCall(FNGetMethod(fn,&meth));
1014: PetscCall(FNSetMethod(*newfn,meth));
1015: PetscCall(FNGetParallel(fn,&ptype));
1016: PetscCall(FNSetParallel(*newfn,ptype));
1017: PetscTryTypeMethod(fn,duplicate,comm,newfn);
1018: PetscFunctionReturn(PETSC_SUCCESS);
1019: }
1021: /*@
1022: FNDestroy - Destroys an `FN` context that was created with `FNCreate()`.
1024: Collective
1026: Input Parameter:
1027: . fn - the math function context
1029: Level: beginner
1031: .seealso: [](sec:fn), `FNCreate()`
1032: @*/
1033: PetscErrorCode FNDestroy(FN *fn)
1034: {
1035: PetscInt i;
1037: PetscFunctionBegin;
1038: if (!*fn) PetscFunctionReturn(PETSC_SUCCESS);
1040: if (--((PetscObject)*fn)->refct > 0) { *fn = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
1041: PetscTryTypeMethod(*fn,destroy);
1042: for (i=0;i<(*fn)->nw;i++) PetscCall(MatDestroy(&(*fn)->W[i]));
1043: PetscCall(PetscHeaderDestroy(fn));
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: /*@C
1048: FNRegister - Adds a mathematical function to the `FN` package.
1050: Not Collective
1052: Input Parameters:
1053: + name - name of a new user-defined `FN`
1054: - function - routine to create the context
1056: Notes:
1057: `FNRegister()` may be called multiple times to add several user-defined functions.
1059: Level: advanced
1061: .seealso: [](sec:fn), `FNRegisterAll()`
1062: @*/
1063: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1064: {
1065: PetscFunctionBegin;
1066: PetscCall(FNInitializePackage());
1067: PetscCall(PetscFunctionListAdd(&FNList,name,function));
1068: PetscFunctionReturn(PETSC_SUCCESS);
1069: }