1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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: fn->alpha = alpha;
306: fn->beta = beta;
307: return(0);
308: }
310: /*@
311: FNGetScale - Gets the scaling parameters that define the matematical function.
313: Not Collective
315: Input Parameter:
316: . fn - the math function context
318: Output Parameters:
319: + alpha - inner scaling (argument)
320: - beta - outer scaling (result)
322: Level: intermediate
324: .seealso: FNSetScale()
325: @*/
326: PetscErrorCode FNGetScale(FN fn,PetscScalar *alpha,PetscScalar *beta)327: {
330: if (alpha) *alpha = fn->alpha;
331: if (beta) *beta = fn->beta;
332: return(0);
333: }
335: /*@
336: FNSetMethod - Selects the method to be used to evaluate functions of matrices.
338: Logically Collective on FN340: Input Parameter:
341: + fn - the math function context
342: - meth - an index indentifying the method
344: Options Database Key:
345: . -fn_method <meth> - Sets the method
347: Notes:
348: In some FN types there are more than one algorithms available for computing
349: matrix functions. In that case, this function allows choosing the wanted method.
351: If meth is currently set to 0 (the default) and the input argument A of
352: FNEvaluateFunctionMat() is a symmetric/Hermitian matrix, then the computation
353: is done via the eigendecomposition of A, rather than with the general algorithm.
355: Level: intermediate
357: .seealso: FNGetMethod(), FNEvaluateFunctionMat()
358: @*/
359: PetscErrorCode FNSetMethod(FN fn,PetscInt meth)360: {
364: if (meth<0) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
365: if (meth>FN_MAX_SOLVE) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
366: fn->method = meth;
367: return(0);
368: }
370: /*@
371: FNGetMethod - Gets the method currently used in the FN.
373: Not Collective
375: Input Parameter:
376: . fn - the math function context
378: Output Parameter:
379: . meth - identifier of the method
381: Level: intermediate
383: .seealso: FNSetMethod()
384: @*/
385: PetscErrorCode FNGetMethod(FN fn,PetscInt *meth)386: {
390: *meth = fn->method;
391: return(0);
392: }
394: /*@
395: FNSetParallel - Selects the mode of operation in parallel runs.
397: Logically Collective on FN399: Input Parameter:
400: + fn - the math function context
401: - pmode - the parallel mode
403: Options Database Key:
404: . -fn_parallel <mode> - Sets the parallel mode, either 'redundant' or 'synchronized'
406: Notes:
407: This is relevant only when the function is evaluated on a matrix, with
408: either FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
410: In the 'redundant' parallel mode, all processes will make the computation
411: redundantly, starting from the same data, and producing the same result.
412: This result may be slightly different in the different processes if using a
413: multithreaded BLAS library, which may cause issues in ill-conditioned problems.
415: In the 'synchronized' parallel mode, only the first MPI process performs the
416: computation and then the computed matrix is broadcast to the other
417: processes in the communicator. This communication is done automatically at
418: the end of FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec().
420: Level: advanced
422: .seealso: FNEvaluateFunctionMat() or FNEvaluateFunctionMatVec(), FNGetParallel()
423: @*/
424: PetscErrorCode FNSetParallel(FN fn,FNParallelType pmode)425: {
429: fn->pmode = pmode;
430: return(0);
431: }
433: /*@
434: FNGetParallel - Gets the mode of operation in parallel runs.
436: Not Collective
438: Input Parameter:
439: . fn - the math function context
441: Output Parameter:
442: . pmode - the parallel mode
444: Level: advanced
446: .seealso: FNSetParallel()
447: @*/
448: PetscErrorCode FNGetParallel(FN fn,FNParallelType *pmode)449: {
453: *pmode = fn->pmode;
454: return(0);
455: }
457: /*@
458: FNEvaluateFunction - Computes the value of the function f(x) for a given x.
460: Not collective
462: Input Parameters:
463: + fn - the math function context
464: - x - the value where the function must be evaluated
466: Output Parameter:
467: . y - the result of f(x)
469: Note:
470: Scaling factors are taken into account, so the actual function evaluation
471: will return beta*f(alpha*x).
473: Level: intermediate
475: .seealso: FNEvaluateDerivative(), FNEvaluateFunctionMat(), FNSetScale()
476: @*/
477: PetscErrorCode FNEvaluateFunction(FN fn,PetscScalar x,PetscScalar *y)478: {
480: PetscScalar xf,yf;
486: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
487: xf = fn->alpha*x;
488: (*fn->ops->evaluatefunction)(fn,xf,&yf);
489: *y = fn->beta*yf;
490: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
491: return(0);
492: }
494: /*@
495: FNEvaluateDerivative - Computes the value of the derivative f'(x) for a given x.
497: Not Collective
499: Input Parameters:
500: + fn - the math function context
501: - x - the value where the derivative must be evaluated
503: Output Parameter:
504: . y - the result of f'(x)
506: Note:
507: Scaling factors are taken into account, so the actual derivative evaluation will
508: return alpha*beta*f'(alpha*x).
510: Level: intermediate
512: .seealso: FNEvaluateFunction()
513: @*/
514: PetscErrorCode FNEvaluateDerivative(FN fn,PetscScalar x,PetscScalar *y)515: {
517: PetscScalar xf,yf;
523: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
524: xf = fn->alpha*x;
525: (*fn->ops->evaluatederivative)(fn,xf,&yf);
526: *y = fn->alpha*fn->beta*yf;
527: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
528: return(0);
529: }
531: static PetscErrorCode FNEvaluateFunctionMat_Sym_Private(FN fn,PetscScalar *As,PetscScalar *Bs,PetscInt m,PetscBool firstonly)532: {
533: #if defined(PETSC_MISSING_LAPACK_SYEV) || defined(SLEPC_MISSING_LAPACK_LACPY)
535: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SYEV/LACPY - Lapack routines are unavailable");
536: #else
538: PetscInt i,j;
539: PetscBLASInt n,k,ld,lwork,info;
540: PetscScalar *Q,*W,*work,a,x,y,one=1.0,zero=0.0;
541: PetscReal *eig,dummy;
542: #if defined(PETSC_USE_COMPLEX)
543: PetscReal *rwork,rdummy;
544: #endif
547: PetscBLASIntCast(m,&n);
548: ld = n;
549: k = firstonly? 1: n;
551: /* workspace query and memory allocation */
552: lwork = -1;
553: #if defined(PETSC_USE_COMPLEX)
554: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,As,&ld,&dummy,&a,&lwork,&rdummy,&info));
555: PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
556: PetscMalloc5(m,&eig,m*m,&Q,m*k,&W,lwork,&work,PetscMax(1,3*m-2),&rwork);
557: #else
558: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,As,&ld,&dummy,&a,&lwork,&info));
559: PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
560: PetscMalloc4(m,&eig,m*m,&Q,m*k,&W,lwork,&work);
561: #endif
563: /* compute eigendecomposition */
564: PetscStackCallBLAS("LAPACKlacpy",LAPACKlacpy_("L",&n,&n,As,&ld,Q,&ld));
565: #if defined(PETSC_USE_COMPLEX)
566: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
567: #else
568: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
569: #endif
570: SlepcCheckLapackInfo("syev",info);
572: /* W = f(Lambda)*Q' */
573: for (i=0;i<n;i++) {
574: x = fn->alpha*eig[i];
575: (*fn->ops->evaluatefunction)(fn,x,&y); /* y = f(x) */
576: for (j=0;j<k;j++) W[i+j*ld] = PetscConj(Q[j+i*ld])*fn->beta*y;
577: }
578: /* Bs = Q*W */
579: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,Bs,&ld));
580: #if defined(PETSC_USE_COMPLEX)
581: PetscFree5(eig,Q,W,work,rwork);
582: #else
583: PetscFree4(eig,Q,W,work);
584: #endif
585: PetscLogFlops(9.0*n*n*n+2.0*n*n*n);
586: return(0);
587: #endif
588: }
590: /*
591: FNEvaluateFunctionMat_Sym_Default - given a symmetric matrix A,
592: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
593: decomposition of A is A=Q*D*Q'
594: */
595: static PetscErrorCode FNEvaluateFunctionMat_Sym_Default(FN fn,Mat A,Mat B)596: {
598: PetscInt m;
599: PetscScalar *As,*Bs;
602: MatDenseGetArray(A,&As);
603: MatDenseGetArray(B,&Bs);
604: MatGetSize(A,&m,NULL);
605: FNEvaluateFunctionMat_Sym_Private(fn,As,Bs,m,PETSC_FALSE);
606: MatDenseRestoreArray(A,&As);
607: MatDenseRestoreArray(B,&Bs);
608: return(0);
609: }
611: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)612: {
614: PetscBool set,flg,symm=PETSC_FALSE;
615: PetscInt m,n;
616: PetscMPIInt size,rank;
617: PetscScalar *pF;
618: Mat M,F;
621: /* destination matrix */
622: F = B?B:A;
624: /* check symmetry of A */
625: MatIsHermitianKnown(A,&set,&flg);
626: symm = set? flg: PETSC_FALSE;
628: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
629: MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
630: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
631: 632: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
633: if (symm && !fn->method) { /* prefer diagonalization */
634: PetscInfo(fn,"Computing matrix function via diagonalization\n");
635: FNEvaluateFunctionMat_Sym_Default(fn,A,F);
636: } else {
637: /* scale argument */
638: if (fn->alpha!=(PetscScalar)1.0) {
639: FN_AllocateWorkMat(fn,A,&M);
640: MatScale(M,fn->alpha);
641: } else M = A;
642: if (fn->ops->evaluatefunctionmat[fn->method]) {
643: (*fn->ops->evaluatefunctionmat[fn->method])(fn,M,F);
644: } else SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN");
645: if (fn->alpha!=(PetscScalar)1.0) {
646: FN_FreeWorkMat(fn,&M);
647: }
648: /* scale result */
649: MatScale(F,fn->beta);
650: }
651: PetscFPTrapPop();
652: }
653: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { /* synchronize */
654: MatGetSize(A,&m,&n);
655: MatDenseGetArray(F,&pF);
656: MPI_Bcast(pF,n*n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
657: MatDenseRestoreArray(F,&pF);
658: }
659: return(0);
660: }
662: /*@
663: FNEvaluateFunctionMat - Computes the value of the function f(A) for a given
664: matrix A, where the result is also a matrix.
666: Logically Collective on FN668: Input Parameters:
669: + fn - the math function context
670: - A - matrix on which the function must be evaluated
672: Output Parameter:
673: . B - (optional) matrix resulting from evaluating f(A)
675: Notes:
676: Matrix A must be a square sequential dense Mat, with all entries equal on
677: all processes (otherwise each process will compute different results).
678: If matrix B is provided, it must also be a square sequential dense Mat, and
679: both matrices must have the same dimensions. If B is NULL (or B=A) then the
680: function will perform an in-place computation, overwriting A with f(A).
682: If A is known to be real symmetric or complex Hermitian then it is
683: recommended to set the appropriate flag with MatSetOption(), because
684: symmetry can sometimes be exploited by the algorithm.
686: Scaling factors are taken into account, so the actual function evaluation
687: will return beta*f(alpha*A).
689: Level: advanced
691: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
692: @*/
693: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)694: {
696: PetscBool match,inplace=PETSC_FALSE;
697: PetscInt m,n,n1;
704: if (B) {
707: } else inplace = PETSC_TRUE;
708: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&match);
709: if (!match) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Mat A must be of type seqdense");
710: MatGetSize(A,&m,&n);
711: if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %D rows, %D cols)",m,n);
712: if (!inplace) {
713: PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&match);
714: if (!match) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Mat B must be of type seqdense");
715: n1 = n;
716: MatGetSize(B,&m,&n);
717: if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat B is not square (has %D rows, %D cols)",m,n);
718: if (n1!=n) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrices A and B must have the same dimension");
719: }
721: /* evaluate matrix function */
722: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
723: FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE);
724: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
725: return(0);
726: }
728: /*
729: FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
730: and then copies the first column.
731: */
732: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)733: {
735: Mat F;
738: FN_AllocateWorkMat(fn,A,&F);
739: if (fn->ops->evaluatefunctionmat[fn->method]) {
740: (*fn->ops->evaluatefunctionmat[fn->method])(fn,A,F);
741: } else SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this FN");
742: MatGetColumnVector(F,v,0);
743: FN_FreeWorkMat(fn,&F);
744: return(0);
745: }
747: /*
748: FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
749: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
750: decomposition of A is A=Q*D*Q'. Only the first column is computed.
751: */
752: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)753: {
755: PetscInt m;
756: PetscScalar *As,*vs;
759: MatDenseGetArray(A,&As);
760: VecGetArray(v,&vs);
761: MatGetSize(A,&m,NULL);
762: FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE);
763: MatDenseRestoreArray(A,&As);
764: VecRestoreArray(v,&vs);
765: return(0);
766: }
768: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)769: {
771: PetscBool set,flg,symm=PETSC_FALSE;
772: PetscInt m,n;
773: Mat M;
774: PetscMPIInt size,rank;
775: PetscScalar *pv;
778: /* check symmetry of A */
779: MatIsHermitianKnown(A,&set,&flg);
780: symm = set? flg: PETSC_FALSE;
782: /* evaluate matrix function */
783: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
784: MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
785: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
786: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
787: if (symm && !fn->method) { /* prefer diagonalization */
788: PetscInfo(fn,"Computing matrix function via diagonalization\n");
789: FNEvaluateFunctionMatVec_Sym_Default(fn,A,v);
790: } else {
791: /* scale argument */
792: if (fn->alpha!=(PetscScalar)1.0) {
793: FN_AllocateWorkMat(fn,A,&M);
794: MatScale(M,fn->alpha);
795: } else M = A;
796: if (fn->ops->evaluatefunctionmatvec[fn->method]) {
797: (*fn->ops->evaluatefunctionmatvec[fn->method])(fn,M,v);
798: } else {
799: FNEvaluateFunctionMatVec_Default(fn,M,v);
800: }
801: if (fn->alpha!=(PetscScalar)1.0) {
802: FN_FreeWorkMat(fn,&M);
803: }
804: /* scale result */
805: VecScale(v,fn->beta);
806: }
807: PetscFPTrapPop();
808: }
810: /* synchronize */
811: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
812: MatGetSize(A,&m,&n);
813: VecGetArray(v,&pv);
814: MPI_Bcast(pv,n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
815: VecRestoreArray(v,&pv);
816: }
817: return(0);
818: }
820: /*@
821: FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
822: for a given matrix A.
824: Logically Collective on FN826: Input Parameters:
827: + fn - the math function context
828: - A - matrix on which the function must be evaluated
830: Output Parameter:
831: . v - vector to hold the first column of f(A)
833: Notes:
834: This operation is similar to FNEvaluateFunctionMat() but returns only
835: the first column of f(A), hence saving computations in most cases.
837: Level: advanced
839: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
840: @*/
841: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)842: {
844: PetscBool match;
845: PetscInt m,n;
854: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&match);
855: if (!match) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_SUP,"Mat A must be of type seqdense");
856: MatGetSize(A,&m,&n);
857: if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Mat A is not square (has %D rows, %D cols)",m,n);
858: VecGetSize(v,&m);
859: if (m!=n) SETERRQ(PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_SIZ,"Matrix A and vector v must have the same size");
860: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
861: FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE);
862: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
863: return(0);
864: }
866: /*@
867: FNSetFromOptions - Sets FN options from the options database.
869: Collective on FN871: Input Parameters:
872: . fn - the math function context
874: Notes:
875: To see all options, run your program with the -help option.
877: Level: beginner
878: @*/
879: PetscErrorCode FNSetFromOptions(FN fn)880: {
882: char type[256];
883: PetscScalar array[2];
884: PetscInt k,meth;
885: PetscBool flg;
886: FNParallelType pmode;
890: FNRegisterAll();
891: PetscObjectOptionsBegin((PetscObject)fn);
892: PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,256,&flg);
893: if (flg) {
894: FNSetType(fn,type);
895: } else if (!((PetscObject)fn)->type_name) {
896: FNSetType(fn,FNRATIONAL);
897: }
899: k = 2;
900: array[0] = 0.0; array[1] = 0.0;
901: PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg);
902: if (flg) {
903: if (k<2) array[1] = 1.0;
904: FNSetScale(fn,array[0],array[1]);
905: }
907: PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg);
908: if (flg) { FNSetMethod(fn,meth); }
910: PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg);
911: if (flg) { FNSetParallel(fn,pmode); }
913: if (fn->ops->setfromoptions) {
914: (*fn->ops->setfromoptions)(PetscOptionsObject,fn);
915: }
916: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)fn);
917: PetscOptionsEnd();
918: return(0);
919: }
921: /*@C
922: FNView - Prints the FN data structure.
924: Collective on FN926: Input Parameters:
927: + fn - the math function context
928: - viewer - optional visualization context
930: Note:
931: The available visualization contexts include
932: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
933: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
934: output where only the first processor opens
935: the file. All other processors send their
936: data to the first processor to print.
938: The user can open an alternative visualization context with
939: PetscViewerASCIIOpen() - output to a specified file.
941: Level: beginner
942: @*/
943: PetscErrorCode FNView(FN fn,PetscViewer viewer)944: {
945: PetscBool isascii;
946: PetscInt tabs;
948: PetscMPIInt size;
952: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)fn));
955: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
956: if (isascii) {
957: PetscViewerASCIIGetTab(viewer,&tabs);
958: PetscViewerASCIISetTab(viewer,((PetscObject)fn)->tablevel);
959: PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer);
960: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
961: if (size>1) {
962: PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",FNParallelTypes[fn->pmode]);
963: }
964: if (fn->ops->view) {
965: PetscViewerASCIIPushTab(viewer);
966: (*fn->ops->view)(fn,viewer);
967: PetscViewerASCIIPopTab(viewer);
968: }
969: PetscViewerASCIISetTab(viewer,tabs);
970: }
971: return(0);
972: }
974: /*@
975: FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
976: different communicator.
978: Collective on FN980: Input Parameters:
981: + fn - the math function context
982: - comm - MPI communicator
984: Output Parameter:
985: . newfn - location to put the new FN context
987: Note:
988: In order to use the same MPI communicator as in the original object,
989: use PetscObjectComm((PetscObject)fn).
991: Level: developer
993: .seealso: FNCreate()
994: @*/
995: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)996: {
998: FNType type;
999: PetscScalar alpha,beta;
1000: PetscInt meth;
1001: FNParallelType ptype;
1007: FNCreate(comm,newfn);
1008: FNGetType(fn,&type);
1009: FNSetType(*newfn,type);
1010: FNGetScale(fn,&alpha,&beta);
1011: FNSetScale(*newfn,alpha,beta);
1012: FNGetMethod(fn,&meth);
1013: FNSetMethod(*newfn,meth);
1014: FNGetParallel(fn,&ptype);
1015: FNSetParallel(*newfn,ptype);
1016: if (fn->ops->duplicate) {
1017: (*fn->ops->duplicate)(fn,comm,newfn);
1018: }
1019: return(0);
1020: }
1022: /*@
1023: FNDestroy - Destroys FN context that was created with FNCreate().
1025: Collective on FN1027: Input Parameter:
1028: . fn - the math function context
1030: Level: beginner
1032: .seealso: FNCreate()
1033: @*/
1034: PetscErrorCode FNDestroy(FN *fn)1035: {
1037: PetscInt i;
1040: if (!*fn) return(0);
1042: if (--((PetscObject)(*fn))->refct > 0) { *fn = 0; return(0); }
1043: if ((*fn)->ops->destroy) { (*(*fn)->ops->destroy)(*fn); }
1044: for (i=0;i<(*fn)->nw;i++) {
1045: MatDestroy(&(*fn)->W[i]);
1046: }
1047: PetscHeaderDestroy(fn);
1048: return(0);
1049: }
1051: /*@C
1052: FNRegister - Adds a mathematical function to the FN package.
1054: Not collective
1056: Input Parameters:
1057: + name - name of a new user-defined FN1058: - function - routine to create context
1060: Notes:
1061: FNRegister() may be called multiple times to add several user-defined functions.
1063: Level: advanced
1065: .seealso: FNRegisterAll()
1066: @*/
1067: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))1068: {
1072: PetscFunctionListAdd(&FNList,name,function);
1073: return(0);
1074: }