Actual source code: fnbasic.c
slepc-main 2024-12-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: Basic FN routines
12: */
14: #include <slepc/private/fnimpl.h>
15: #include <slepcblaslapack.h>
17: PetscFunctionList FNList = NULL;
18: PetscBool FNRegisterAllCalled = PETSC_FALSE;
19: PetscClassId FN_CLASSID = 0;
20: PetscLogEvent FN_Evaluate = 0;
21: static PetscBool FNPackageInitialized = PETSC_FALSE;
23: const char *FNParallelTypes[] = {"REDUNDANT","SYNCHRONIZED","FNParallelType","FN_PARALLEL_",NULL};
25: /*@C
26: FNFinalizePackage - This function destroys everything in the Slepc interface
27: to the FN package. It is called from SlepcFinalize().
29: Level: developer
31: .seealso: SlepcFinalize()
32: @*/
33: PetscErrorCode FNFinalizePackage(void)
34: {
35: PetscFunctionBegin;
36: PetscCall(PetscFunctionListDestroy(&FNList));
37: FNPackageInitialized = PETSC_FALSE;
38: FNRegisterAllCalled = PETSC_FALSE;
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: /*@C
43: FNInitializePackage - This function initializes everything in the FN package.
44: It is called from PetscDLLibraryRegister() when using dynamic libraries, and
45: on the first call to FNCreate() when using static libraries.
47: Level: developer
49: .seealso: SlepcInitialize()
50: @*/
51: PetscErrorCode FNInitializePackage(void)
52: {
53: char logList[256];
54: PetscBool opt,pkg;
55: PetscClassId classids[1];
57: PetscFunctionBegin;
58: if (FNPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
59: FNPackageInitialized = PETSC_TRUE;
60: /* Register Classes */
61: PetscCall(PetscClassIdRegister("Math Function",&FN_CLASSID));
62: /* Register Constructors */
63: PetscCall(FNRegisterAll());
64: /* Register Events */
65: PetscCall(PetscLogEventRegister("FNEvaluate",FN_CLASSID,&FN_Evaluate));
66: /* Process Info */
67: classids[0] = FN_CLASSID;
68: PetscCall(PetscInfoProcessClass("fn",1,&classids[0]));
69: /* Process summary exclusions */
70: PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
71: if (opt) {
72: PetscCall(PetscStrInList("fn",logList,',',&pkg));
73: if (pkg) PetscCall(PetscLogEventDeactivateClass(FN_CLASSID));
74: }
75: /* Register package finalizer */
76: PetscCall(PetscRegisterFinalize(FNFinalizePackage));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: /*@
81: FNCreate - Creates an FN context.
83: Collective
85: Input Parameter:
86: . comm - MPI communicator
88: Output Parameter:
89: . newfn - location to put the FN context
91: Level: beginner
93: .seealso: FNDestroy(), FN
94: @*/
95: PetscErrorCode FNCreate(MPI_Comm comm,FN *newfn)
96: {
97: FN fn;
99: PetscFunctionBegin;
100: PetscAssertPointer(newfn,2);
101: PetscCall(FNInitializePackage());
102: PetscCall(SlepcHeaderCreate(fn,FN_CLASSID,"FN","Math Function","FN",comm,FNDestroy,FNView));
104: fn->alpha = 1.0;
105: fn->beta = 1.0;
106: fn->method = 0;
108: fn->nw = 0;
109: fn->cw = 0;
110: fn->data = NULL;
112: *newfn = fn;
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /*@
117: FNSetOptionsPrefix - Sets the prefix used for searching for all
118: FN options in the database.
120: Logically Collective
122: Input Parameters:
123: + fn - the math function context
124: - prefix - the prefix string to prepend to all FN option requests
126: Notes:
127: A hyphen (-) must NOT be given at the beginning of the prefix name.
128: The first character of all runtime options is AUTOMATICALLY the
129: hyphen.
131: Level: advanced
133: .seealso: FNAppendOptionsPrefix()
134: @*/
135: PetscErrorCode FNSetOptionsPrefix(FN fn,const char *prefix)
136: {
137: PetscFunctionBegin;
139: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fn,prefix));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: /*@
144: FNAppendOptionsPrefix - Appends to the prefix used for searching for all
145: FN options in the database.
147: Logically Collective
149: Input Parameters:
150: + fn - the math function context
151: - prefix - the prefix string to prepend to all FN option requests
153: Notes:
154: A hyphen (-) must NOT be given at the beginning of the prefix name.
155: The first character of all runtime options is AUTOMATICALLY the hyphen.
157: Level: advanced
159: .seealso: FNSetOptionsPrefix()
160: @*/
161: PetscErrorCode FNAppendOptionsPrefix(FN fn,const char *prefix)
162: {
163: PetscFunctionBegin;
165: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)fn,prefix));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: /*@
170: FNGetOptionsPrefix - Gets the prefix used for searching for all
171: FN options in the database.
173: Not Collective
175: Input Parameters:
176: . fn - the math function context
178: Output Parameters:
179: . prefix - pointer to the prefix string used is returned
181: Note:
182: On the Fortran side, the user should pass in a string 'prefix' of
183: sufficient length to hold the prefix.
185: Level: advanced
187: .seealso: FNSetOptionsPrefix(), FNAppendOptionsPrefix()
188: @*/
189: PetscErrorCode FNGetOptionsPrefix(FN fn,const char *prefix[])
190: {
191: PetscFunctionBegin;
193: PetscAssertPointer(prefix,2);
194: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)fn,prefix));
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: /*@
199: FNSetType - Selects the type for the FN object.
201: Logically Collective
203: Input Parameters:
204: + fn - the math function context
205: - type - a known type
207: Notes:
208: The default is FNRATIONAL, which includes polynomials as a particular
209: case as well as simple functions such as f(x)=x and f(x)=constant.
211: Level: intermediate
213: .seealso: FNGetType()
214: @*/
215: PetscErrorCode FNSetType(FN fn,FNType type)
216: {
217: PetscErrorCode (*r)(FN);
218: PetscBool match;
220: PetscFunctionBegin;
222: PetscAssertPointer(type,2);
224: PetscCall(PetscObjectTypeCompare((PetscObject)fn,type,&match));
225: if (match) PetscFunctionReturn(PETSC_SUCCESS);
227: PetscCall(PetscFunctionListFind(FNList,type,&r));
228: PetscCheck(r,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested FN type %s",type);
230: PetscTryTypeMethod(fn,destroy);
231: PetscCall(PetscMemzero(fn->ops,sizeof(struct _FNOps)));
233: PetscCall(PetscObjectChangeTypeName((PetscObject)fn,type));
234: PetscCall((*r)(fn));
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@
239: FNGetType - Gets the FN type name (as a string) from the FN context.
241: Not Collective
243: Input Parameter:
244: . fn - the math function context
246: Output Parameter:
247: . type - name of the math function
249: Level: intermediate
251: .seealso: FNSetType()
252: @*/
253: PetscErrorCode FNGetType(FN fn,FNType *type)
254: {
255: PetscFunctionBegin;
257: PetscAssertPointer(type,2);
258: *type = ((PetscObject)fn)->type_name;
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: /*@
263: FNSetScale - Sets the scaling parameters that define the matematical function.
265: Logically Collective
267: Input Parameters:
268: + fn - the math function context
269: . alpha - inner scaling (argument)
270: - beta - outer scaling (result)
272: Notes:
273: Given a function f(x) specified by the FN type, the scaling parameters can
274: be used to realize the function beta*f(alpha*x). So when these values are given,
275: the procedure for function evaluation will first multiply the argument by alpha,
276: then evaluate the function itself, and finally scale the result by beta.
277: Likewise, these values are also considered when evaluating the derivative.
279: If you want to provide only one of the two scaling factors, set the other
280: one to 1.0.
282: Level: intermediate
284: .seealso: FNGetScale(), FNEvaluateFunction()
285: @*/
286: PetscErrorCode FNSetScale(FN fn,PetscScalar alpha,PetscScalar beta)
287: {
288: PetscFunctionBegin;
292: PetscCheck(PetscAbsScalar(alpha)!=0.0 && PetscAbsScalar(beta)!=0.0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_WRONG,"Scaling factors must be nonzero");
293: fn->alpha = alpha;
294: fn->beta = beta;
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: /*@
299: FNGetScale - Gets the scaling parameters that define the matematical function.
301: Not Collective
303: Input Parameter:
304: . fn - the math function context
306: Output Parameters:
307: + alpha - inner scaling (argument)
308: - beta - outer scaling (result)
310: Level: intermediate
312: .seealso: FNSetScale()
313: @*/
314: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)
315: {
316: PetscFunctionBegin;
318: if (alpha) *alpha = fn->alpha;
319: if (beta) *beta = fn->beta;
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /*@
324: FNSetMethod - Selects the method to be used to evaluate functions of matrices.
326: Logically Collective
328: Input Parameters:
329: + fn - the math function context
330: - meth - an index identifying the method
332: Options Database Key:
333: . -fn_method <meth> - Sets the method
335: Notes:
336: In some FN types there are more than one algorithms available for computing
337: matrix functions. In that case, this function allows choosing the wanted method.
339: If meth is currently set to 0 (the default) and the input argument A of
340: FNEvaluateFunctionMat() is a symmetric/Hermitian matrix, then the computation
341: is done via the eigendecomposition of A, rather than with the general algorithm.
343: Level: intermediate
345: .seealso: FNGetMethod(), FNEvaluateFunctionMat()
346: @*/
347: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)
348: {
349: PetscFunctionBegin;
352: PetscCheck(meth>=0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
353: PetscCheck(meth<=FN_MAX_SOLVE,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
354: fn->method = meth;
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: /*@
359: FNGetMethod - Gets the method currently used in the FN.
361: Not Collective
363: Input Parameter:
364: . fn - the math function context
366: Output Parameter:
367: . meth - identifier of the method
369: Level: intermediate
371: .seealso: FNSetMethod()
372: @*/
373: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)
374: {
375: PetscFunctionBegin;
377: PetscAssertPointer(meth,2);
378: *meth = fn->method;
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: /*@
383: FNSetParallel - Selects the mode of operation in parallel runs.
385: Logically Collective
387: Input Parameters:
388: + fn - the math function context
389: - pmode - the parallel mode
391: Options Database Key:
392: . -fn_parallel <mode> - Sets the parallel mode, either 'redundant' or 'synchronized'
394: Notes:
395: This is relevant only when the function is evaluated on a matrix, with
396: either FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
398: In the 'redundant' parallel mode, all processes will make the computation
399: redundantly, starting from the same data, and producing the same result.
400: This result may be slightly different in the different processes if using a
401: multithreaded BLAS library, which may cause issues in ill-conditioned problems.
403: In the 'synchronized' parallel mode, only the first MPI process performs the
404: computation and then the computed matrix is broadcast to the other
405: processes in the communicator. This communication is done automatically at
406: the end of FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
408: Level: advanced
410: .seealso: FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec(), FNGetParallel()
411: @*/
412: PetscErrorCode FNSetParallel(FN fn,FNParallelType pmode)
413: {
414: PetscFunctionBegin;
417: fn->pmode = pmode;
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: /*@
422: FNGetParallel - Gets the mode of operation in parallel runs.
424: Not Collective
426: Input Parameter:
427: . fn - the math function context
429: Output Parameter:
430: . pmode - the parallel mode
432: Level: advanced
434: .seealso: FNSetParallel()
435: @*/
436: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)
437: {
438: PetscFunctionBegin;
440: PetscAssertPointer(pmode,2);
441: *pmode = fn->pmode;
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: /*@
446: FNEvaluateFunction - Computes the value of the function f(x) for a given x.
448: Not Collective
450: Input Parameters:
451: + fn - the math function context
452: - x - the value where the function must be evaluated
454: Output Parameter:
455: . y - the result of f(x)
457: Note:
458: Scaling factors are taken into account, so the actual function evaluation
459: will return beta*f(alpha*x).
461: Level: intermediate
463: .seealso: FNEvaluateDerivative(), FNEvaluateFunctionMat(), FNSetScale()
464: @*/
465: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)
466: {
467: PetscScalar xf,yf;
469: PetscFunctionBegin;
472: PetscAssertPointer(y,3);
473: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
474: xf = fn->alpha*x;
475: PetscUseTypeMethod(fn,evaluatefunction,xf,&yf);
476: *y = fn->beta*yf;
477: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@
482: FNEvaluateDerivative - Computes the value of the derivative f'(x) for a given x.
484: Not Collective
486: Input Parameters:
487: + fn - the math function context
488: - x - the value where the derivative must be evaluated
490: Output Parameter:
491: . y - the result of f'(x)
493: Note:
494: Scaling factors are taken into account, so the actual derivative evaluation will
495: return alpha*beta*f'(alpha*x).
497: Level: intermediate
499: .seealso: FNEvaluateFunction()
500: @*/
501: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)
502: {
503: PetscScalar xf,yf;
505: PetscFunctionBegin;
508: PetscAssertPointer(y,3);
509: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
510: xf = fn->alpha*x;
511: PetscUseTypeMethod(fn,evaluatederivative,xf,&yf);
512: *y = fn->alpha*fn->beta*yf;
513: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,const PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)
518: {
519: PetscInt i,j;
520: PetscBLASInt n,k,ld,lwork,info;
521: PetscScalar *Q,*W,*work,adummy,a,x,y,one=1.0,zero=0.0;
522: PetscReal *eig,dummy;
523: #if defined(PETSC_USE_COMPLEX)
524: PetscReal *rwork,rdummy;
525: #endif
527: PetscFunctionBegin;
528: PetscCall(PetscBLASIntCast(m,&n));
529: ld = n;
530: k = firstonly? 1: n;
532: /* workspace query and memory allocation */
533: lwork = -1;
534: #if defined(PETSC_USE_COMPLEX)
535: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&rdummy,&info));
536: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
537: PetscCall(PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork));
538: #else
539: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&info));
540: PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
541: PetscCall(PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work));
542: #endif
544: /* compute eigendecomposition */
545: for (j=0;j<n;j++) for (i=j;i<n;i++) Q[i+j*ld] = As[i+j*ld];
546: #if defined(PETSC_USE_COMPLEX)
547: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
548: #else
549: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
550: #endif
551: SlepcCheckLapackInfo("syev",info);
553: /* W = f(Lambda)*Q' */
554: for (i=0;i<n;i++) {
555: x = fn->alpha*eig[i];
556: PetscUseTypeMethod(fn,evaluatefunction,x,&y); /* y = f(x) */
557: for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
558: }
559: /* Bs = Q*W */
560: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
561: #if defined(PETSC_USE_COMPLEX)
562: PetscCall(PetscFree5(eig,Q,W,work,rwork));
563: #else
564: PetscCall(PetscFree4(eig,Q,W,work));
565: #endif
566: PetscCall(PetscLogFlops(9.0*n*n*n+2.0*n*n*n));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: /*
571: FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
572: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
573: decomposition of A is A=Q*D*Q'
574: */
575: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)
576: {
577: PetscInt m;
578: const PetscScalar *As;
579: PetscScalar *Bs;
581: PetscFunctionBegin;
582: PetscCall(MatDenseGetArrayRead(A,&As));
583: PetscCall(MatDenseGetArray(B,&Bs));
584: PetscCall(MatGetSize(A,&m,NULL));
585: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE));
586: PetscCall(MatDenseRestoreArrayRead(A,&As));
587: PetscCall(MatDenseRestoreArray(B,&Bs));
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: static PetscErrorCode FNEvaluateFunctionMat_Basic(FN fn,Mat A,Mat F)
592: {
593: PetscBool iscuda;
595: PetscFunctionBegin;
596: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
597: if (iscuda && !fn->ops->evaluatefunctionmatcuda[fn->method]) PetscCall(PetscInfo(fn,"The method %" PetscInt_FMT " is not implemented for CUDA, falling back to CPU version\n",fn->method));
598: if (iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatcuda[fn->method],A,F);
599: else if (fn->ops->evaluatefunctionmat[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmat[fn->method],A,F);
600: else {
601: PetscCheck(fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
602: PetscCheck(!fn->method,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
603: }
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
608: {
609: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
610: PetscInt m,n;
611: PetscMPIInt size,rank,n2;
612: PetscScalar *pF;
613: Mat M,F;
615: PetscFunctionBegin;
616: /* destination matrix */
617: F = B?B:A;
619: /* check symmetry of A */
620: PetscCall(MatIsHermitianKnown(A,&set,&flg));
621: symm = set? flg: PETSC_FALSE;
623: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
624: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
625: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
626: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
627: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
628: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
629: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
630: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
631: PetscCall(FNEvaluateFunctionMat_Sym_Default(fn,A,F));
632: } else {
633: /* scale argument */
634: if (fn->alpha!=(PetscScalar)1.0) {
635: PetscCall(FN_AllocateWorkMat(fn,A,&M));
636: PetscCall(MatScale(M,fn->alpha));
637: } else M = A;
638: PetscCall(FNEvaluateFunctionMat_Basic(fn,M,F));
639: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
640: /* scale result */
641: PetscCall(MatScale(F,fn->beta));
642: }
643: PetscCall(PetscFPTrapPop());
644: }
645: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { /* synchronize */
646: PetscCall(MatGetSize(A,&m,&n));
647: PetscCall(MatDenseGetArray(F,&pF));
648: PetscCall(PetscMPIIntCast(n*n,&n2));
649: PetscCallMPI(MPI_Bcast(pF,n2,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
650: PetscCall(MatDenseRestoreArray(F,&pF));
651: }
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: /*@
656: FNEvaluateFunctionMat - Computes the value of the function f(A) for a given
657: matrix A, where the result is also a matrix.
659: Logically Collective
661: Input Parameters:
662: + fn - the math function context
663: - A - matrix on which the function must be evaluated
665: Output Parameter:
666: . B - (optional) matrix resulting from evaluating f(A)
668: Notes:
669: Matrix A must be a square sequential dense Mat, with all entries equal on
670: all processes (otherwise each process will compute different results).
671: If matrix B is provided, it must also be a square sequential dense Mat, and
672: both matrices must have the same dimensions. If B is NULL (or B=A) then the
673: function will perform an in-place computation, overwriting A with f(A).
675: If A is known to be real symmetric or complex Hermitian then it is
676: recommended to set the appropriate flag with MatSetOption(), because
677: symmetry can sometimes be exploited by the algorithm.
679: Scaling factors are taken into account, so the actual function evaluation
680: will return beta*f(alpha*A).
682: Level: advanced
684: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
685: @*/
686: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
687: {
688: PetscBool inplace=PETSC_FALSE;
689: PetscInt m,n,n1;
690: MatType type;
692: PetscFunctionBegin;
697: if (B) {
700: } else inplace = PETSC_TRUE;
701: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA); //SlepcMatCheckSeq(A);
702: PetscCall(MatGetSize(A,&m,&n));
703: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
704: if (!inplace) {
705: PetscCall(MatGetType(A,&type));
706: PetscCheckTypeName(B,type);
707: n1 = n;
708: PetscCall(MatGetSize(B,&m,&n));
709: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat B is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
710: PetscCheck(n1==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrices A and B must have the same dimension");
711: }
713: /* evaluate matrix function */
714: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
715: PetscCall(FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE));
716: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: /*
721: FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
722: and then copies the first column.
723: */
724: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
725: {
726: Mat F;
728: PetscFunctionBegin;
729: PetscCall(FN_AllocateWorkMat(fn,A,&F));
730: PetscCall(FNEvaluateFunctionMat_Basic(fn,A,F));
731: PetscCall(MatGetColumnVector(F,v,0));
732: PetscCall(FN_FreeWorkMat(fn,&F));
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: /*
737: FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
738: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
739: decomposition of A is A=Q*D*Q'. Only the first column is computed.
740: */
741: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
742: {
743: PetscInt m;
744: const PetscScalar *As;
745: PetscScalar *vs;
747: PetscFunctionBegin;
748: PetscCall(MatDenseGetArrayRead(A,&As));
749: PetscCall(VecGetArray(v,&vs));
750: PetscCall(MatGetSize(A,&m,NULL));
751: PetscCall(FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE));
752: PetscCall(MatDenseRestoreArrayRead(A,&As));
753: PetscCall(VecRestoreArray(v,&vs));
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
758: {
759: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
760: PetscInt m,n;
761: Mat M;
762: PetscMPIInt size,rank,n_;
763: PetscScalar *pv;
765: PetscFunctionBegin;
766: /* check symmetry of A */
767: PetscCall(MatIsHermitianKnown(A,&set,&flg));
768: symm = set? flg: PETSC_FALSE;
770: /* evaluate matrix function */
771: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
772: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank));
773: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
774: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
775: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
776: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
777: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
778: PetscCall(PetscInfo(fn,"Computing matrix function via diagonalization\n"));
779: PetscCall(FNEvaluateFunctionMatVec_Sym_Default(fn,A,v));
780: } else {
781: /* scale argument */
782: if (fn->alpha!=(PetscScalar)1.0) {
783: PetscCall(FN_AllocateWorkMat(fn,A,&M));
784: PetscCall(MatScale(M,fn->alpha));
785: } else M = A;
786: if (iscuda && fn->ops->evaluatefunctionmatveccuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatveccuda[fn->method],M,v);
787: else if (fn->ops->evaluatefunctionmatvec[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatvec[fn->method],M,v);
788: else PetscCall(FNEvaluateFunctionMatVec_Default(fn,M,v));
789: if (fn->alpha!=(PetscScalar)1.0) PetscCall(FN_FreeWorkMat(fn,&M));
790: /* scale result */
791: PetscCall(VecScale(v,fn->beta));
792: }
793: PetscCall(PetscFPTrapPop());
794: }
796: /* synchronize */
797: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
798: PetscCall(MatGetSize(A,&m,&n));
799: PetscCall(VecGetArray(v,&pv));
800: PetscCall(PetscMPIIntCast(n,&n_));
801: PetscCallMPI(MPI_Bcast(pv,n_,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn)));
802: PetscCall(VecRestoreArray(v,&pv));
803: }
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: /*@
808: FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
809: for a given matrix A.
811: Logically Collective
813: Input Parameters:
814: + fn - the math function context
815: - A - matrix on which the function must be evaluated
817: Output Parameter:
818: . v - vector to hold the first column of f(A)
820: Notes:
821: This operation is similar to FNEvaluateFunctionMat() but returns only
822: the first column of f(A), hence saving computations in most cases.
824: Level: advanced
826: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
827: @*/
828: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
829: {
830: PetscInt m,n;
831: PetscBool iscuda;
833: PetscFunctionBegin;
840: PetscCheckTypeNames(A,MATSEQDENSE,MATSEQDENSECUDA); //SlepcMatCheckSeq(A);
841: PetscCall(MatGetSize(A,&m,&n));
842: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
843: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
844: PetscCheckTypeName(v,iscuda?VECSEQCUDA:VECSEQ);
845: PetscCall(VecGetSize(v,&m));
846: PetscCheck(m==n,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrix A and vector v must have the same size");
847: PetscCall(PetscLogEventBegin(FN_Evaluate,fn,0,0,0));
848: PetscCall(FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE));
849: PetscCall(PetscLogEventEnd(FN_Evaluate,fn,0,0,0));
850: PetscFunctionReturn(PETSC_SUCCESS);
851: }
853: /*@
854: FNSetFromOptions - Sets FN options from the options database.
856: Collective
858: Input Parameters:
859: . fn - the math function context
861: Notes:
862: To see all options, run your program with the -help option.
864: Level: beginner
866: .seealso: FNSetOptionsPrefix()
867: @*/
868: PetscErrorCode FNSetFromOptions(FN fn)
869: {
870: char type[256];
871: PetscScalar array[2];
872: PetscInt k,meth;
873: PetscBool flg;
874: FNParallelType pmode;
876: PetscFunctionBegin;
878: PetscCall(FNRegisterAll());
879: PetscObjectOptionsBegin((PetscObject)fn);
880: PetscCall(PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg));
881: if (flg) PetscCall(FNSetType(fn,type));
882: else if (!((PetscObject)fn)->type_name) PetscCall(FNSetType(fn,FNRATIONAL));
884: k = 2;
885: array[0] = 0.0; array[1] = 0.0;
886: PetscCall(PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg));
887: if (flg) {
888: if (k<2) array[1] = 1.0;
889: PetscCall(FNSetScale(fn,array[0],array[1]));
890: }
892: PetscCall(PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg));
893: if (flg) PetscCall(FNSetMethod(fn,meth));
895: PetscCall(PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg));
896: if (flg) PetscCall(FNSetParallel(fn,pmode));
898: PetscTryTypeMethod(fn,setfromoptions,PetscOptionsObject);
899: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fn,PetscOptionsObject));
900: PetscOptionsEnd();
901: PetscFunctionReturn(PETSC_SUCCESS);
902: }
904: /*@
905: FNView - Prints the FN data structure.
907: Collective
909: Input Parameters:
910: + fn - the math function context
911: - viewer - optional visualization context
913: Note:
914: The available visualization contexts include
915: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
916: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
917: output where only the first processor opens
918: the file. All other processors send their
919: data to the first processor to print.
921: The user can open an alternative visualization context with
922: PetscViewerASCIIOpen() - output to a specified file.
924: Level: beginner
926: .seealso: FNCreate()
927: @*/
928: PetscErrorCode FNView(FN fn,PetscViewer viewer)
929: {
930: PetscBool isascii;
931: PetscMPIInt size;
933: PetscFunctionBegin;
935: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer));
937: PetscCheckSameComm(fn,1,viewer,2);
938: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
939: if (isascii) {
940: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer));
941: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size));
942: if (size>1) PetscCall(PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",FNParallelTypes[fn->pmode]));
943: PetscCall(PetscViewerASCIIPushTab(viewer));
944: PetscTryTypeMethod(fn,view,viewer);
945: PetscCall(PetscViewerASCIIPopTab(viewer));
946: }
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
950: /*@
951: FNViewFromOptions - View from options
953: Collective
955: Input Parameters:
956: + fn - the math function context
957: . obj - optional object
958: - name - command line option
960: Level: intermediate
962: .seealso: FNView(), FNCreate()
963: @*/
964: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])
965: {
966: PetscFunctionBegin;
968: PetscCall(PetscObjectViewFromOptions((PetscObject)fn,obj,name));
969: PetscFunctionReturn(PETSC_SUCCESS);
970: }
972: /*@
973: FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
974: different communicator.
976: Collective
978: Input Parameters:
979: + fn - the math function context
980: - comm - MPI communicator
982: Output Parameter:
983: . newfn - location to put the new FN context
985: Note:
986: In order to use the same MPI communicator as in the original object,
987: use PetscObjectComm((PetscObject)fn).
989: Level: developer
991: .seealso: FNCreate()
992: @*/
993: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)
994: {
995: FNType type;
996: PetscScalar alpha,beta;
997: PetscInt meth;
998: FNParallelType ptype;
1000: PetscFunctionBegin;
1003: PetscAssertPointer(newfn,3);
1004: PetscCall(FNCreate(comm,newfn));
1005: PetscCall(FNGetType(fn,&type));
1006: PetscCall(FNSetType(*newfn,type));
1007: PetscCall(FNGetScale(fn,&alpha,&beta));
1008: PetscCall(FNSetScale(*newfn,alpha,beta));
1009: PetscCall(FNGetMethod(fn,&meth));
1010: PetscCall(FNSetMethod(*newfn,meth));
1011: PetscCall(FNGetParallel(fn,&ptype));
1012: PetscCall(FNSetParallel(*newfn,ptype));
1013: PetscTryTypeMethod(fn,duplicate,comm,newfn);
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: }
1017: /*@
1018: FNDestroy - Destroys FN context that was created with FNCreate().
1020: Collective
1022: Input Parameter:
1023: . fn - the math function context
1025: Level: beginner
1027: .seealso: FNCreate()
1028: @*/
1029: PetscErrorCode FNDestroy(FN *fn)
1030: {
1031: PetscInt i;
1033: PetscFunctionBegin;
1034: if (!*fn) PetscFunctionReturn(PETSC_SUCCESS);
1036: if (--((PetscObject)*fn)->refct > 0) { *fn = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
1037: PetscTryTypeMethod(*fn,destroy);
1038: for (i=0;i<(*fn)->nw;i++) PetscCall(MatDestroy(&(*fn)->W[i]));
1039: PetscCall(PetscHeaderDestroy(fn));
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: /*@C
1044: FNRegister - Adds a mathematical function to the FN package.
1046: Not Collective
1048: Input Parameters:
1049: + name - name of a new user-defined FN
1050: - function - routine to create context
1052: Notes:
1053: FNRegister() may be called multiple times to add several user-defined functions.
1055: Level: advanced
1057: .seealso: FNRegisterAll()
1058: @*/
1059: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1060: {
1061: PetscFunctionBegin;
1062: PetscCall(FNInitializePackage());
1063: PetscCall(PetscFunctionListAdd(&FNList,name,function));
1064: PetscFunctionReturn(PETSC_SUCCESS);
1065: }