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