Actual source code: fnbasic.c
slepc-3.18.3 2023-03-24
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: PetscFunctionListDestroy(&FNList);
36: FNPackageInitialized = PETSC_FALSE;
37: FNRegisterAllCalled = PETSC_FALSE;
38: return 0;
39: }
41: /*@C
42: FNInitializePackage - This function initializes everything in the FN package.
43: It is called from PetscDLLibraryRegister() when using dynamic libraries, and
44: on the first call to FNCreate() when using static libraries.
46: Level: developer
48: .seealso: SlepcInitialize()
49: @*/
50: PetscErrorCode FNInitializePackage(void)
51: {
52: char logList[256];
53: PetscBool opt,pkg;
54: PetscClassId classids[1];
56: if (FNPackageInitialized) return 0;
57: FNPackageInitialized = PETSC_TRUE;
58: /* Register Classes */
59: PetscClassIdRegister("Math Function",&FN_CLASSID);
60: /* Register Constructors */
61: FNRegisterAll();
62: /* Register Events */
63: PetscLogEventRegister("FNEvaluate",FN_CLASSID,&FN_Evaluate);
64: /* Process Info */
65: classids[0] = FN_CLASSID;
66: PetscInfoProcessClass("fn",1,&classids[0]);
67: /* Process summary exclusions */
68: PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
69: if (opt) {
70: PetscStrInList("fn",logList,',',&pkg);
71: if (pkg) PetscLogEventDeactivateClass(FN_CLASSID);
72: }
73: /* Register package finalizer */
74: PetscRegisterFinalize(FNFinalizePackage);
75: return 0;
76: }
78: /*@
79: FNCreate - Creates an FN context.
81: Collective
83: Input Parameter:
84: . comm - MPI communicator
86: Output Parameter:
87: . newfn - location to put the FN context
89: Level: beginner
91: .seealso: FNDestroy(), FN
92: @*/
93: PetscErrorCode FNCreate(MPI_Comm comm,FN *newfn)
94: {
95: FN fn;
98: *newfn = NULL;
99: FNInitializePackage();
100: SlepcHeaderCreate(fn,FN_CLASSID,"FN","Math Function","FN",comm,FNDestroy,FNView);
102: fn->alpha = 1.0;
103: fn->beta = 1.0;
104: fn->method = 0;
106: fn->nw = 0;
107: fn->cw = 0;
108: fn->data = NULL;
110: *newfn = fn;
111: return 0;
112: }
114: /*@C
115: FNSetOptionsPrefix - Sets the prefix used for searching for all
116: FN options in the database.
118: Logically Collective on fn
120: Input Parameters:
121: + fn - the math function context
122: - prefix - the prefix string to prepend to all FN option requests
124: Notes:
125: A hyphen (-) must NOT be given at the beginning of the prefix name.
126: The first character of all runtime options is AUTOMATICALLY the
127: hyphen.
129: Level: advanced
131: .seealso: FNAppendOptionsPrefix()
132: @*/
133: PetscErrorCode FNSetOptionsPrefix(FN fn,const char *prefix)
134: {
136: PetscObjectSetOptionsPrefix((PetscObject)fn,prefix);
137: return 0;
138: }
140: /*@C
141: FNAppendOptionsPrefix - Appends to the prefix used for searching for all
142: FN options in the database.
144: Logically Collective on fn
146: Input Parameters:
147: + fn - the math function context
148: - prefix - the prefix string to prepend to all FN option requests
150: Notes:
151: A hyphen (-) must NOT be given at the beginning of the prefix name.
152: The first character of all runtime options is AUTOMATICALLY the hyphen.
154: Level: advanced
156: .seealso: FNSetOptionsPrefix()
157: @*/
158: PetscErrorCode FNAppendOptionsPrefix(FN fn,const char *prefix)
159: {
161: PetscObjectAppendOptionsPrefix((PetscObject)fn,prefix);
162: return 0;
163: }
165: /*@C
166: FNGetOptionsPrefix - Gets the prefix used for searching for all
167: FN options in the database.
169: Not Collective
171: Input Parameters:
172: . fn - the math function context
174: Output Parameters:
175: . prefix - pointer to the prefix string used is returned
177: Note:
178: On the Fortran side, the user should pass in a string 'prefix' of
179: sufficient length to hold the prefix.
181: Level: advanced
183: .seealso: FNSetOptionsPrefix(), FNAppendOptionsPrefix()
184: @*/
185: PetscErrorCode FNGetOptionsPrefix(FN fn,const char *prefix[])
186: {
189: PetscObjectGetOptionsPrefix((PetscObject)fn,prefix);
190: return 0;
191: }
193: /*@C
194: FNSetType - Selects the type for the FN object.
196: Logically Collective on fn
198: Input Parameters:
199: + fn - the math function context
200: - type - a known type
202: Notes:
203: The default is FNRATIONAL, which includes polynomials as a particular
204: case as well as simple functions such as f(x)=x and f(x)=constant.
206: Level: intermediate
208: .seealso: FNGetType()
209: @*/
210: PetscErrorCode FNSetType(FN fn,FNType type)
211: {
212: PetscErrorCode (*r)(FN);
213: PetscBool match;
218: PetscObjectTypeCompare((PetscObject)fn,type,&match);
219: if (match) return 0;
221: PetscFunctionListFind(FNList,type,&r);
224: PetscTryTypeMethod(fn,destroy);
225: PetscMemzero(fn->ops,sizeof(struct _FNOps));
227: PetscObjectChangeTypeName((PetscObject)fn,type);
228: (*r)(fn);
229: return 0;
230: }
232: /*@C
233: FNGetType - Gets the FN type name (as a string) from the FN context.
235: Not Collective
237: Input Parameter:
238: . fn - the math function context
240: Output Parameter:
241: . type - name of the math function
243: Level: intermediate
245: .seealso: FNSetType()
246: @*/
247: PetscErrorCode FNGetType(FN fn,FNType *type)
248: {
251: *type = ((PetscObject)fn)->type_name;
252: return 0;
253: }
255: /*@
256: FNSetScale - Sets the scaling parameters that define the matematical function.
258: Logically Collective on fn
260: Input Parameters:
261: + fn - the math function context
262: . alpha - inner scaling (argument)
263: - beta - outer scaling (result)
265: Notes:
266: Given a function f(x) specified by the FN type, the scaling parameters can
267: be used to realize the function beta*f(alpha*x). So when these values are given,
268: the procedure for function evaluation will first multiply the argument by alpha,
269: then evaluate the function itself, and finally scale the result by beta.
270: Likewise, these values are also considered when evaluating the derivative.
272: If you want to provide only one of the two scaling factors, set the other
273: one to 1.0.
275: Level: intermediate
277: .seealso: FNGetScale(), FNEvaluateFunction()
278: @*/
279: PetscErrorCode FNSetScale(FN fn,PetscScalar alpha,PetscScalar beta)
280: {
285: fn->alpha = alpha;
286: fn->beta = beta;
287: return 0;
288: }
290: /*@
291: FNGetScale - Gets the scaling parameters that define the matematical function.
293: Not Collective
295: Input Parameter:
296: . fn - the math function context
298: Output Parameters:
299: + alpha - inner scaling (argument)
300: - beta - outer scaling (result)
302: Level: intermediate
304: .seealso: FNSetScale()
305: @*/
306: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)
307: {
309: if (alpha) *alpha = fn->alpha;
310: if (beta) *beta = fn->beta;
311: return 0;
312: }
314: /*@
315: FNSetMethod - Selects the method to be used to evaluate functions of matrices.
317: Logically Collective on fn
319: Input Parameters:
320: + fn - the math function context
321: - meth - an index identifying the method
323: Options Database Key:
324: . -fn_method <meth> - Sets the method
326: Notes:
327: In some FN types there are more than one algorithms available for computing
328: matrix functions. In that case, this function allows choosing the wanted method.
330: If meth is currently set to 0 (the default) and the input argument A of
331: FNEvaluateFunctionMat() is a symmetric/Hermitian matrix, then the computation
332: is done via the eigendecomposition of A, rather than with the general algorithm.
334: Level: intermediate
336: .seealso: FNGetMethod(), FNEvaluateFunctionMat()
337: @*/
338: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)
339: {
344: fn->method = meth;
345: return 0;
346: }
348: /*@
349: FNGetMethod - Gets the method currently used in the FN.
351: Not Collective
353: Input Parameter:
354: . fn - the math function context
356: Output Parameter:
357: . meth - identifier of the method
359: Level: intermediate
361: .seealso: FNSetMethod()
362: @*/
363: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)
364: {
367: *meth = fn->method;
368: return 0;
369: }
371: /*@
372: FNSetParallel - Selects the mode of operation in parallel runs.
374: Logically Collective on fn
376: Input Parameters:
377: + fn - the math function context
378: - pmode - the parallel mode
380: Options Database Key:
381: . -fn_parallel <mode> - Sets the parallel mode, either 'redundant' or 'synchronized'
383: Notes:
384: This is relevant only when the function is evaluated on a matrix, with
385: either FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
387: In the 'redundant' parallel mode, all processes will make the computation
388: redundantly, starting from the same data, and producing the same result.
389: This result may be slightly different in the different processes if using a
390: multithreaded BLAS library, which may cause issues in ill-conditioned problems.
392: In the 'synchronized' parallel mode, only the first MPI process performs the
393: computation and then the computed matrix is broadcast to the other
394: processes in the communicator. This communication is done automatically at
395: the end of FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
397: Level: advanced
399: .seealso: FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec(), FNGetParallel()
400: @*/
401: PetscErrorCode FNSetParallel(FN fn,FNParallelType pmode)
402: {
405: fn->pmode = pmode;
406: return 0;
407: }
409: /*@
410: FNGetParallel - Gets the mode of operation in parallel runs.
412: Not Collective
414: Input Parameter:
415: . fn - the math function context
417: Output Parameter:
418: . pmode - the parallel mode
420: Level: advanced
422: .seealso: FNSetParallel()
423: @*/
424: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)
425: {
428: *pmode = fn->pmode;
429: return 0;
430: }
432: /*@
433: FNEvaluateFunction - Computes the value of the function f(x) for a given x.
435: Not collective
437: Input Parameters:
438: + fn - the math function context
439: - x - the value where the function must be evaluated
441: Output Parameter:
442: . y - the result of f(x)
444: Note:
445: Scaling factors are taken into account, so the actual function evaluation
446: will return beta*f(alpha*x).
448: Level: intermediate
450: .seealso: FNEvaluateDerivative(), FNEvaluateFunctionMat(), FNSetScale()
451: @*/
452: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)
453: {
454: PetscScalar xf,yf;
459: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
460: xf = fn->alpha*x;
461: PetscUseTypeMethod(fn,evaluatefunction,xf,&yf);
462: *y = fn->beta*yf;
463: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
464: return 0;
465: }
467: /*@
468: FNEvaluateDerivative - Computes the value of the derivative f'(x) for a given x.
470: Not Collective
472: Input Parameters:
473: + fn - the math function context
474: - x - the value where the derivative must be evaluated
476: Output Parameter:
477: . y - the result of f'(x)
479: Note:
480: Scaling factors are taken into account, so the actual derivative evaluation will
481: return alpha*beta*f'(alpha*x).
483: Level: intermediate
485: .seealso: FNEvaluateFunction()
486: @*/
487: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)
488: {
489: PetscScalar xf,yf;
494: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
495: xf = fn->alpha*x;
496: PetscUseTypeMethod(fn,evaluatederivative,xf,&yf);
497: *y = fn->alpha*fn->beta*yf;
498: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
499: return 0;
500: }
502: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,const PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)
503: {
504: PetscInt i,j;
505: PetscBLASInt n,k,ld,lwork,info;
506: PetscScalar *Q,*W,*work,adummy,a,x,y,one=1.0,zero=0.0;
507: PetscReal *eig,dummy;
508: #if defined(PETSC_USE_COMPLEX)
509: PetscReal *rwork,rdummy;
510: #endif
512: PetscBLASIntCast(m,&n);
513: ld = n;
514: k = firstonly? 1: n;
516: /* workspace query and memory allocation */
517: lwork = -1;
518: #if defined(PETSC_USE_COMPLEX)
519: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&rdummy,&info));
520: PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
521: PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork);
522: #else
523: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,&adummy,&ld,&dummy,&a,&lwork,&info));
524: PetscBLASIntCast((PetscInt)a,&lwork);
525: PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work);
526: #endif
528: /* compute eigendecomposition */
529: for (j=0;j<n;j++) for (i=j;i<n;i++) Q[i+j*ld] = As[i+j*ld];
530: #if defined(PETSC_USE_COMPLEX)
531: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
532: #else
533: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
534: #endif
535: SlepcCheckLapackInfo("syev",info);
537: /* W = f(Lambda)*Q' */
538: for (i=0;i<n;i++) {
539: x = fn->alpha*eig[i];
540: PetscUseTypeMethod(fn,evaluatefunction,x,&y); /* y = f(x) */
541: for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
542: }
543: /* Bs = Q*W */
544: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
545: #if defined(PETSC_USE_COMPLEX)
546: PetscFree5(eig,Q,W,work,rwork);
547: #else
548: PetscFree4(eig,Q,W,work);
549: #endif
550: PetscLogFlops(9.0*n*n*n+2.0*n*n*n);
551: return 0;
552: }
554: /*
555: FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
556: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
557: decomposition of A is A=Q*D*Q'
558: */
559: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)
560: {
561: PetscInt m;
562: const PetscScalar *As;
563: PetscScalar *Bs;
565: MatDenseGetArrayRead(A,&As);
566: MatDenseGetArray(B,&Bs);
567: MatGetSize(A,&m,NULL);
568: FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE);
569: MatDenseRestoreArrayRead(A,&As);
570: MatDenseRestoreArray(B,&Bs);
571: return 0;
572: }
574: PetscErrorCode FNEvaluateFunctionMat_Basic(FN fn,Mat A,Mat F)
575: {
576: PetscBool iscuda;
578: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda);
579: if (iscuda && !fn->ops->evaluatefunctionmatcuda[fn->method]) PetscInfo(fn,"The method %" PetscInt_FMT " is not implemented for CUDA, falling back to CPU version\n",fn->method);
580: if (iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatcuda[fn->method],A,F);
581: else if (fn->ops->evaluatefunctionmat[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmat[fn->method],A,F);
582: else {
585: }
586: return 0;
587: }
589: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)
590: {
591: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
592: PetscInt m,n;
593: PetscMPIInt size,rank;
594: PetscScalar *pF;
595: Mat M,F;
597: /* destination matrix */
598: F = B?B:A;
600: /* check symmetry of A */
601: MatIsHermitianKnown(A,&set,&flg);
602: symm = set? flg: PETSC_FALSE;
604: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
605: MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
606: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
607: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
608: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda);
609: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
610: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
611: PetscInfo(fn,"Computing matrix function via diagonalization\n");
612: FNEvaluateFunctionMat_Sym_Default(fn,A,F);
613: } else {
614: /* scale argument */
615: if (fn->alpha!=(PetscScalar)1.0) {
616: FN_AllocateWorkMat(fn,A,&M);
617: MatScale(M,fn->alpha);
618: } else M = A;
619: FNEvaluateFunctionMat_Basic(fn,M,F);
620: if (fn->alpha!=(PetscScalar)1.0) FN_FreeWorkMat(fn,&M);
621: /* scale result */
622: MatScale(F,fn->beta);
623: }
624: PetscFPTrapPop();
625: }
626: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { /* synchronize */
627: MatGetSize(A,&m,&n);
628: MatDenseGetArray(F,&pF);
629: MPI_Bcast(pF,n*n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
630: MatDenseRestoreArray(F,&pF);
631: }
632: return 0;
633: }
635: /*@
636: FNEvaluateFunctionMat - Computes the value of the function f(A) for a given
637: matrix A, where the result is also a matrix.
639: Logically Collective on fn
641: Input Parameters:
642: + fn - the math function context
643: - A - matrix on which the function must be evaluated
645: Output Parameter:
646: . B - (optional) matrix resulting from evaluating f(A)
648: Notes:
649: Matrix A must be a square sequential dense Mat, with all entries equal on
650: all processes (otherwise each process will compute different results).
651: If matrix B is provided, it must also be a square sequential dense Mat, and
652: both matrices must have the same dimensions. If B is NULL (or B=A) then the
653: function will perform an in-place computation, overwriting A with f(A).
655: If A is known to be real symmetric or complex Hermitian then it is
656: recommended to set the appropriate flag with MatSetOption(), because
657: symmetry can sometimes be exploited by the algorithm.
659: Scaling factors are taken into account, so the actual function evaluation
660: will return beta*f(alpha*A).
662: Level: advanced
664: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
665: @*/
666: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)
667: {
668: PetscBool inplace=PETSC_FALSE;
669: PetscInt m,n,n1;
670: MatType type;
676: if (B) {
679: } else inplace = PETSC_TRUE;
681: MatGetSize(A,&m,&n);
683: if (!inplace) {
684: MatGetType(A,&type);
686: n1 = n;
687: MatGetSize(B,&m,&n);
690: }
692: /* evaluate matrix function */
693: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
694: FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE);
695: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
696: return 0;
697: }
699: /*
700: FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
701: and then copies the first column.
702: */
703: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)
704: {
705: Mat F;
707: FN_AllocateWorkMat(fn,A,&F);
708: FNEvaluateFunctionMat_Basic(fn,A,F);
709: MatGetColumnVector(F,v,0);
710: FN_FreeWorkMat(fn,&F);
711: return 0;
712: }
714: /*
715: FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
716: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
717: decomposition of A is A=Q*D*Q'. Only the first column is computed.
718: */
719: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)
720: {
721: PetscInt m;
722: const PetscScalar *As;
723: PetscScalar *vs;
725: MatDenseGetArrayRead(A,&As);
726: VecGetArray(v,&vs);
727: MatGetSize(A,&m,NULL);
728: FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE);
729: MatDenseRestoreArrayRead(A,&As);
730: VecRestoreArray(v,&vs);
731: return 0;
732: }
734: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)
735: {
736: PetscBool set,flg,symm=PETSC_FALSE,iscuda,hasspecificmeth;
737: PetscInt m,n;
738: Mat M;
739: PetscMPIInt size,rank;
740: PetscScalar *pv;
742: /* check symmetry of A */
743: MatIsHermitianKnown(A,&set,&flg);
744: symm = set? flg: PETSC_FALSE;
746: /* evaluate matrix function */
747: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
748: MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
749: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
750: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
751: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda);
752: hasspecificmeth = ((iscuda && fn->ops->evaluatefunctionmatcuda[fn->method]) || (!iscuda && fn->method && fn->ops->evaluatefunctionmat[fn->method]))? PETSC_TRUE: PETSC_FALSE;
753: if (!hasspecificmeth && symm && !fn->method) { /* prefer diagonalization */
754: PetscInfo(fn,"Computing matrix function via diagonalization\n");
755: FNEvaluateFunctionMatVec_Sym_Default(fn,A,v);
756: } else {
757: /* scale argument */
758: if (fn->alpha!=(PetscScalar)1.0) {
759: FN_AllocateWorkMat(fn,A,&M);
760: MatScale(M,fn->alpha);
761: } else M = A;
762: if (iscuda && fn->ops->evaluatefunctionmatveccuda[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatveccuda[fn->method],M,v);
763: else if (fn->ops->evaluatefunctionmatvec[fn->method]) PetscUseTypeMethod(fn,evaluatefunctionmatvec[fn->method],M,v);
764: else FNEvaluateFunctionMatVec_Default(fn,M,v);
765: if (fn->alpha!=(PetscScalar)1.0) FN_FreeWorkMat(fn,&M);
766: /* scale result */
767: VecScale(v,fn->beta);
768: }
769: PetscFPTrapPop();
770: }
772: /* synchronize */
773: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
774: MatGetSize(A,&m,&n);
775: VecGetArray(v,&pv);
776: MPI_Bcast(pv,n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
777: VecRestoreArray(v,&pv);
778: }
779: return 0;
780: }
782: /*@
783: FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
784: for a given matrix A.
786: Logically Collective on fn
788: Input Parameters:
789: + fn - the math function context
790: - A - matrix on which the function must be evaluated
792: Output Parameter:
793: . v - vector to hold the first column of f(A)
795: Notes:
796: This operation is similar to FNEvaluateFunctionMat() but returns only
797: the first column of f(A), hence saving computations in most cases.
799: Level: advanced
801: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
802: @*/
803: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)
804: {
805: PetscInt m,n;
806: PetscBool iscuda;
815: MatGetSize(A,&m,&n);
817: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda);
819: VecGetSize(v,&m);
821: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
822: FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE);
823: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
824: return 0;
825: }
827: /*@
828: FNSetFromOptions - Sets FN options from the options database.
830: Collective on fn
832: Input Parameters:
833: . fn - the math function context
835: Notes:
836: To see all options, run your program with the -help option.
838: Level: beginner
840: .seealso: FNSetOptionsPrefix()
841: @*/
842: PetscErrorCode FNSetFromOptions(FN fn)
843: {
844: char type[256];
845: PetscScalar array[2];
846: PetscInt k,meth;
847: PetscBool flg;
848: FNParallelType pmode;
851: FNRegisterAll();
852: PetscObjectOptionsBegin((PetscObject)fn);
853: PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg);
854: if (flg) FNSetType(fn,type);
855: else if (!((PetscObject)fn)->type_name) FNSetType(fn,FNRATIONAL);
857: k = 2;
858: array[0] = 0.0; array[1] = 0.0;
859: PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg);
860: if (flg) {
861: if (k<2) array[1] = 1.0;
862: FNSetScale(fn,array[0],array[1]);
863: }
865: PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg);
866: if (flg) FNSetMethod(fn,meth);
868: PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg);
869: if (flg) FNSetParallel(fn,pmode);
871: PetscTryTypeMethod(fn,setfromoptions,PetscOptionsObject);
872: PetscObjectProcessOptionsHandlers((PetscObject)fn,PetscOptionsObject);
873: PetscOptionsEnd();
874: return 0;
875: }
877: /*@C
878: FNView - Prints the FN data structure.
880: Collective on fn
882: Input Parameters:
883: + fn - the math function context
884: - viewer - optional visualization context
886: Note:
887: The available visualization contexts include
888: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
889: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
890: output where only the first processor opens
891: the file. All other processors send their
892: data to the first processor to print.
894: The user can open an alternative visualization context with
895: PetscViewerASCIIOpen() - output to a specified file.
897: Level: beginner
899: .seealso: FNCreate()
900: @*/
901: PetscErrorCode FNView(FN fn,PetscViewer viewer)
902: {
903: PetscBool isascii;
904: PetscMPIInt size;
907: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer);
910: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
911: if (isascii) {
912: PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer);
913: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
914: if (size>1) PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",FNParallelTypes[fn->pmode]);
915: PetscViewerASCIIPushTab(viewer);
916: PetscTryTypeMethod(fn,view,viewer);
917: PetscViewerASCIIPopTab(viewer);
918: }
919: return 0;
920: }
922: /*@C
923: FNViewFromOptions - View from options
925: Collective on FN
927: Input Parameters:
928: + fn - the math function context
929: . obj - optional object
930: - name - command line option
932: Level: intermediate
934: .seealso: FNView(), FNCreate()
935: @*/
936: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])
937: {
939: PetscObjectViewFromOptions((PetscObject)fn,obj,name);
940: return 0;
941: }
943: /*@
944: FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
945: different communicator.
947: Collective on fn
949: Input Parameters:
950: + fn - the math function context
951: - comm - MPI communicator
953: Output Parameter:
954: . newfn - location to put the new FN context
956: Note:
957: In order to use the same MPI communicator as in the original object,
958: use PetscObjectComm((PetscObject)fn).
960: Level: developer
962: .seealso: FNCreate()
963: @*/
964: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)
965: {
966: FNType type;
967: PetscScalar alpha,beta;
968: PetscInt meth;
969: FNParallelType ptype;
974: FNCreate(comm,newfn);
975: FNGetType(fn,&type);
976: FNSetType(*newfn,type);
977: FNGetScale(fn,&alpha,&beta);
978: FNSetScale(*newfn,alpha,beta);
979: FNGetMethod(fn,&meth);
980: FNSetMethod(*newfn,meth);
981: FNGetParallel(fn,&ptype);
982: FNSetParallel(*newfn,ptype);
983: PetscTryTypeMethod(fn,duplicate,comm,newfn);
984: return 0;
985: }
987: /*@C
988: FNDestroy - Destroys FN context that was created with FNCreate().
990: Collective on fn
992: Input Parameter:
993: . fn - the math function context
995: Level: beginner
997: .seealso: FNCreate()
998: @*/
999: PetscErrorCode FNDestroy(FN *fn)
1000: {
1001: PetscInt i;
1003: if (!*fn) return 0;
1005: if (--((PetscObject)(*fn))->refct > 0) { *fn = NULL; return 0; }
1006: PetscTryTypeMethod(*fn,destroy);
1007: for (i=0;i<(*fn)->nw;i++) MatDestroy(&(*fn)->W[i]);
1008: PetscHeaderDestroy(fn);
1009: return 0;
1010: }
1012: /*@C
1013: FNRegister - Adds a mathematical function to the FN package.
1015: Not collective
1017: Input Parameters:
1018: + name - name of a new user-defined FN
1019: - function - routine to create context
1021: Notes:
1022: FNRegister() may be called multiple times to add several user-defined functions.
1024: Level: advanced
1026: .seealso: FNRegisterAll()
1027: @*/
1028: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))
1029: {
1030: FNInitializePackage();
1031: PetscFunctionListAdd(&FNList,name,function);
1032: return 0;
1033: }