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 : Rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials
12 : */
13 :
14 : #include <slepc/private/fnimpl.h> /*I "slepcfn.h" I*/
15 :
16 : typedef struct {
17 : PetscScalar *pcoeff; /* numerator coefficients */
18 : PetscInt np; /* length of array pcoeff, p(x) has degree np-1 */
19 : PetscScalar *qcoeff; /* denominator coefficients */
20 : PetscInt nq; /* length of array qcoeff, q(x) has degree nq-1 */
21 : } FN_RATIONAL;
22 :
23 11392 : static PetscErrorCode FNEvaluateFunction_Rational(FN fn,PetscScalar x,PetscScalar *y)
24 : {
25 11392 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
26 11392 : PetscInt i;
27 11392 : PetscScalar p,q;
28 :
29 11392 : PetscFunctionBegin;
30 11392 : if (!ctx->np) p = 1.0;
31 : else {
32 11391 : p = ctx->pcoeff[0];
33 18560 : for (i=1;i<ctx->np;i++)
34 7169 : p = ctx->pcoeff[i]+x*p;
35 : }
36 11392 : if (!ctx->nq) *y = p;
37 : else {
38 1242 : q = ctx->qcoeff[0];
39 2572 : for (i=1;i<ctx->nq;i++)
40 1330 : q = ctx->qcoeff[i]+x*q;
41 1242 : PetscCheck(q!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Function not defined in the requested value");
42 1242 : *y = p/q;
43 : }
44 11392 : PetscFunctionReturn(PETSC_SUCCESS);
45 : }
46 :
47 : /*
48 : Horner evaluation of P=p(A)
49 : d = degree of polynomial; coeff = coefficients of polynomial; W = workspace
50 : */
51 1513 : static PetscErrorCode EvaluatePoly(Mat A,Mat P,Mat W,PetscInt d,PetscScalar *coeff)
52 : {
53 1513 : PetscInt j;
54 :
55 1513 : PetscFunctionBegin;
56 1513 : PetscCall(MatZeroEntries(P));
57 1513 : if (!d) PetscCall(MatShift(P,1.0));
58 : else {
59 1513 : PetscCall(MatShift(P,coeff[0]));
60 2541 : for (j=1;j<d;j++) {
61 1028 : PetscCall(MatMatMult(P,A,MAT_REUSE_MATRIX,PETSC_DEFAULT,&W));
62 1028 : PetscCall(MatCopy(W,P,SAME_NONZERO_PATTERN));
63 1028 : PetscCall(MatShift(P,coeff[j]));
64 : }
65 : }
66 1513 : PetscFunctionReturn(PETSC_SUCCESS);
67 : }
68 :
69 1242 : static PetscErrorCode FNEvaluateFunctionMat_Rational(FN fn,Mat A,Mat B)
70 : {
71 1242 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
72 1242 : Mat P,Q,W,F;
73 1242 : PetscBool iscuda;
74 :
75 1242 : PetscFunctionBegin;
76 1242 : if (A==B) PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&P));
77 1240 : else P = B;
78 1242 : PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&W));
79 :
80 1242 : PetscCall(EvaluatePoly(A,P,W,ctx->np,ctx->pcoeff));
81 1242 : if (ctx->nq) {
82 245 : PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Q));
83 245 : PetscCall(EvaluatePoly(A,Q,W,ctx->nq,ctx->qcoeff));
84 245 : PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
85 490 : PetscCall(MatGetFactor(Q,iscuda?MATSOLVERCUDA:MATSOLVERPETSC,MAT_FACTOR_LU,&F));
86 245 : PetscCall(MatLUFactorSymbolic(F,Q,NULL,NULL,NULL));
87 245 : PetscCall(MatLUFactorNumeric(F,Q,NULL));
88 245 : PetscCall(MatMatSolve(F,P,P));
89 245 : PetscCall(MatDestroy(&F));
90 245 : PetscCall(MatDestroy(&Q));
91 : }
92 :
93 1242 : if (A==B) {
94 2 : PetscCall(MatCopy(P,B,SAME_NONZERO_PATTERN));
95 2 : PetscCall(MatDestroy(&P));
96 : }
97 1242 : PetscCall(MatDestroy(&W));
98 1242 : PetscFunctionReturn(PETSC_SUCCESS);
99 : }
100 :
101 24 : static PetscErrorCode FNEvaluateFunctionMatVec_Rational(FN fn,Mat A,Vec v)
102 : {
103 24 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
104 24 : Mat P,Q,W,F;
105 24 : Vec b;
106 24 : PetscBool iscuda;
107 :
108 24 : PetscFunctionBegin;
109 24 : PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&P));
110 24 : PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&W));
111 :
112 24 : PetscCall(EvaluatePoly(A,P,W,ctx->np,ctx->pcoeff));
113 24 : if (ctx->nq) {
114 2 : PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Q));
115 2 : PetscCall(EvaluatePoly(A,Q,W,ctx->nq,ctx->qcoeff));
116 2 : PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSECUDA,&iscuda));
117 4 : PetscCall(MatGetFactor(Q,iscuda?MATSOLVERCUDA:MATSOLVERPETSC,MAT_FACTOR_LU,&F));
118 2 : PetscCall(MatLUFactorSymbolic(F,Q,NULL,NULL,NULL));
119 2 : PetscCall(MatLUFactorNumeric(F,Q,NULL));
120 2 : PetscCall(MatCreateVecs(P,&b,NULL));
121 2 : PetscCall(MatGetColumnVector(P,b,0));
122 2 : PetscCall(MatSolve(F,b,v));
123 2 : PetscCall(VecDestroy(&b));
124 2 : PetscCall(MatDestroy(&F));
125 2 : PetscCall(MatDestroy(&Q));
126 22 : } else PetscCall(MatGetColumnVector(P,v,0));
127 :
128 24 : PetscCall(MatDestroy(&P));
129 24 : PetscCall(MatDestroy(&W));
130 24 : PetscFunctionReturn(PETSC_SUCCESS);
131 : }
132 :
133 3863 : static PetscErrorCode FNEvaluateDerivative_Rational(FN fn,PetscScalar x,PetscScalar *yp)
134 : {
135 3863 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
136 3863 : PetscInt i;
137 3863 : PetscScalar p,q,pp,qp;
138 :
139 3863 : PetscFunctionBegin;
140 3863 : if (!ctx->np) {
141 : p = 1.0;
142 : pp = 0.0;
143 : } else {
144 3862 : p = ctx->pcoeff[0];
145 3862 : pp = 0.0;
146 6164 : for (i=1;i<ctx->np;i++) {
147 2302 : pp = p+x*pp;
148 2302 : p = ctx->pcoeff[i]+x*p;
149 : }
150 : }
151 3863 : if (!ctx->nq) *yp = pp;
152 : else {
153 459 : q = ctx->qcoeff[0];
154 459 : qp = 0.0;
155 922 : for (i=1;i<ctx->nq;i++) {
156 463 : qp = q+x*qp;
157 463 : q = ctx->qcoeff[i]+x*q;
158 : }
159 459 : PetscCheck(q!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
160 459 : *yp = (pp*q-p*qp)/(q*q);
161 : }
162 3863 : PetscFunctionReturn(PETSC_SUCCESS);
163 : }
164 :
165 20 : static PetscErrorCode FNView_Rational(FN fn,PetscViewer viewer)
166 : {
167 20 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
168 20 : PetscBool isascii;
169 20 : PetscInt i;
170 20 : char str[50];
171 :
172 20 : PetscFunctionBegin;
173 20 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
174 20 : if (isascii) {
175 20 : if (fn->alpha!=(PetscScalar)1.0 || fn->beta!=(PetscScalar)1.0) {
176 2 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_FALSE));
177 2 : PetscCall(PetscViewerASCIIPrintf(viewer," scale factors: alpha=%s,",str));
178 2 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
179 2 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_FALSE));
180 2 : PetscCall(PetscViewerASCIIPrintf(viewer," beta=%s\n",str));
181 2 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
182 : }
183 20 : if (!ctx->nq) {
184 13 : if (!ctx->np) PetscCall(PetscViewerASCIIPrintf(viewer," constant: 1.0\n"));
185 13 : else if (ctx->np==1) {
186 5 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[0],PETSC_FALSE));
187 5 : PetscCall(PetscViewerASCIIPrintf(viewer," constant: %s\n",str));
188 : } else {
189 8 : PetscCall(PetscViewerASCIIPrintf(viewer," polynomial: "));
190 8 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
191 21 : for (i=0;i<ctx->np-1;i++) {
192 13 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[i],PETSC_TRUE));
193 13 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s*x^%1" PetscInt_FMT,str,ctx->np-i-1));
194 : }
195 8 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[ctx->np-1],PETSC_TRUE));
196 8 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s\n",str));
197 8 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
198 : }
199 7 : } else if (!ctx->np) {
200 1 : PetscCall(PetscViewerASCIIPrintf(viewer," inverse polynomial: 1 / ("));
201 1 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
202 3 : for (i=0;i<ctx->nq-1;i++) {
203 2 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->qcoeff[i],PETSC_TRUE));
204 2 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s*x^%1" PetscInt_FMT,str,ctx->nq-i-1));
205 : }
206 1 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->qcoeff[ctx->nq-1],PETSC_TRUE));
207 1 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s)\n",str));
208 1 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
209 : } else {
210 6 : PetscCall(PetscViewerASCIIPrintf(viewer," rational function: ("));
211 6 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
212 12 : for (i=0;i<ctx->np-1;i++) {
213 6 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[i],PETSC_TRUE));
214 6 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s*x^%1" PetscInt_FMT,str,ctx->np-i-1));
215 : }
216 6 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->pcoeff[ctx->np-1],PETSC_TRUE));
217 6 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s) / (",str));
218 18 : for (i=0;i<ctx->nq-1;i++) {
219 12 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->qcoeff[i],PETSC_TRUE));
220 12 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s*x^%1" PetscInt_FMT,str,ctx->nq-i-1));
221 : }
222 6 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->qcoeff[ctx->nq-1],PETSC_TRUE));
223 6 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s)\n",str));
224 6 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
225 : }
226 : }
227 20 : PetscFunctionReturn(PETSC_SUCCESS);
228 : }
229 :
230 241 : static PetscErrorCode FNRationalSetNumerator_Rational(FN fn,PetscInt np,PetscScalar *pcoeff)
231 : {
232 241 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
233 241 : PetscInt i;
234 :
235 241 : PetscFunctionBegin;
236 241 : PetscCheck(np>=0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Argument np cannot be negative");
237 241 : ctx->np = np;
238 241 : PetscCall(PetscFree(ctx->pcoeff));
239 241 : if (np) {
240 240 : PetscCall(PetscMalloc1(np,&ctx->pcoeff));
241 626 : for (i=0;i<np;i++) ctx->pcoeff[i] = pcoeff[i];
242 : }
243 241 : PetscFunctionReturn(PETSC_SUCCESS);
244 : }
245 :
246 : /*@
247 : FNRationalSetNumerator - Sets the parameters defining the numerator of the
248 : rational function.
249 :
250 : Logically Collective
251 :
252 : Input Parameters:
253 : + fn - the math function context
254 : . np - number of coefficients
255 : - pcoeff - coefficients (array of scalar values)
256 :
257 : Notes:
258 : Let the rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials.
259 : This function provides the coefficients of the numerator p(x).
260 : Hence, p(x) is of degree np-1.
261 : If np is zero, then the numerator is assumed to be p(x)=1.
262 :
263 : In polynomials, high order coefficients are stored in the first positions
264 : of the array, e.g. to represent x^2-3 use {1,0,-3}.
265 :
266 : Level: intermediate
267 :
268 : .seealso: FNRationalSetDenominator(), FNRationalGetNumerator()
269 : @*/
270 241 : PetscErrorCode FNRationalSetNumerator(FN fn,PetscInt np,PetscScalar pcoeff[])
271 : {
272 241 : PetscFunctionBegin;
273 241 : PetscValidHeaderSpecific(fn,FN_CLASSID,1);
274 723 : PetscValidLogicalCollectiveInt(fn,np,2);
275 241 : if (np) PetscAssertPointer(pcoeff,3);
276 241 : PetscTryMethod(fn,"FNRationalSetNumerator_C",(FN,PetscInt,PetscScalar*),(fn,np,pcoeff));
277 241 : PetscFunctionReturn(PETSC_SUCCESS);
278 : }
279 :
280 1 : static PetscErrorCode FNRationalGetNumerator_Rational(FN fn,PetscInt *np,PetscScalar *pcoeff[])
281 : {
282 1 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
283 1 : PetscInt i;
284 :
285 1 : PetscFunctionBegin;
286 1 : if (np) *np = ctx->np;
287 1 : if (pcoeff) {
288 1 : if (!ctx->np) *pcoeff = NULL;
289 : else {
290 1 : PetscCall(PetscMalloc1(ctx->np,pcoeff));
291 3 : for (i=0;i<ctx->np;i++) (*pcoeff)[i] = ctx->pcoeff[i];
292 : }
293 : }
294 1 : PetscFunctionReturn(PETSC_SUCCESS);
295 : }
296 :
297 : /*@C
298 : FNRationalGetNumerator - Gets the parameters that define the numerator of the
299 : rational function.
300 :
301 : Not Collective
302 :
303 : Input Parameter:
304 : . fn - the math function context
305 :
306 : Output Parameters:
307 : + np - number of coefficients
308 : - pcoeff - coefficients (array of scalar values, length nq)
309 :
310 : Notes:
311 : The values passed by user with FNRationalSetNumerator() are returned (or null
312 : pointers otherwise).
313 : The pcoeff array should be freed by the user when no longer needed.
314 :
315 : Level: intermediate
316 :
317 : .seealso: FNRationalSetNumerator()
318 : @*/
319 1 : PetscErrorCode FNRationalGetNumerator(FN fn,PetscInt *np,PetscScalar *pcoeff[])
320 : {
321 1 : PetscFunctionBegin;
322 1 : PetscValidHeaderSpecific(fn,FN_CLASSID,1);
323 1 : PetscUseMethod(fn,"FNRationalGetNumerator_C",(FN,PetscInt*,PetscScalar**),(fn,np,pcoeff));
324 1 : PetscFunctionReturn(PETSC_SUCCESS);
325 : }
326 :
327 33 : static PetscErrorCode FNRationalSetDenominator_Rational(FN fn,PetscInt nq,PetscScalar *qcoeff)
328 : {
329 33 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
330 33 : PetscInt i;
331 :
332 33 : PetscFunctionBegin;
333 33 : PetscCheck(nq>=0,PetscObjectComm((PetscObject)fn),PETSC_ERR_ARG_OUTOFRANGE,"Argument nq cannot be negative");
334 33 : ctx->nq = nq;
335 33 : PetscCall(PetscFree(ctx->qcoeff));
336 33 : if (nq) {
337 32 : PetscCall(PetscMalloc1(nq,&ctx->qcoeff));
338 103 : for (i=0;i<nq;i++) ctx->qcoeff[i] = qcoeff[i];
339 : }
340 33 : PetscFunctionReturn(PETSC_SUCCESS);
341 : }
342 :
343 : /*@
344 : FNRationalSetDenominator - Sets the parameters defining the denominator of the
345 : rational function.
346 :
347 : Logically Collective
348 :
349 : Input Parameters:
350 : + fn - the math function context
351 : . nq - number of coefficients
352 : - qcoeff - coefficients (array of scalar values)
353 :
354 : Notes:
355 : Let the rational function r(x) = p(x)/q(x), where p(x) and q(x) are polynomials.
356 : This function provides the coefficients of the denominator q(x).
357 : Hence, q(x) is of degree nq-1.
358 : If nq is zero, then the function is assumed to be polynomial, r(x) = p(x).
359 :
360 : In polynomials, high order coefficients are stored in the first positions
361 : of the array, e.g. to represent x^2-3 use {1,0,-3}.
362 :
363 : Level: intermediate
364 :
365 : .seealso: FNRationalSetNumerator(), FNRationalGetDenominator()
366 : @*/
367 33 : PetscErrorCode FNRationalSetDenominator(FN fn,PetscInt nq,PetscScalar qcoeff[])
368 : {
369 33 : PetscFunctionBegin;
370 33 : PetscValidHeaderSpecific(fn,FN_CLASSID,1);
371 99 : PetscValidLogicalCollectiveInt(fn,nq,2);
372 33 : if (nq) PetscAssertPointer(qcoeff,3);
373 33 : PetscTryMethod(fn,"FNRationalSetDenominator_C",(FN,PetscInt,PetscScalar*),(fn,nq,qcoeff));
374 33 : PetscFunctionReturn(PETSC_SUCCESS);
375 : }
376 :
377 11 : static PetscErrorCode FNRationalGetDenominator_Rational(FN fn,PetscInt *nq,PetscScalar *qcoeff[])
378 : {
379 11 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
380 11 : PetscInt i;
381 :
382 11 : PetscFunctionBegin;
383 11 : if (nq) *nq = ctx->nq;
384 11 : if (qcoeff) {
385 11 : if (!ctx->nq) *qcoeff = NULL;
386 : else {
387 4 : PetscCall(PetscMalloc1(ctx->nq,qcoeff));
388 13 : for (i=0;i<ctx->nq;i++) (*qcoeff)[i] = ctx->qcoeff[i];
389 : }
390 : }
391 11 : PetscFunctionReturn(PETSC_SUCCESS);
392 : }
393 :
394 : /*@C
395 : FNRationalGetDenominator - Gets the parameters that define the denominator of the
396 : rational function.
397 :
398 : Not Collective
399 :
400 : Input Parameter:
401 : . fn - the math function context
402 :
403 : Output Parameters:
404 : + nq - number of coefficients
405 : - qcoeff - coefficients (array of scalar values, length nq)
406 :
407 : Notes:
408 : The values passed by user with FNRationalSetDenominator() are returned (or a null
409 : pointer otherwise).
410 : The qcoeff array should be freed by the user when no longer needed.
411 :
412 : Level: intermediate
413 :
414 : .seealso: FNRationalSetDenominator()
415 : @*/
416 11 : PetscErrorCode FNRationalGetDenominator(FN fn,PetscInt *nq,PetscScalar *qcoeff[])
417 : {
418 11 : PetscFunctionBegin;
419 11 : PetscValidHeaderSpecific(fn,FN_CLASSID,1);
420 11 : PetscUseMethod(fn,"FNRationalGetDenominator_C",(FN,PetscInt*,PetscScalar**),(fn,nq,qcoeff));
421 11 : PetscFunctionReturn(PETSC_SUCCESS);
422 : }
423 :
424 223 : static PetscErrorCode FNSetFromOptions_Rational(FN fn,PetscOptionItems *PetscOptionsObject)
425 : {
426 : #define PARMAX 10
427 223 : PetscScalar array[PARMAX];
428 223 : PetscInt i,k;
429 223 : PetscBool flg;
430 :
431 223 : PetscFunctionBegin;
432 223 : PetscOptionsHeadBegin(PetscOptionsObject,"FN Rational Options");
433 :
434 223 : k = PARMAX;
435 2453 : for (i=0;i<k;i++) array[i] = 0;
436 223 : PetscCall(PetscOptionsScalarArray("-fn_rational_numerator","Numerator coefficients (one or more scalar values separated with a comma without spaces)","FNRationalSetNumerator",array,&k,&flg));
437 223 : if (flg) PetscCall(FNRationalSetNumerator(fn,k,array));
438 :
439 223 : k = PARMAX;
440 2453 : for (i=0;i<k;i++) array[i] = 0;
441 223 : PetscCall(PetscOptionsScalarArray("-fn_rational_denominator","Denominator coefficients (one or more scalar values separated with a comma without spaces)","FNRationalSetDenominator",array,&k,&flg));
442 223 : if (flg) PetscCall(FNRationalSetDenominator(fn,k,array));
443 :
444 223 : PetscOptionsHeadEnd();
445 223 : PetscFunctionReturn(PETSC_SUCCESS);
446 : }
447 :
448 28 : static PetscErrorCode FNDuplicate_Rational(FN fn,MPI_Comm comm,FN *newfn)
449 : {
450 28 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data,*ctx2 = (FN_RATIONAL*)(*newfn)->data;
451 28 : PetscInt i;
452 :
453 28 : PetscFunctionBegin;
454 28 : ctx2->np = ctx->np;
455 28 : if (ctx->np) {
456 28 : PetscCall(PetscMalloc1(ctx->np,&ctx2->pcoeff));
457 78 : for (i=0;i<ctx->np;i++) ctx2->pcoeff[i] = ctx->pcoeff[i];
458 : }
459 28 : ctx2->nq = ctx->nq;
460 28 : if (ctx->nq) {
461 10 : PetscCall(PetscMalloc1(ctx->nq,&ctx2->qcoeff));
462 32 : for (i=0;i<ctx->nq;i++) ctx2->qcoeff[i] = ctx->qcoeff[i];
463 : }
464 28 : PetscFunctionReturn(PETSC_SUCCESS);
465 : }
466 :
467 267 : static PetscErrorCode FNDestroy_Rational(FN fn)
468 : {
469 267 : FN_RATIONAL *ctx = (FN_RATIONAL*)fn->data;
470 :
471 267 : PetscFunctionBegin;
472 267 : PetscCall(PetscFree(ctx->pcoeff));
473 267 : PetscCall(PetscFree(ctx->qcoeff));
474 267 : PetscCall(PetscFree(fn->data));
475 267 : PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetNumerator_C",NULL));
476 267 : PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetNumerator_C",NULL));
477 267 : PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetDenominator_C",NULL));
478 267 : PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetDenominator_C",NULL));
479 267 : PetscFunctionReturn(PETSC_SUCCESS);
480 : }
481 :
482 267 : SLEPC_EXTERN PetscErrorCode FNCreate_Rational(FN fn)
483 : {
484 267 : FN_RATIONAL *ctx;
485 :
486 267 : PetscFunctionBegin;
487 267 : PetscCall(PetscNew(&ctx));
488 267 : fn->data = (void*)ctx;
489 :
490 267 : fn->ops->evaluatefunction = FNEvaluateFunction_Rational;
491 267 : fn->ops->evaluatederivative = FNEvaluateDerivative_Rational;
492 267 : fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Rational;
493 267 : fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Rational;
494 : #if defined(PETSC_HAVE_CUDA)
495 : fn->ops->evaluatefunctionmatcuda[0] = FNEvaluateFunctionMat_Rational;
496 : fn->ops->evaluatefunctionmatveccuda[0] = FNEvaluateFunctionMatVec_Rational;
497 : #endif
498 267 : fn->ops->setfromoptions = FNSetFromOptions_Rational;
499 267 : fn->ops->view = FNView_Rational;
500 267 : fn->ops->duplicate = FNDuplicate_Rational;
501 267 : fn->ops->destroy = FNDestroy_Rational;
502 267 : PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetNumerator_C",FNRationalSetNumerator_Rational));
503 267 : PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetNumerator_C",FNRationalGetNumerator_Rational));
504 267 : PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNRationalSetDenominator_C",FNRationalSetDenominator_Rational));
505 267 : PetscCall(PetscObjectComposeFunction((PetscObject)fn,"FNRationalGetDenominator_C",FNRationalGetDenominator_Rational));
506 267 : PetscFunctionReturn(PETSC_SUCCESS);
507 : }
|
