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 = 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: {
35: PetscFunctionListDestroy(&FNList);
36: FNPackageInitialized = PETSC_FALSE;
37: FNRegisterAllCalled = PETSC_FALSE;
38: PetscFunctionReturn(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) PetscFunctionReturn(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: PetscFunctionReturn(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 = 0;
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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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) PetscFunctionReturn(0);
221: PetscFunctionListFind(FNList,type,&r);
224: if (fn->ops->destroy) (*fn->ops->destroy)(fn);
225: PetscMemzero(fn->ops,sizeof(struct _FNOps));
227: PetscObjectChangeTypeName((PetscObject)fn,type);
228: (*r)(fn);
229: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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: (*fn->ops->evaluatefunction)(fn,xf,&yf);
462: *y = fn->beta*yf;
463: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
464: PetscFunctionReturn(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: (*fn->ops->evaluatederivative)(fn,xf,&yf);
497: *y = fn->alpha*fn->beta*yf;
498: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
499: PetscFunctionReturn(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: PetscStackCallBLAS("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: PetscStackCallBLAS("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: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
532: #else
533: PetscStackCallBLAS("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: (*fn->ops->evaluatefunction)(fn,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: PetscStackCallBLAS("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: PetscFunctionReturn(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: PetscFunctionReturn(0);
572: }
574: PetscErrorCode FNEvaluateFunctionMat_Private(FN fn,Mat A,Mat B,PetscBool sync)575: {
576: PetscBool set,flg,symm=PETSC_FALSE;
577: PetscInt m,n;
578: PetscMPIInt size,rank;
579: PetscScalar *pF;
580: Mat M,F;
582: /* destination matrix */
583: F = B?B:A;
585: /* check symmetry of A */
586: MatIsHermitianKnown(A,&set,&flg);
587: symm = set? flg: PETSC_FALSE;
589: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
590: MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
591: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
593: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
594: if (symm && !fn->method) { /* prefer diagonalization */
595: PetscInfo(fn,"Computing matrix function via diagonalization\n");
596: FNEvaluateFunctionMat_Sym_Default(fn,A,F);
597: } else {
598: /* scale argument */
599: if (fn->alpha!=(PetscScalar)1.0) {
600: FN_AllocateWorkMat(fn,A,&M);
601: MatScale(M,fn->alpha);
602: } else M = A;
603: if (fn->ops->evaluatefunctionmat[fn->method]) (*fn->ops->evaluatefunctionmat[fn->method])(fn,M,F);
604: else {
607: }
608: if (fn->alpha!=(PetscScalar)1.0) FN_FreeWorkMat(fn,&M);
609: /* scale result */
610: MatScale(F,fn->beta);
611: }
612: PetscFPTrapPop();
613: }
614: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) { /* synchronize */
615: MatGetSize(A,&m,&n);
616: MatDenseGetArray(F,&pF);
617: MPI_Bcast(pF,n*n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
618: MatDenseRestoreArray(F,&pF);
619: }
620: PetscFunctionReturn(0);
621: }
623: /*@
624: FNEvaluateFunctionMat - Computes the value of the function f(A) for a given
625: matrix A, where the result is also a matrix.
627: Logically Collective on fn
629: Input Parameters:
630: + fn - the math function context
631: - A - matrix on which the function must be evaluated
633: Output Parameter:
634: . B - (optional) matrix resulting from evaluating f(A)
636: Notes:
637: Matrix A must be a square sequential dense Mat, with all entries equal on
638: all processes (otherwise each process will compute different results).
639: If matrix B is provided, it must also be a square sequential dense Mat, and
640: both matrices must have the same dimensions. If B is NULL (or B=A) then the
641: function will perform an in-place computation, overwriting A with f(A).
643: If A is known to be real symmetric or complex Hermitian then it is
644: recommended to set the appropriate flag with MatSetOption(), because
645: symmetry can sometimes be exploited by the algorithm.
647: Scaling factors are taken into account, so the actual function evaluation
648: will return beta*f(alpha*A).
650: Level: advanced
652: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMatVec(), FNSetMethod()
653: @*/
654: PetscErrorCode FNEvaluateFunctionMat(FN fn,Mat A,Mat B)655: {
656: PetscBool inplace=PETSC_FALSE;
657: PetscInt m,n,n1;
663: if (B) {
666: } else inplace = PETSC_TRUE;
668: MatGetSize(A,&m,&n);
670: if (!inplace) {
672: n1 = n;
673: MatGetSize(B,&m,&n);
676: }
678: /* evaluate matrix function */
679: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
680: FNEvaluateFunctionMat_Private(fn,A,B,PETSC_TRUE);
681: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
682: PetscFunctionReturn(0);
683: }
685: /*
686: FNEvaluateFunctionMatVec_Default - computes the full matrix f(A)
687: and then copies the first column.
688: */
689: static PetscErrorCode FNEvaluateFunctionMatVec_Default(FN fn,Mat A,Vec v)690: {
691: Mat F;
693: FN_AllocateWorkMat(fn,A,&F);
694: if (fn->ops->evaluatefunctionmat[fn->method]) (*fn->ops->evaluatefunctionmat[fn->method])(fn,A,F);
695: else {
698: }
699: MatGetColumnVector(F,v,0);
700: FN_FreeWorkMat(fn,&F);
701: PetscFunctionReturn(0);
702: }
704: /*
705: FNEvaluateFunctionMatVec_Sym_Default - given a symmetric matrix A,
706: compute the matrix function as f(A)=Q*f(D)*Q' where the spectral
707: decomposition of A is A=Q*D*Q'. Only the first column is computed.
708: */
709: static PetscErrorCode FNEvaluateFunctionMatVec_Sym_Default(FN fn,Mat A,Vec v)710: {
711: PetscInt m;
712: const PetscScalar *As;
713: PetscScalar *vs;
715: MatDenseGetArrayRead(A,&As);
716: VecGetArray(v,&vs);
717: MatGetSize(A,&m,NULL);
718: FNEvaluateFunctionMat_Sym_Private(fn,As,vs,m,PETSC_TRUE);
719: MatDenseRestoreArrayRead(A,&As);
720: VecRestoreArray(v,&vs);
721: PetscFunctionReturn(0);
722: }
724: PetscErrorCode FNEvaluateFunctionMatVec_Private(FN fn,Mat A,Vec v,PetscBool sync)725: {
726: PetscBool set,flg,symm=PETSC_FALSE;
727: PetscInt m,n;
728: Mat M;
729: PetscMPIInt size,rank;
730: PetscScalar *pv;
732: /* check symmetry of A */
733: MatIsHermitianKnown(A,&set,&flg);
734: symm = set? flg: PETSC_FALSE;
736: /* evaluate matrix function */
737: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
738: MPI_Comm_rank(PetscObjectComm((PetscObject)fn),&rank);
739: if (size==1 || fn->pmode==FN_PARALLEL_REDUNDANT || (fn->pmode==FN_PARALLEL_SYNCHRONIZED && !rank)) {
740: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
741: if (symm && !fn->method) { /* prefer diagonalization */
742: PetscInfo(fn,"Computing matrix function via diagonalization\n");
743: FNEvaluateFunctionMatVec_Sym_Default(fn,A,v);
744: } else {
745: /* scale argument */
746: if (fn->alpha!=(PetscScalar)1.0) {
747: FN_AllocateWorkMat(fn,A,&M);
748: MatScale(M,fn->alpha);
749: } else M = A;
750: if (fn->ops->evaluatefunctionmatvec[fn->method]) (*fn->ops->evaluatefunctionmatvec[fn->method])(fn,M,v);
751: else FNEvaluateFunctionMatVec_Default(fn,M,v);
752: if (fn->alpha!=(PetscScalar)1.0) FN_FreeWorkMat(fn,&M);
753: /* scale result */
754: VecScale(v,fn->beta);
755: }
756: PetscFPTrapPop();
757: }
759: /* synchronize */
760: if (size>1 && fn->pmode==FN_PARALLEL_SYNCHRONIZED && sync) {
761: MatGetSize(A,&m,&n);
762: VecGetArray(v,&pv);
763: MPI_Bcast(pv,n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)fn));
764: VecRestoreArray(v,&pv);
765: }
766: PetscFunctionReturn(0);
767: }
769: /*@
770: FNEvaluateFunctionMatVec - Computes the first column of the matrix f(A)
771: for a given matrix A.
773: Logically Collective on fn
775: Input Parameters:
776: + fn - the math function context
777: - A - matrix on which the function must be evaluated
779: Output Parameter:
780: . v - vector to hold the first column of f(A)
782: Notes:
783: This operation is similar to FNEvaluateFunctionMat() but returns only
784: the first column of f(A), hence saving computations in most cases.
786: Level: advanced
788: .seealso: FNEvaluateFunction(), FNEvaluateFunctionMat(), FNSetMethod()
789: @*/
790: PetscErrorCode FNEvaluateFunctionMatVec(FN fn,Mat A,Vec v)791: {
792: PetscInt m,n;
801: MatGetSize(A,&m,&n);
803: VecGetSize(v,&m);
805: PetscLogEventBegin(FN_Evaluate,fn,0,0,0);
806: FNEvaluateFunctionMatVec_Private(fn,A,v,PETSC_TRUE);
807: PetscLogEventEnd(FN_Evaluate,fn,0,0,0);
808: PetscFunctionReturn(0);
809: }
811: /*@
812: FNSetFromOptions - Sets FN options from the options database.
814: Collective on fn
816: Input Parameters:
817: . fn - the math function context
819: Notes:
820: To see all options, run your program with the -help option.
822: Level: beginner
824: .seealso: FNSetOptionsPrefix()
825: @*/
826: PetscErrorCode FNSetFromOptions(FN fn)827: {
829: char type[256];
830: PetscScalar array[2];
831: PetscInt k,meth;
832: PetscBool flg;
833: FNParallelType pmode;
836: FNRegisterAll();
837: ierr = PetscObjectOptionsBegin((PetscObject)fn);
838: PetscOptionsFList("-fn_type","Math function type","FNSetType",FNList,(char*)(((PetscObject)fn)->type_name?((PetscObject)fn)->type_name:FNRATIONAL),type,sizeof(type),&flg);
839: if (flg) FNSetType(fn,type);
840: else if (!((PetscObject)fn)->type_name) FNSetType(fn,FNRATIONAL);
842: k = 2;
843: array[0] = 0.0; array[1] = 0.0;
844: PetscOptionsScalarArray("-fn_scale","Scale factors (one or two scalar values separated with a comma without spaces)","FNSetScale",array,&k,&flg);
845: if (flg) {
846: if (k<2) array[1] = 1.0;
847: FNSetScale(fn,array[0],array[1]);
848: }
850: PetscOptionsInt("-fn_method","Method to be used for computing matrix functions","FNSetMethod",fn->method,&meth,&flg);
851: if (flg) FNSetMethod(fn,meth);
853: PetscOptionsEnum("-fn_parallel","Operation mode in parallel runs","FNSetParallel",FNParallelTypes,(PetscEnum)fn->pmode,(PetscEnum*)&pmode,&flg);
854: if (flg) FNSetParallel(fn,pmode);
856: if (fn->ops->setfromoptions) (*fn->ops->setfromoptions)(PetscOptionsObject,fn);
857: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)fn);
858: ierr = PetscOptionsEnd();
859: PetscFunctionReturn(0);
860: }
862: /*@C
863: FNView - Prints the FN data structure.
865: Collective on fn
867: Input Parameters:
868: + fn - the math function context
869: - viewer - optional visualization context
871: Note:
872: The available visualization contexts include
873: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
874: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
875: output where only the first processor opens
876: the file. All other processors send their
877: data to the first processor to print.
879: The user can open an alternative visualization context with
880: PetscViewerASCIIOpen() - output to a specified file.
882: Level: beginner
884: .seealso: FNCreate()
885: @*/
886: PetscErrorCode FNView(FN fn,PetscViewer viewer)887: {
888: PetscBool isascii;
889: PetscMPIInt size;
892: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fn),&viewer);
895: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
896: if (isascii) {
897: PetscObjectPrintClassNamePrefixType((PetscObject)fn,viewer);
898: MPI_Comm_size(PetscObjectComm((PetscObject)fn),&size);
899: if (size>1) PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",FNParallelTypes[fn->pmode]);
900: if (fn->ops->view) {
901: PetscViewerASCIIPushTab(viewer);
902: (*fn->ops->view)(fn,viewer);
903: PetscViewerASCIIPopTab(viewer);
904: }
905: }
906: PetscFunctionReturn(0);
907: }
909: /*@C
910: FNViewFromOptions - View from options
912: Collective on FN914: Input Parameters:
915: + fn - the math function context
916: . obj - optional object
917: - name - command line option
919: Level: intermediate
921: .seealso: FNView(), FNCreate()
922: @*/
923: PetscErrorCode FNViewFromOptions(FN fn,PetscObject obj,const char name[])924: {
926: PetscObjectViewFromOptions((PetscObject)fn,obj,name);
927: PetscFunctionReturn(0);
928: }
930: /*@
931: FNDuplicate - Duplicates a math function, copying all parameters, possibly with a
932: different communicator.
934: Collective on fn
936: Input Parameters:
937: + fn - the math function context
938: - comm - MPI communicator
940: Output Parameter:
941: . newfn - location to put the new FN context
943: Note:
944: In order to use the same MPI communicator as in the original object,
945: use PetscObjectComm((PetscObject)fn).
947: Level: developer
949: .seealso: FNCreate()
950: @*/
951: PetscErrorCode FNDuplicate(FN fn,MPI_Comm comm,FN *newfn)952: {
953: FNType type;
954: PetscScalar alpha,beta;
955: PetscInt meth;
956: FNParallelType ptype;
961: FNCreate(comm,newfn);
962: FNGetType(fn,&type);
963: FNSetType(*newfn,type);
964: FNGetScale(fn,&alpha,&beta);
965: FNSetScale(*newfn,alpha,beta);
966: FNGetMethod(fn,&meth);
967: FNSetMethod(*newfn,meth);
968: FNGetParallel(fn,&ptype);
969: FNSetParallel(*newfn,ptype);
970: if (fn->ops->duplicate) (*fn->ops->duplicate)(fn,comm,newfn);
971: PetscFunctionReturn(0);
972: }
974: /*@C
975: FNDestroy - Destroys FN context that was created with FNCreate().
977: Collective on fn
979: Input Parameter:
980: . fn - the math function context
982: Level: beginner
984: .seealso: FNCreate()
985: @*/
986: PetscErrorCode FNDestroy(FN *fn)987: {
988: PetscInt i;
990: if (!*fn) PetscFunctionReturn(0);
992: if (--((PetscObject)(*fn))->refct > 0) { *fn = 0; PetscFunctionReturn(0); }
993: if ((*fn)->ops->destroy) (*(*fn)->ops->destroy)(*fn);
994: for (i=0;i<(*fn)->nw;i++) MatDestroy(&(*fn)->W[i]);
995: PetscHeaderDestroy(fn);
996: PetscFunctionReturn(0);
997: }
999: /*@C
1000: FNRegister - Adds a mathematical function to the FN package.
1002: Not collective
1004: Input Parameters:
1005: + name - name of a new user-defined FN1006: - function - routine to create context
1008: Notes:
1009: FNRegister() may be called multiple times to add several user-defined functions.
1011: Level: advanced
1013: .seealso: FNRegisterAll()
1014: @*/
1015: PetscErrorCode FNRegister(const char *name,PetscErrorCode (*function)(FN))1016: {
1017: FNInitializePackage();
1018: PetscFunctionListAdd(&FNList,name,function);
1019: PetscFunctionReturn(0);
1020: }