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