1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcfn.h" I*/
15: #include <slepcblaslapack.h>
17: PetscFunctionList FNList = 0;
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_",0};
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: {
38: PetscFunctionListDestroy(&FNList);
39: FNPackageInitialized = PETSC_FALSE;
40: FNRegisterAllCalled = PETSC_FALSE;
41: return(0);
42: }
44: /*@C
45: FNInitializePackage - This function initializes everything in the FN package.
46: It is called from PetscDLLibraryRegister() when using dynamic libraries, and
47: on the first call to FNCreate() when using static libraries.
49: Level: developer
51: .seealso: SlepcInitialize()
52: @*/
53: PetscErrorCode FNInitializePackage(void) 54: {
55: char logList[256];
56: PetscBool opt,pkg;
60: if (FNPackageInitialized) return(0);
61: FNPackageInitialized = PETSC_TRUE;
62: /* Register Classes */
63: PetscClassIdRegister("Math Function",&FN_CLASSID);
64: /* Register Constructors */
65: FNRegisterAll();
66: /* Register Events */
67: PetscLogEventRegister("FNEvaluate",FN_CLASSID,&FN_Evaluate);
68: /* Process info exclusions */
69: PetscOptionsGetString(NULL,NULL,"-info_exclude",logList,sizeof(logList),&opt);
70: if (opt) {
71: PetscStrInList("fn",logList,',',&pkg);
72: if (pkg) { PetscInfoDeactivateClass(FN_CLASSID); }
73: }
74: /* Process summary exclusions */
75: PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
76: if (opt) {
77: PetscStrInList("fn",logList,',',&pkg);
78: if (pkg) { PetscLogEventDeactivateClass(FN_CLASSID); }
79: }
80: /* Register package finalizer */
81: PetscRegisterFinalize(FNFinalizePackage);
82: return(0);
83: }
85: /*@
86: FNCreate - Creates an FN context.
88: Collective on MPI_Comm
90: Input Parameter:
91: . comm - MPI communicator
93: Output Parameter:
94: . newfn - location to put the FN context
96: Level: beginner
98: .seealso: FNDestroy(), FN 99: @*/
100: PetscErrorCode FNCreate(MPI_Comm comm,FN *newfn)101: {
102: FN fn;
107: *newfn = 0;
108: FNInitializePackage();
109: SlepcHeaderCreate(fn,FN_CLASSID,"FN","Math Function","FN",comm,FNDestroy,FNView);
111: fn->alpha = 1.0;
112: fn->beta = 1.0;
113: fn->method = 0;
115: fn->nw = 0;
116: fn->cw = 0;
117: fn->data = NULL;
119: *newfn = fn;
120: return(0);
121: }
123: /*@C
124: FNSetOptionsPrefix - Sets the prefix used for searching for all
125: FN options in the database.
127: Logically Collective on FN129: Input Parameters:
130: + fn - the math function context
131: - prefix - the prefix string to prepend to all FN option requests
133: Notes:
134: A hyphen (-) must NOT be given at the beginning of the prefix name.
135: The first character of all runtime options is AUTOMATICALLY the
136: hyphen.
138: Level: advanced
140: .seealso: FNAppendOptionsPrefix()
141: @*/
142: PetscErrorCode FNSetOptionsPrefix(FN fn,const char *prefix)143: {
148: PetscObjectSetOptionsPrefix((PetscObject)fn,prefix);
149: return(0);
150: }
152: /*@C
153: FNAppendOptionsPrefix - Appends to the prefix used for searching for all
154: FN options in the database.
156: Logically Collective on FN158: Input Parameters:
159: + fn - the math function context
160: - prefix - the prefix string to prepend to all FN option requests
162: Notes:
163: A hyphen (-) must NOT be given at the beginning of the prefix name.
164: The first character of all runtime options is AUTOMATICALLY the hyphen.
166: Level: advanced
168: .seealso: FNSetOptionsPrefix()
169: @*/
170: PetscErrorCode FNAppendOptionsPrefix(FN fn,const char *prefix)171: {
176: PetscObjectAppendOptionsPrefix((PetscObject)fn,prefix);
177: return(0);
178: }
180: /*@C
181: FNGetOptionsPrefix - Gets the prefix used for searching for all
182: FN options in the database.
184: Not Collective
186: Input Parameters:
187: . fn - the math function context
189: Output Parameters:
190: . prefix - pointer to the prefix string used is returned
192: Note:
193: On the Fortran side, the user should pass in a string 'prefix' of
194: sufficient length to hold the prefix.
196: Level: advanced
198: .seealso: FNSetOptionsPrefix(), FNAppendOptionsPrefix()
199: @*/
200: PetscErrorCode FNGetOptionsPrefix(FN fn,const char *prefix[])201: {
207: PetscObjectGetOptionsPrefix((PetscObject)fn,prefix);
208: return(0);
209: }
211: /*@C
212: FNSetType - Selects the type for the FN object.
214: Logically Collective on FN216: Input Parameter:
217: + fn - the math function context
218: - type - a known type
220: Notes:
221: The default is FNRATIONAL, which includes polynomials as a particular
222: case as well as simple functions such as f(x)=x and f(x)=constant.
224: Level: intermediate
226: .seealso: FNGetType()
227: @*/
228: PetscErrorCode FNSetType(FN fn,FNType type)229: {
230: PetscErrorCode ierr,(*r)(FN);
231: PetscBool match;
237: PetscObjectTypeCompare((PetscObject)fn,type,&match);
238: if (match) return(0);
240: PetscFunctionListFind(FNList,type,&r);
241: if (!r) SETERRQ1(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested FN type %s",type);
243: if (fn->ops->destroy) { (*fn->ops->destroy)(fn); }
244: PetscMemzero(fn->ops,sizeof(struct _FNOps));
246: PetscObjectChangeTypeName((PetscObject)fn,type);
247: (*r)(fn);
248: return(0);
249: }
251: /*@C
252: FNGetType - Gets the FN type name (as a string) from the FN context.
254: Not Collective
256: Input Parameter:
257: . fn - the math function context
259: Output Parameter:
260: . name - name of the math function
262: Level: intermediate
264: .seealso: FNSetType()
265: @*/
266: PetscErrorCode FNGetType(FN fn,FNType *type)267: {
271: *type = ((PetscObject)fn)->type_name;
272: return(0);
273: }
275: /*@
276: FNSetScale - Sets the scaling parameters that define the matematical function.
278: Logically Collective on FN280: Input Parameters:
281: + fn - the math function context
282: . alpha - inner scaling (argument)
283: - beta - outer scaling (result)
285: Notes:
286: Given a function f(x) specified by the FN type, the scaling parameters can
287: be used to realize the function beta*f(alpha*x). So when these values are given,
288: the procedure for function evaluation will first multiply the argument by alpha,
289: then evaluate the function itself, and finally scale the result by beta.
290: Likewise, these values are also considered when evaluating the derivative.
292: If you want to provide only one of the two scaling factors, set the other
293: one to 1.0.
295: Level: intermediate
297: .seealso: FNGetScale(), FNEvaluateFunction()
298: @*/
299: PetscErrorCode FNSetScale(FN fn,PetscScalar alpha,PetscScalar beta)300: {
305: if (PetscAbsScalar(alpha)==0.0 || PetscAbsScalar(beta)==0.0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_WRONG,"Scaling factors must be nonzero");
306: fn->alpha = alpha;
307: fn->beta = beta;
308: return(0);
309: }
311: /*@
312: FNGetScale - Gets the scaling parameters that define the matematical function.
314: Not Collective
316: Input Parameter:
317: . fn - the math function context
319: Output Parameters:
320: + alpha - inner scaling (argument)
321: - beta - outer scaling (result)
323: Level: intermediate
325: .seealso: FNSetScale()
326: @*/
327: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)328: {
331: if (alpha) *alpha = fn->alpha;
332: if (beta) *beta = fn->beta;
333: return(0);
334: }
336: /*@
337: FNSetMethod - Selects the method to be used to evaluate functions of matrices.
339: Logically Collective on FN341: Input Parameter:
342: + fn - the math function context
343: - meth - an index indentifying the method
345: Options Database Key:
346: . -fn_method <meth> - Sets the method
348: Notes:
349: In some FN types there are more than one algorithms available for computing
350: matrix functions. In that case, this function allows choosing the wanted method.
352: If meth is currently set to 0 (the default) and the input argument A of
353: FNEvaluateFunctionMat() is a symmetric/Hermitian matrix, then the computation
354: is done via the eigendecomposition of A, rather than with the general algorithm.
356: Level: intermediate
358: .seealso: FNGetMethod(), FNEvaluateFunctionMat()
359: @*/
360: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)361: {
365: if (meth<0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
366: if (meth>FN_MAX_SOLVE) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
367: fn->method = meth;
368: return(0);
369: }
371: /*@
372: FNGetMethod - Gets the method currently used in the FN.
374: Not Collective
376: Input Parameter:
377: . fn - the math function context
379: Output Parameter:
380: . meth - identifier of the method
382: Level: intermediate
384: .seealso: FNSetMethod()
385: @*/
386: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)387: {
391: *meth = fn->method;
392: return(0);
393: }
395: /*@
396: FNSetParallel - Selects the mode of operation in parallel runs.
398: Logically Collective on FN400: Input Parameter:
401: + fn - the math function context
402: - pmode - the parallel mode
404: Options Database Key:
405: . -fn_parallel <mode> - Sets the parallel mode, either 'redundant' or 'synchronized'
407: Notes:
408: This is relevant only when the function is evaluated on a matrix, with
409: either FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
411: In the 'redundant' parallel mode, all processes will make the computation
412: redundantly, starting from the same data, and producing the same result.
413: This result may be slightly different in the different processes if using a
414: multithreaded BLAS library, which may cause issues in ill-conditioned problems.
416: In the 'synchronized' parallel mode, only the first MPI process performs the
417: computation and then the computed matrix is broadcast to the other
418: processes in the communicator. This communication is done automatically at
419: the end of FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
421: Level: advanced
423: .seealso: FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec(), FNGetParallel()
424: @*/
425: PetscErrorCode FNSetParallel(FN fn,FNParallelType pmode)426: {
430: fn->pmode = pmode;
431: return(0);
432: }
434: /*@
435: FNGetParallel - Gets the mode of operation in parallel runs.
437: Not Collective
439: Input Parameter:
440: . fn - the math function context
442: Output Parameter:
443: . pmode - the parallel mode
445: Level: advanced
447: .seealso: FNSetParallel()
448: @*/
449: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)450: {
454: *pmode = fn->pmode;
455: return(0);
456: }
458: /*@
459: FNEvaluateFunction - Computes the value of the function f(x) for a given x.
461: Not collective
463: Input Parameters:
464: + fn - the math function context
465: - x - the value where the function must be evaluated
467: Output Parameter:
468: . y - the result of f(x)
470: Note:
471: Scaling factors are taken into account, so the actual function evaluation
472: will return beta*f(alpha*x).
474: Level: intermediate
476: .seealso: FNEvaluateDerivative(), FNEvaluateFunctionMat(), FNSetScale()
477: @*/
478: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)479: {
481: PetscScalar xf,yf;
487: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
488: xf = fn->alpha*x;
489: (*fn->ops->evaluatefunction)(fn,xf,&yf);
490: *y = fn->beta*yf;
491: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
492: return(0);
493: }
495: /*@
496: FNEvaluateDerivative - Computes the value of the derivative f'(x) for a given x.
498: Not Collective
500: Input Parameters:
501: + fn - the math function context
502: - x - the value where the derivative must be evaluated
504: Output Parameter:
505: . y - the result of f'(x)
507: Note:
508: Scaling factors are taken into account, so the actual derivative evaluation will
509: return alpha*beta*f'(alpha*x).
511: Level: intermediate
513: .seealso: FNEvaluateFunction()
514: @*/
515: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)516: {
518: PetscScalar xf,yf;
524: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
525: xf = fn->alpha*x;
526: (*fn->ops->evaluatederivative)(fn,xf,&yf);
527: *y = fn->alpha*fn->beta*yf;
528: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
529: return(0);
530: }
532: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)533: {
534: #if defined(PETSC_MISSING_LAPACK_SYEV) || defined(SLEPC_MISSING_LAPACK_LACPY)
536: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SYEV/LACPY - Lapack routines are unavailable");
537: #else
539: PetscInt i,j;
540: PetscBLASInt n,k,ld,lwork,info;
541: PetscScalar *Q,*W,*work,a,x,y,one=1.0,zero=0.0;
542: PetscReal *eig,dummy;
543: #if defined(PETSC_USE_COMPLEX)
544: PetscReal *rwork,rdummy;
545: #endif
548: PetscBLASIntCast(m,&n);
549: ld = n;
550: k = firstonly? 1: n;
552: /* workspace query and memory allocation */
553: lwork = -1;
554: #if defined(PETSC_USE_COMPLEX)
555: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,As,&ld,&dummy,&a,&lwork,&rdummy,&info));
556: PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
557: PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork);
558: #else
559: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,As,&ld,&dummy,&a,&lwork,&info));
560: PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
561: PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work);
562: #endif
564: /* compute eigendecomposition */
565: PetscStackCallBLAS("LAPACKlacpy",LAPACKlacpy_("L",&n,&n,As,&ld,Q,&ld));
566: #if defined(PETSC_USE_COMPLEX)
567: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
568: #else
569: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
570: #endif
571: SlepcCheckLapackInfo("syev",info);
573: /* W = f(Lambda)*Q' */
574: for (i=0;i<n;i++) {
575: x = fn->alpha*eig[i];
576: (*fn->ops->evaluatefunction)(fn,x,&y); /* y = f(x) */
577: for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
578: }
579: /* Bs = Q*W */
580: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
581: #if defined(PETSC_USE_COMPLEX)
582: PetscFree5(eig,Q,W,work,rwork);
583: #else
584: PetscFree4(eig,Q,W,work);
585: #endif
586: PetscLogFlops(9.0*n*n*n+2.0*n*n*n);
587: return(0);
588: #endif
589: }
591: /*
592: FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
593: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
594: decomposition of A is A=Q*D*Q'
595: */
596: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)597: {
599: PetscInt m;
600: PetscScalar *As,*Bs;
603: MatDenseGetArray(A,&As);
604: MatDenseGetArray(B,&Bs);
605: MatGetSize(A,&m,NULL);
606: FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE);
607: MatDenseRestoreArray(A,&As);
608: MatDenseRestoreArray(B,&Bs);
609: return(0);
610: }
612: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)613: {
615: PetscBool set,flg,symm=PETSC_FALSE;
616: PetscInt m,n;
617: PetscMPIInt size,rank;
618: PetscScalar *pF;
619: Mat M,F;
622: /* destination matrix */
623: F = B?B:A;
625: /* check symmetry of A */
626: MatIsHermitianKnown(A,&set,&flg);
627: symm = set? flg: PETSC_FALSE;
629: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
630: MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
631: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
632: 633: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
634: if (symm && !fn->method) { /* prefer diagonalization */
635: PetscInfo(fn,"Computing matrix function via diagonalization\n");
636: FNEvaluateFunctionMat_Sym_Default(fn,A,F);
637: } else {
638: /* scale argument */
639: if (fn->alpha!=(PetscScalar)1.0) {
640: FN_AllocateWorkMat(fn,A,&M);
641: MatScale(M,fn->alpha);
642: } else M = A;
643: if (fn->ops->evaluatefunctionmat[fn->method]) {
644: (*fn->ops->evaluatefunctionmat[fn->method])(fn,M,F);
645: } else if (!fn->method) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
646: else SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
647: if (fn->alpha!=(PetscScalar)1.0) {
648: FN_FreeWorkMat(fn,&M);
649: }
650: /* scale result */
651: MatScale(F,fn->beta);
652: }
653: PetscFPTrapPop();
654: }
655: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { /* synchronize */
656: MatGetSize(A,&m,&n);
657: MatDenseGetArray(F,&pF);
658: MPI_Bcast(pF,n*n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
659: MatDenseRestoreArray(F,&pF);
660: }
661: return(0);
662: }
664: /*@
665: FNEvaluateFunctionMat - Computes the value of the function f(A) for a given
666: matrix A, where the result is also a matrix.
668: Logically Collective on FN670: Input Parameters:
671: + fn - the math function context
672: - A - matrix on which the function must be evaluated
674: Output Parameter:
675: . B - (optional) matrix resulting from evaluating f(A)
677: Notes:
678: Matrix A must be a square sequential dense Mat, with all entries equal on
679: all processes (otherwise each process will compute different results).
680: If matrix B is provided, it must also be a square sequential dense Mat, and
681: both matrices must have the same dimensions. If B is NULL (or B=A) then the
682: function will perform an in-place computation, overwriting A with f(A).
684: If A is known to be real symmetric or complex Hermitian then it is
685: recommended to set the appropriate flag with MatSetOption(), because
686: symmetry can sometimes be exploited by the algorithm.
688: Scaling factors are taken into account, so the actual function evaluation
689: will return beta*f(alpha*A).
691: Level: advanced
693: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
694: @*/
695: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)696: {
698: PetscBool match,inplace=PETSC_FALSE;
699: PetscInt m,n,n1;
706: if (B) {
709: } else inplace = PETSC_TRUE;
710: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&match);
711: if (!match) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Mat A must be of type seqdense");
712: MatGetSize(A,&m,&n);
713: if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %D rows, %D cols)",m,n);
714: if (!inplace) {
715: PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&match);
716: if (!match) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Mat B must be of type seqdense");
717: n1 = n;
718: MatGetSize(B,&m,&n);
719: if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat B is not square (has %D rows, %D cols)",m,n);
720: if (n1!=n) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrices A and B must have the same dimension");
721: }
723: /* evaluate matrix function */
724: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
725: FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE);
726: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
727: return(0);
728: }
730: /*
731: FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
732: and then copies the first column.
733: */
734: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)735: {
737: Mat F;
740: FN_AllocateWorkMat(fn,A,&F);
741: if (fn->ops->evaluatefunctionmat[fn->method]) {
742: (*fn->ops->evaluatefunctionmat[fn->method])(fn,A,F);
743: } else if (!fn->method) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Matrix functions not implemented in this FN type");
744: else SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN type");
745: MatGetColumnVector(F,v,0);
746: FN_FreeWorkMat(fn,&F);
747: return(0);
748: }
750: /*
751: FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
752: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
753: decomposition of A is A=Q*D*Q'. Only the first column is computed.
754: */
755: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)756: {
758: PetscInt m;
759: PetscScalar *As,*vs;
762: MatDenseGetArray(A,&As);
763: VecGetArray(v,&vs);
764: MatGetSize(A,&m,NULL);
765: FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE);
766: MatDenseRestoreArray(A,&As);
767: VecRestoreArray(v,&vs);
768: return(0);
769: }
771: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)772: {
774: PetscBool set,flg,symm=PETSC_FALSE;
775: PetscInt m,n;
776: Mat M;
777: PetscMPIInt size,rank;
778: PetscScalar *pv;
781: /* check symmetry of A */
782: MatIsHermitianKnown(A,&set,&flg);
783: symm = set? flg: PETSC_FALSE;
785: /* evaluate matrix function */
786: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
787: MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
788: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
789: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
790: if (symm && !fn->method) { /* prefer diagonalization */
791: PetscInfo(fn,"Computing matrix function via diagonalization\n");
792: FNEvaluateFunctionMatVec_Sym_Default(fn,A,v);
793: } else {
794: /* scale argument */
795: if (fn->alpha!=(PetscScalar)1.0) {
796: FN_AllocateWorkMat(fn,A,&M);
797: MatScale(M,fn->alpha);
798: } else M = A;
799: if (fn->ops->evaluatefunctionmatvec[fn->method]) {
800: (*fn->ops->evaluatefunctionmatvec[fn->method])(fn,M,v);
801: } else {
802: FNEvaluateFunctionMatVec_Default(fn,M,v);
803: }
804: if (fn->alpha!=(PetscScalar)1.0) {
805: FN_FreeWorkMat(fn,&M);
806: }
807: /* scale result */
808: VecScale(v,fn->beta);
809: }
810: PetscFPTrapPop();
811: }
813: /* synchronize */
814: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
815: MatGetSize(A,&m,&n);
816: VecGetArray(v,&pv);
817: MPI_Bcast(pv,n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
818: VecRestoreArray(v,&pv);
819: }
820: return(0);
821: }
823: /*@
824: FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
825: for a given matrix A.
827: Logically Collective on FN829: Input Parameters:
830: + fn - the math function context
831: - A - matrix on which the function must be evaluated
833: Output Parameter:
834: . v - vector to hold the first column of f(A)
836: Notes:
837: This operation is similar to FNEvaluateFunctionMat() but returns only
838: the first column of f(A), hence saving computations in most cases.
840: Level: advanced
842: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
843: @*/
844: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)845: {
847: PetscBool match;
848: PetscInt m,n;
857: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&match);
858: if (!match) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Mat A must be of type seqdense");
859: MatGetSize(A,&m,&n);
860: if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %D rows, %D cols)",m,n);
861: VecGetSize(v,&m);
862: if (m!=n) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrix A and vector v must have the same size");
863: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
864: FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE);
865: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
866: return(0);
867: }
869: /*@
870: FNSetFromOptions - Sets FN options from the options database.
872: Collective on FN874: Input Parameters:
875: . fn - the math function context
877: Notes:
878: To see all options, run your program with the -help option.
880: Level: beginner
881: @*/
882: PetscErrorCode FNSetFromOptions(FN fn)883: {
885: char type[256];
886: PetscScalar array[2];
887: PetscInt k,meth;
888: PetscBool flg;
889: FNParallelType pmode;
893: FNRegisterAll();
894: PetscObjectOptionsBegin((PetscObject)fn);
895: PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,256,&flg);
896: if (flg) {
897: FNSetType(fn,type);
898: } else if (!((PetscObject)fn)->type_name) {
899: FNSetType(fn,FNRATIONAL);
900: }
902: k = 2;
903: array[0] = 0.0; array[1] = 0.0;
904: PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg);
905: if (flg) {
906: if (k<2) array[1] = 1.0;
907: FNSetScale(fn,array[0],array[1]);
908: }
910: PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg);
911: if (flg) { FNSetMethod(fn,meth); }
913: PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg);
914: if (flg) { FNSetParallel(fn,pmode); }
916: if (fn->ops->setfromoptions) {
917: (*fn->ops->setfromoptions)(PetscOptionsObject,fn);
918: }
919: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)fn);
920: PetscOptionsEnd();
921: return(0);
922: }
924: /*@C
925: FNView - Prints the FN data structure.
927: Collective on FN929: Input Parameters:
930: + fn - the math function context
931: - viewer - optional visualization context
933: Note:
934: The available visualization contexts include
935: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
936: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
937: output where only the first processor opens
938: the file. All other processors send their
939: data to the first processor to print.
941: The user can open an alternative visualization context with
942: PetscViewerASCIIOpen() - output to a specified file.
944: Level: beginner
945: @*/
946: PetscErrorCode FNView(FN fn,PetscViewer viewer)947: {
948: PetscBool isascii;
950: PetscMPIInt size;
954: if (!viewer) {
955: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer);
956: }
959: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
960: if (isascii) {
961: PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer);
962: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
963: if (size>1) {
964: PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",FNParallelTypes[fn->pmode]);
965: }
966: if (fn->ops->view) {
967: PetscViewerASCIIPushTab(viewer);
968: (*fn->ops->view)(fn,viewer);
969: PetscViewerASCIIPopTab(viewer);
970: }
971: }
972: return(0);
973: }
975: /*@
976: FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
977: different communicator.
979: Collective on FN981: Input Parameters:
982: + fn - the math function context
983: - comm - MPI communicator
985: Output Parameter:
986: . newfn - location to put the new FN context
988: Note:
989: In order to use the same MPI communicator as in the original object,
990: use PetscObjectComm((PetscObject)fn).
992: Level: developer
994: .seealso: FNCreate()
995: @*/
996: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)997: {
999: FNType type;
1000: PetscScalar alpha,beta;
1001: PetscInt meth;
1002: FNParallelType ptype;
1008: FNCreate(comm,newfn);
1009: FNGetType(fn,&type);
1010: FNSetType(*newfn,type);
1011: FNGetScale(fn,&alpha,&beta);
1012: FNSetScale(*newfn,alpha,beta);
1013: FNGetMethod(fn,&meth);
1014: FNSetMethod(*newfn,meth);
1015: FNGetParallel(fn,&ptype);
1016: FNSetParallel(*newfn,ptype);
1017: if (fn->ops->duplicate) {
1018: (*fn->ops->duplicate)(fn,comm,newfn);
1019: }
1020: return(0);
1021: }
1023: /*@
1024: FNDestroy - Destroys FN context that was created with FNCreate().
1026: Collective on FN1028: Input Parameter:
1029: . fn - the math function context
1031: Level: beginner
1033: .seealso: FNCreate()
1034: @*/
1035: PetscErrorCode FNDestroy(FN *fn)1036: {
1038: PetscInt i;
1041: if (!*fn) return(0);
1043: if (--((PetscObject)(*fn))->refct > 0) { *fn = 0; return(0); }
1044: if ((*fn)->ops->destroy) { (*(*fn)->ops->destroy)(*fn); }
1045: for (i=0;i<(*fn)->nw;i++) {
1046: MatDestroy(&(*fn)->W[i]);
1047: }
1048: PetscHeaderDestroy(fn);
1049: return(0);
1050: }
1052: /*@C
1053: FNRegister - Adds a mathematical function to the FN package.
1055: Not collective
1057: Input Parameters:
1058: + name - name of a new user-defined FN1059: - function - routine to create context
1061: Notes:
1062: FNRegister() may be called multiple times to add several user-defined functions.
1064: Level: advanced
1066: .seealso: FNRegisterAll()
1067: @*/
1068: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))1069: {
1073: FNInitializePackage();
1074: PetscFunctionListAdd(&FNList,name,function);
1075: return(0);
1076: }