Actual source code: fnbasic.c
slepc-3.21.0 2024-03-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: Basic FN routines
12: */
14: #include <slepc/private/fnimpl.h>
15: #include <slepcblaslapack.h>
17: PetscFunctionList FNList = NULL;
18: PetscBool FNRegisterAllCalled = PETSC_FALSE;
19: PetscClassId FN_CLASSID = 0;
20: PetscLogEvent FN_Evaluate = 0;
21: static PetscBool FNPackageInitialized = PETSC_FALSE;
23: const char *FNParallelTypes[] = {"REDUNDANT","SYNCHRONIZED","FNParallelType","FN_PARALLEL_",NULL};
25: /*@C
26: FNFinalizePackage - This function destroys everything in the Slepc interface
27: to the FN package. It is called from SlepcFinalize().
29: Level: developer
31: .seealso: SlepcFinalize()
32: @*/
33: PetscErrorCode FNFinalizePackage(void)
34: {
35: PetscFunctionBegin;
36: PetscCall(PetscFunctionListDestroy(&FNList));
37: FNPackageInitialized = PETSC_FALSE;
38: FNRegisterAllCalled = PETSC_FALSE;
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: /*@C
43: FNInitializePackage - This function initializes everything in the FN package.
44: It is called from PetscDLLibraryRegister() when using dynamic libraries, and
45: on the first call to FNCreate() when using static libraries.
47: Level: developer
49: .seealso: SlepcInitialize()
50: @*/
51: PetscErrorCode FNInitializePackage(void)
52: {
53: char logList[256];
54: PetscBool opt,pkg;
55: PetscClassId classids[1];
57: PetscFunctionBegin;
58: if (FNPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
59: FNPackageInitialized = PETSC_TRUE;
60: /* Register Classes */
61: PetscCall(PetscClassIdRegister("Math Function",&FN_CLASSID));
62: /* Register Constructors */
63: PetscCall(FNRegisterAll());
64: /* Register Events */
65: PetscCall(PetscLogEventRegister("FNEvaluate",FN_CLASSID,&FN_Evaluate));
66: /* Process Info */
67: classids[0] = FN_CLASSID;
68: PetscCall(PetscInfoProcessClass("fn",1,&classids[0]));
69: /* Process summary exclusions */
70: PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
71: if (opt) {
72: PetscCall(PetscStrInList("fn",logList,',',&pkg));
73: if (pkg) PetscCall(PetscLogEventDeactivateClass(FN_CLASSID));
74: }
75: /* Register package finalizer */
76: PetscCall(PetscRegisterFinalize(FNFinalizePackage));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: /*@
81: FNCreate - Creates an FN context.
83: Collective
85: Input Parameter:
86: . comm - MPI communicator
88: Output Parameter:
89: . newfn - location to put the FN context
91: Level: beginner
93: .seealso: FNDestroy(), FN
94: @*/
95: PetscErrorCode FNCreate(MPI_Comm comm,FN *newfn)
96: {
97: FN fn;
99: PetscFunctionBegin;
100: PetscAssertPointer(newfn,2);
101: *newfn = NULL;
102: PetscCall(FNInitializePackage());
103: PetscCall(SlepcHeaderCreate(fn,FN_CLASSID,"FN","Math Function","FN",comm,FNDestroy,FNView));
105: fn->alpha = 1.0;
106: fn->beta = 1.0;
107: fn->method = 0;
109: fn->nw = 0;
110: fn->cw = 0;
111: fn->data = NULL;
113: *newfn = fn;
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: /*@C
118: FNSetOptionsPrefix - Sets the prefix used for searching for all
119: FN options in the database.
121: Logically Collective
123: Input Parameters:
124: + fn - the math function context
125: - prefix - the prefix string to prepend to all FN option requests
127: Notes:
128: A hyphen (-) must NOT be given at the beginning of the prefix name.
129: The first character of all runtime options is AUTOMATICALLY the
130: hyphen.
132: Level: advanced
134: .seealso: FNAppendOptionsPrefix()
135: @*/
136: PetscErrorCode FNSetOptionsPrefix(FN fn,const char *prefix)
137: {
138: PetscFunctionBegin;
140: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fn,prefix));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /*@C
145: FNAppendOptionsPrefix - Appends to the prefix used for searching for all
146: FN options in the database.
148: Logically Collective
150: Input Parameters:
151: + fn - the math function context
152: - prefix - the prefix string to prepend to all FN option requests
154: Notes:
155: A hyphen (-) must NOT be given at the beginning of the prefix name.
156: The first character of all runtime options is AUTOMATICALLY the hyphen.
158: Level: advanced
160: .seealso: FNSetOptionsPrefix()
161: @*/
162: PetscErrorCode FNAppendOptionsPrefix(FN fn,const char *prefix)
163: {
164: PetscFunctionBegin;
166: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)fn,prefix));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: /*@C
171: FNGetOptionsPrefix - Gets the prefix used for searching for all
172: FN options in the database.
174: Not Collective
176: Input Parameters:
177: . fn - the math function context
179: Output Parameters:
180: . prefix - pointer to the prefix string used is returned
182: Note:
183: On the Fortran side, the user should pass in a string 'prefix' of
184: sufficient length to hold the prefix.
186: Level: advanced
188: .seealso: FNSetOptionsPrefix(), FNAppendOptionsPrefix()
189: @*/
190: PetscErrorCode FNGetOptionsPrefix(FN fn,const char *prefix[])
191: {
192: PetscFunctionBegin;
194: PetscAssertPointer(prefix,2);
195: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)fn,prefix));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@C
200: FNSetType - Selects the type for the FN object.
202: Logically Collective
204: Input Parameters:
205: + fn - the math function context
206: - type - a known type
208: Notes:
209: The default is FNRATIONAL, which includes polynomials as a particular
210: case as well as simple functions such as f(x)=x and f(x)=constant.
212: Level: intermediate
214: .seealso: FNGetType()
215: @*/
216: PetscErrorCode FNSetType(FN fn,FNType type)
217: {
218: PetscErrorCode (*r)(FN);
219: PetscBool match;
221: PetscFunctionBegin;
223: PetscAssertPointer(type,2);
225: PetscCall(PetscObjectTypeCompare((PetscObject)fn,type,&match));
226: if (match) PetscFunctionReturn(PETSC_SUCCESS);
228: PetscCall(PetscFunctionListFind(FNList,type,&r));
229: PetscCheck(r,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested FN type %s",type);
231: PetscTryTypeMethod(fn,destroy);
232: PetscCall(PetscMemzero(fn->ops,sizeof(struct _FNOps)));
234: PetscCall(PetscObjectChangeTypeName((PetscObject)fn,type));
235: PetscCall((*r)(fn));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@C
240: FNGetType - Gets the FN type name (as a string) from the FN context.
242: Not Collective
244: Input Parameter:
245: . fn - the math function context
247: Output Parameter:
248: . type - name of the math function
250: Level: intermediate
252: .seealso: FNSetType()
253: @*/
254: PetscErrorCode FNGetType(FN fn,FNType *type)
255: {
256: PetscFunctionBegin;
258: PetscAssertPointer(type,2);
259: *type = ((PetscObject)fn)->type_name;
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: /*@
264: FNSetScale - Sets the scaling parameters that define the matematical function.
266: Logically Collective
268: Input Parameters:
269: + fn - the math function context
270: . alpha - inner scaling (argument)
271: - beta - outer scaling (result)
273: Notes:
274: Given a function f(x) specified by the FN type, the scaling parameters can
275: be used to realize the function beta*f(alpha*x). So when these values are given,
276: the procedure for function evaluation will first multiply the argument by alpha,
277: then evaluate the function itself, and finally scale the result by beta.
278: Likewise, these values are also considered when evaluating the derivative.
280: If you want to provide only one of the two scaling factors, set the other
281: one to 1.0.
283: Level: intermediate
285: .seealso: FNGetScale(), FNEvaluateFunction()
286: @*/
287: PetscErrorCode FNSetScale(FN fn,PetscScalar alpha,PetscScalar beta)
288: {
289: PetscFunctionBegin;
293: PetscCheck(PetscAbsScalar(alpha)!=0.0 && PetscAbsScalar(beta)!=0.0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_WRONG,"Scaling factors must be nonzero");
294: fn->alpha = alpha;
295: fn->beta = beta;
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: /*@
300: FNGetScale - Gets the scaling parameters that define the matematical function.
302: Not Collective
304: Input Parameter:
305: . fn - the math function context
307: Output Parameters:
308: + alpha - inner scaling (argument)
309: - beta - outer scaling (result)
311: Level: intermediate
313: .seealso: FNSetScale()
314: @*/
315: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)
316: {
317: PetscFunctionBegin;
319: if (alpha) *alpha = fn->alpha;
320: if (beta) *beta = fn->beta;
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: FNSetMethod - Selects the method to be used to evaluate functions of matrices.
327: Logically Collective
329: Input Parameters:
330: + fn - the math function context
331: - meth - an index identifying the method
333: Options Database Key:
334: . -fn_method <meth> - Sets the method
336: Notes:
337: In some FN types there are more than one algorithms available for computing
338: matrix functions. In that case, this function allows choosing the wanted method.
340: If meth is currently set to 0 (the default) and the input argument A of
341: FNEvaluateFunctionMat() is a symmetric/Hermitian matrix, then the computation
342: is done via the eigendecomposition of A, rather than with the general algorithm.
344: Level: intermediate
346: .seealso: FNGetMethod(), FNEvaluateFunctionMat()
347: @*/
348: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)
349: {
350: PetscFunctionBegin;
353: PetscCheck(meth>=0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
354: PetscCheck(meth<=FN_MAX_SOLVE,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
355: fn->method = meth;
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: /*@
360: FNGetMethod - Gets the method currently used in the FN.
362: Not Collective
364: Input Parameter:
365: . fn - the math function context
367: Output Parameter:
368: . meth - identifier of the method
370: Level: intermediate
372: .seealso: FNSetMethod()
373: @*/
374: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)
375: {
376: PetscFunctionBegin;
378: PetscAssertPointer(meth,2);
379: *meth = fn->method;
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: /*@
384: FNSetParallel - Selects the mode of operation in parallel runs.
386: Logically Collective
388: Input Parameters:
389: + fn - the math function context
390: - pmode - the parallel mode
392: Options Database Key:
393: . -fn_parallel <mode> - Sets the parallel mode, either 'redundant' or 'synchronized'
395: Notes:
396: This is relevant only when the function is evaluated on a matrix, with
397: either FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
399: In the 'redundant' parallel mode, all processes will make the computation
400: redundantly, starting from the same data, and producing the same result.
401: This result may be slightly different in the different processes if using a
402: multithreaded BLAS library, which may cause issues in ill-conditioned problems.
404: In the 'synchronized' parallel mode, only the first MPI process performs the
405: computation and then the computed matrix is broadcast to the other
406: processes in the communicator. This communication is done automatically at
407: the end of FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
409: Level: advanced
411: .seealso: FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec(), FNGetParallel()
412: @*/
413: PetscErrorCode FNSetParallel(FN fn,FNParallelType pmode)
414: {
415: PetscFunctionBegin;
418: fn->pmode = pmode;
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*@
423: FNGetParallel - Gets the mode of operation in parallel runs.
425: Not Collective
427: Input Parameter:
428: . fn - the math function context
430: Output Parameter:
431: . pmode - the parallel mode
433: Level: advanced
435: .seealso: FNSetParallel()
436: @*/
437: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)
438: {
439: PetscFunctionBegin;
441: PetscAssertPointer(pmode,2);
442: *pmode = fn->pmode;
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /*@
447: FNEvaluateFunction - Computes the value of the function f(x) for a given x.
449: Not Collective
451: Input Parameters:
452: + fn - the math function context
453: - x - the value where the function must be evaluated
455: Output Parameter:
456: . y - the result of f(x)
458: Note:
459: Scaling factors are taken into account, so the actual function evaluation
460: will return beta*f(alpha*x).
462: Level: intermediate
464: .seealso: FNEvaluateDerivative(), FNEvaluateFunctionMat(), FNSetScale()
465: @*/
466: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)
467: {
468: PetscScalar xf,yf;
470: PetscFunctionBegin;
473: PetscAssertPointer(y,3);
474: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
475: xf = fn->alpha*x;
476: PetscUseTypeMethod(fn,evaluatefunction,xf,&yf);
477: *y = fn->beta*yf;
478: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: /*@
483: FNEvaluateDerivative - Computes the value of the derivative f'(x) for a given x.
485: Not Collective
487: Input Parameters:
488: + fn - the math function context
489: - x - the value where the derivative must be evaluated
491: Output Parameter:
492: . y - the result of f'(x)
494: Note:
495: Scaling factors are taken into account, so the actual derivative evaluation will
496: return alpha*beta*f'(alpha*x).
498: Level: intermediate
500: .seealso: FNEvaluateFunction()
501: @*/
502: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)
503: {
504: PetscScalar xf,yf;
506: PetscFunctionBegin;
509: PetscAssertPointer(y,3);
510: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
511: xf = fn->alpha*x;
512: PetscUseTypeMethod(fn,evaluatederivative,xf,&yf);
513: *y = fn->alpha*fn->beta*yf;
514: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,const PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)
519: {
520: PetscInt i,j;
521: PetscBLASInt n,k,ld,lwork,info;
522: PetscScalar *Q,*W,*work,adummy,a,x,y,one=1.0,zero=0.0;
523: PetscReal *eig,dummy;
524: #if defined(PETSC_USE_COMPLEX)
525: PetscReal *rwork,rdummy;
526: #endif
528: PetscFunctionBegin;
529: PetscCall(PetscBLASIntCast(m,&n));
530: ld = n;
531: k = firstonly? 1: n;
533: /* workspace query and memory allocation */
534: lwork = -1;
535: #if defined(PETSC_USE_COMPLEX)
536: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&rdummy,&info));
537: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
538: PetscCall(PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork));
539: #else
540: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&info));
541: PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
542: PetscCall(PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work));
543: #endif
545: /* compute eigendecomposition */
546: for (j=0;j<n;j++) for (i=j;i<n;i++) Q[i+j*ld] = As[i+j*ld];
547: #if defined(PETSC_USE_COMPLEX)
548: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
549: #else
550: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
551: #endif
552: SlepcCheckLapackInfo("syev",info);
554: /* W = f(Lambda)*Q' */
555: for (i=0;i<n;i++) {
556: x = fn->alpha*eig[i];
557: PetscUseTypeMethod(fn,evaluatefunction,x,&y); /* y = f(x) */
558: for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
559: }
560: /* Bs = Q*W */
561: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
562: #if defined(PETSC_USE_COMPLEX)
563: PetscCall(PetscFree5(eig,Q,W,work,rwork));
564: #else
565: PetscCall(PetscFree4(eig,Q,W,work));
566: #endif
567: PetscCall(PetscLogFlops(9.0*n*n*n+2.0*n*n*n));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: /*
572: FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
573: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
574: decomposition of A is A=Q*D*Q'
575: */
576: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)
577: {
578: PetscInt m;
579: const PetscScalar *As;
580: PetscScalar *Bs;
582: PetscFunctionBegin;
583: PetscCall(MatDenseGetArrayRead(A,&As));
584: PetscCall(MatDenseGetArray(B,&Bs));
585: PetscCall(MatGetSize(A,&m,NULL));
586: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE));
587: PetscCall(MatDenseRestoreArrayRead(A,&As));
588: PetscCall(MatDenseRestoreArray(B,&Bs));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: static PetscErrorCode FNEvaluateFunctionMat_Basic(FN fn,Mat A,Mat F)
593: {
594: PetscBool iscuda;
596: PetscFunctionBegin;
597: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
598: 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));
599: if (iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatcuda[fn->method],A,F);
600: else if (fn->ops->evaluatefunctionmat[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmat[fn->method],A,F);
601: else {
602: PetscCheck(fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
603: PetscCheck(!fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
604: }
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
609: {
610: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
611: PetscInt m,n;
612: PetscMPIInt size,rank;
613: PetscScalar *pF;
614: Mat M,F;
616: PetscFunctionBegin;
617: /* destination matrix */
618: F = B?B:A;
620: /* check symmetry of A */
621: PetscCall(MatIsHermitianKnown(A,&set,&flg));
622: symm = set? flg: PETSC_FALSE;
624: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
625: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
626: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
627: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
628: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
629: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
630: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
631: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
632: PetscCall(FNEvaluateFunctionMat_Sym_Default(fn,A,F));
633: } else {
634: /* scale argument */
635: if (fn->alpha!=(PetscScalar)1.0) {
636: PetscCall(FN_AllocateWorkMat(fn,A,&M));
637: PetscCall(MatScale(M,fn->alpha));
638: } else M = A;
639: PetscCall(FNEvaluateFunctionMat_Basic(fn,M,F));
640: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
641: /* scale result */
642: PetscCall(MatScale(F,fn->beta));
643: }
644: PetscCall(PetscFPTrapPop());
645: }
646: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { /* synchronize */
647: PetscCall(MatGetSize(A,&m,&n));
648: PetscCall(MatDenseGetArray(F,&pF));
649: PetscCallMPI(MPI_Bcast(pF,n*n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
650: PetscCall(MatDenseRestoreArray(F,&pF));
651: }
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: /*@
656: FNEvaluateFunctionMat - Computes the value of the function f(A) for a given
657: matrix A, where the result is also a matrix.
659: Logically Collective
661: Input Parameters:
662: + fn - the math function context
663: - A - matrix on which the function must be evaluated
665: Output Parameter:
666: . B - (optional) matrix resulting from evaluating f(A)
668: Notes:
669: Matrix A must be a square sequential dense Mat, with all entries equal on
670: all processes (otherwise each process will compute different results).
671: If matrix B is provided, it must also be a square sequential dense Mat, and
672: both matrices must have the same dimensions. If B is NULL (or B=A) then the
673: function will perform an in-place computation, overwriting A with f(A).
675: If A is known to be real symmetric or complex Hermitian then it is
676: recommended to set the appropriate flag with MatSetOption(), because
677: symmetry can sometimes be exploited by the algorithm.
679: Scaling factors are taken into account, so the actual function evaluation
680: will return beta*f(alpha*A).
682: Level: advanced
684: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
685: @*/
686: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
687: {
688: PetscBool inplace=PETSC_FALSE;
689: PetscInt m,n,n1;
690: MatType type;
692: PetscFunctionBegin;
697: if (B) {
700: } else inplace = PETSC_TRUE;
701: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA);
702: PetscCall(MatGetSize(A,&m,&n));
703: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
704: if (!inplace) {
705: PetscCall(MatGetType(A,&type));
706: PetscCheckTypeName(B,type);
707: n1 = n;
708: PetscCall(MatGetSize(B,&m,&n));
709: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat B is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
710: PetscCheck(n1==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrices A and B must have the same dimension");
711: }
713: /* evaluate matrix function */
714: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
715: PetscCall(FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE));
716: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: /*
721: FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
722: and then copies the first column.
723: */
724: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
725: {
726: Mat F;
728: PetscFunctionBegin;
729: PetscCall(FN_AllocateWorkMat(fn,A,&F));
730: PetscCall(FNEvaluateFunctionMat_Basic(fn,A,F));
731: PetscCall(MatGetColumnVector(F,v,0));
732: PetscCall(FN_FreeWorkMat(fn,&F));
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: /*
737: FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
738: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
739: decomposition of A is A=Q*D*Q'. Only the first column is computed.
740: */
741: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
742: {
743: PetscInt m;
744: const PetscScalar *As;
745: PetscScalar *vs;
747: PetscFunctionBegin;
748: PetscCall(MatDenseGetArrayRead(A,&As));
749: PetscCall(VecGetArray(v,&vs));
750: PetscCall(MatGetSize(A,&m,NULL));
751: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE));
752: PetscCall(MatDenseRestoreArrayRead(A,&As));
753: PetscCall(VecRestoreArray(v,&vs));
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
758: {
759: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
760: PetscInt m,n;
761: Mat M;
762: PetscMPIInt size,rank;
763: PetscScalar *pv;
765: PetscFunctionBegin;
766: /* check symmetry of A */
767: PetscCall(MatIsHermitianKnown(A,&set,&flg));
768: symm = set? flg: PETSC_FALSE;
770: /* evaluate matrix function */
771: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
772: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
773: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
774: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
775: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
776: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
777: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
778: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
779: PetscCall(FNEvaluateFunctionMatVec_Sym_Default(fn,A,v));
780: } else {
781: /* scale argument */
782: if (fn->alpha!=(PetscScalar)1.0) {
783: PetscCall(FN_AllocateWorkMat(fn,A,&M));
784: PetscCall(MatScale(M,fn->alpha));
785: } else M = A;
786: if (iscuda && fn->ops->evaluatefunctionmatveccuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatveccuda[fn->method],M,v);
787: else if (fn->ops->evaluatefunctionmatvec[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatvec[fn->method],M,v);
788: else PetscCall(FNEvaluateFunctionMatVec_Default(fn,M,v));
789: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
790: /* scale result */
791: PetscCall(VecScale(v,fn->beta));
792: }
793: PetscCall(PetscFPTrapPop());
794: }
796: /* synchronize */
797: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
798: PetscCall(MatGetSize(A,&m,&n));
799: PetscCall(VecGetArray(v,&pv));
800: PetscCallMPI(MPI_Bcast(pv,n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
801: PetscCall(VecRestoreArray(v,&pv));
802: }
803: PetscFunctionReturn(PETSC_SUCCESS);
804: }
806: /*@
807: FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
808: for a given matrix A.
810: Logically Collective
812: Input Parameters:
813: + fn - the math function context
814: - A - matrix on which the function must be evaluated
816: Output Parameter:
817: . v - vector to hold the first column of f(A)
819: Notes:
820: This operation is similar to FNEvaluateFunctionMat() but returns only
821: the first column of f(A), hence saving computations in most cases.
823: Level: advanced
825: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
826: @*/
827: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
828: {
829: PetscInt m,n;
830: PetscBool iscuda;
832: PetscFunctionBegin;
839: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA);
840: PetscCall(MatGetSize(A,&m,&n));
841: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
842: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
843: PetscCheckTypeName(v,iscuda?VECSEQCUDA:VECSEQ);
844: PetscCall(VecGetSize(v,&m));
845: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrix A and vector v must have the same size");
846: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
847: PetscCall(FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE));
848: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: /*@
853: FNSetFromOptions - Sets FN options from the options database.
855: Collective
857: Input Parameters:
858: . fn - the math function context
860: Notes:
861: To see all options, run your program with the -help option.
863: Level: beginner
865: .seealso: FNSetOptionsPrefix()
866: @*/
867: PetscErrorCode FNSetFromOptions(FN fn)
868: {
869: char type[256];
870: PetscScalar array[2];
871: PetscInt k,meth;
872: PetscBool flg;
873: FNParallelType pmode;
875: PetscFunctionBegin;
877: PetscCall(FNRegisterAll());
878: PetscObjectOptionsBegin((PetscObject)fn);
879: PetscCall(PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg));
880: if (flg) PetscCall(FNSetType(fn,type));
881: else if (!((PetscObject)fn)->type_name) PetscCall(FNSetType(fn,FNRATIONAL));
883: k = 2;
884: array[0] = 0.0; array[1] = 0.0;
885: PetscCall(PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg));
886: if (flg) {
887: if (k<2) array[1] = 1.0;
888: PetscCall(FNSetScale(fn,array[0],array[1]));
889: }
891: PetscCall(PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg));
892: if (flg) PetscCall(FNSetMethod(fn,meth));
894: PetscCall(PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg));
895: if (flg) PetscCall(FNSetParallel(fn,pmode));
897: PetscTryTypeMethod(fn,setfromoptions,PetscOptionsObject);
898: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fn,PetscOptionsObject));
899: PetscOptionsEnd();
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: /*@C
904: FNView - Prints the FN data structure.
906: Collective
908: Input Parameters:
909: + fn - the math function context
910: - viewer - optional visualization context
912: Note:
913: The available visualization contexts include
914: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
915: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
916: output where only the first processor opens
917: the file. All other processors send their
918: data to the first processor to print.
920: The user can open an alternative visualization context with
921: PetscViewerASCIIOpen() - output to a specified file.
923: Level: beginner
925: .seealso: FNCreate()
926: @*/
927: PetscErrorCode FNView(FN fn,PetscViewer viewer)
928: {
929: PetscBool isascii;
930: PetscMPIInt size;
932: PetscFunctionBegin;
934: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer));
936: PetscCheckSameComm(fn,1,viewer,2);
937: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
938: if (isascii) {
939: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer));
940: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
941: if (size>1) PetscCall(PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",FNParallelTypes[fn->pmode]));
942: PetscCall(PetscViewerASCIIPushTab(viewer));
943: PetscTryTypeMethod(fn,view,viewer);
944: PetscCall(PetscViewerASCIIPopTab(viewer));
945: }
946: PetscFunctionReturn(PETSC_SUCCESS);
947: }
949: /*@C
950: FNViewFromOptions - View from options
952: Collective
954: Input Parameters:
955: + fn - the math function context
956: . obj - optional object
957: - name - command line option
959: Level: intermediate
961: .seealso: FNView(), FNCreate()
962: @*/
963: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])
964: {
965: PetscFunctionBegin;
967: PetscCall(PetscObjectViewFromOptions((PetscObject)fn,obj,name));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@
972: FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
973: different communicator.
975: Collective
977: Input Parameters:
978: + fn - the math function context
979: - comm - MPI communicator
981: Output Parameter:
982: . newfn - location to put the new FN context
984: Note:
985: In order to use the same MPI communicator as in the original object,
986: use PetscObjectComm((PetscObject)fn).
988: Level: developer
990: .seealso: FNCreate()
991: @*/
992: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)
993: {
994: FNType type;
995: PetscScalar alpha,beta;
996: PetscInt meth;
997: FNParallelType ptype;
999: PetscFunctionBegin;
1002: PetscAssertPointer(newfn,3);
1003: PetscCall(FNCreate(comm,newfn));
1004: PetscCall(FNGetType(fn,&type));
1005: PetscCall(FNSetType(*newfn,type));
1006: PetscCall(FNGetScale(fn,&alpha,&beta));
1007: PetscCall(FNSetScale(*newfn,alpha,beta));
1008: PetscCall(FNGetMethod(fn,&meth));
1009: PetscCall(FNSetMethod(*newfn,meth));
1010: PetscCall(FNGetParallel(fn,&ptype));
1011: PetscCall(FNSetParallel(*newfn,ptype));
1012: PetscTryTypeMethod(fn,duplicate,comm,newfn);
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1016: /*@C
1017: FNDestroy - Destroys FN context that was created with FNCreate().
1019: Collective
1021: Input Parameter:
1022: . fn - the math function context
1024: Level: beginner
1026: .seealso: FNCreate()
1027: @*/
1028: PetscErrorCode FNDestroy(FN *fn)
1029: {
1030: PetscInt i;
1032: PetscFunctionBegin;
1033: if (!*fn) PetscFunctionReturn(PETSC_SUCCESS);
1035: if (--((PetscObject)*fn)->refct > 0) { *fn = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
1036: PetscTryTypeMethod(*fn,destroy);
1037: for (i=0;i<(*fn)->nw;i++) PetscCall(MatDestroy(&(*fn)->W[i]));
1038: PetscCall(PetscHeaderDestroy(fn));
1039: PetscFunctionReturn(PETSC_SUCCESS);
1040: }
1042: /*@C
1043: FNRegister - Adds a mathematical function to the FN package.
1045: Not Collective
1047: Input Parameters:
1048: + name - name of a new user-defined FN
1049: - function - routine to create context
1051: Notes:
1052: FNRegister() may be called multiple times to add several user-defined functions.
1054: Level: advanced
1056: .seealso: FNRegisterAll()
1057: @*/
1058: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1059: {
1060: PetscFunctionBegin;
1061: PetscCall(FNInitializePackage());
1062: PetscCall(PetscFunctionListAdd(&FNList,name,function));
1063: PetscFunctionReturn(PETSC_SUCCESS);
1064: }