1 static char help[] = "Tests quadrature.\n\n";
2
3 #include <petscdt.h>
4
func1(const PetscReal a[],void * unused,PetscReal * val)5 static void func1(const PetscReal a[], void *unused, PetscReal *val)
6 {
7 const PetscReal x = a[0];
8 *val = x * PetscLogReal(1 + x);
9 }
10
func2(const PetscReal a[],void * unused,PetscReal * val)11 static void func2(const PetscReal a[], void *unused, PetscReal *val)
12 {
13 const PetscReal x = a[0];
14 *val = x * x * PetscAtanReal(x);
15 }
16
func3(const PetscReal a[],void * unused,PetscReal * val)17 static void func3(const PetscReal a[], void *unused, PetscReal *val)
18 {
19 const PetscReal x = a[0];
20 *val = PetscExpReal(x) * PetscCosReal(x);
21 }
22
func4(const PetscReal a[],void * unused,PetscReal * val)23 static void func4(const PetscReal a[], void *unused, PetscReal *val)
24 {
25 const PetscReal x = a[0];
26 const PetscReal u = PetscSqrtReal(2.0 + x * x);
27 *val = PetscAtanReal(u) / ((1.0 + x * x) * u);
28 }
29
func5(const PetscReal a[],void * unused,PetscReal * val)30 static void func5(const PetscReal a[], void *unused, PetscReal *val)
31 {
32 const PetscReal x = a[0];
33 if (x == 0.0) *val = 0.0;
34 else *val = PetscSqrtReal(x) * PetscLogReal(x);
35 }
36
func6(const PetscReal a[],void * unused,PetscReal * val)37 static void func6(const PetscReal a[], void *unused, PetscReal *val)
38 {
39 const PetscReal x = a[0];
40 *val = PetscSqrtReal(1 - x * x);
41 }
42
func7(const PetscReal a[],void * unused,PetscReal * val)43 static void func7(const PetscReal a[], void *unused, PetscReal *val)
44 {
45 const PetscReal x = a[0];
46 if (x == 1.0) *val = PETSC_INFINITY;
47 else *val = PetscSqrtReal(x) / PetscSqrtReal(1 - x * x);
48 }
49
func8(const PetscReal a[],void * unused,PetscReal * val)50 static void func8(const PetscReal a[], void *unused, PetscReal *val)
51 {
52 const PetscReal x = a[0];
53 if (x == 0.0) *val = PETSC_INFINITY;
54 else *val = PetscLogReal(x) * PetscLogReal(x);
55 }
56
func9(const PetscReal x[],void * unused,PetscReal * val)57 static void func9(const PetscReal x[], void *unused, PetscReal *val)
58 {
59 *val = PetscLogReal(PetscCosReal(x[0]));
60 }
61
func10(const PetscReal a[],void * unused,PetscReal * val)62 static void func10(const PetscReal a[], void *unused, PetscReal *val)
63 {
64 const PetscReal x = a[0];
65 if (x == 0.0) *val = 0.0;
66 else if (x == 1.0) *val = PETSC_INFINITY;
67 *val = PetscSqrtReal(PetscTanReal(x));
68 }
69
func11(const PetscReal a[],void * unused,PetscReal * val)70 static void func11(const PetscReal a[], void *unused, PetscReal *val)
71 {
72 const PetscReal x = a[0];
73 *val = 1 / (1 - 2 * x + 2 * x * x);
74 }
75
func12(const PetscReal a[],void * unused,PetscReal * val)76 static void func12(const PetscReal a[], void *unused, PetscReal *val)
77 {
78 const PetscReal x = a[0];
79 if (x == 0.0) *val = 0.0;
80 else if (x == 1.0) *val = PETSC_INFINITY;
81 else *val = PetscExpReal(1 - 1 / x) / PetscSqrtReal(x * x * x - x * x * x * x);
82 }
83
func13(const PetscReal a[],void * unused,PetscReal * val)84 static void func13(const PetscReal a[], void *unused, PetscReal *val)
85 {
86 const PetscReal x = a[0];
87 if (x == 0.0) *val = 0.0;
88 else if (x == 1.0) *val = 1.0;
89 else *val = PetscExpReal(-(1 / x - 1) * (1 / x - 1) / 2) / (x * x);
90 }
91
func14(const PetscReal a[],void * unused,PetscReal * val)92 static void func14(const PetscReal a[], void *unused, PetscReal *val)
93 {
94 const PetscReal x = a[0];
95 if (x == 0.0) *val = 0.0;
96 else if (x == 1.0) *val = 1.0;
97 else *val = PetscExpReal(1 - 1 / x) * PetscCosReal(1 / x - 1) / (x * x);
98 }
99
main(int argc,char ** argv)100 int main(int argc, char **argv)
101 {
102 #if defined(PETSC_USE_REAL_SINGLE)
103 PetscInt digits = 7;
104 #else
105 PetscInt digits = 14;
106 #endif
107 /* for some reason in __float128 precision it cannot get more accuracy for some of the integrals */
108 #if defined(PETSC_USE_REAL___FLOAT128)
109 const PetscReal epsilon = 2.2204460492503131e-16;
110 #else
111 const PetscReal epsilon = 2500. * PETSC_MACHINE_EPSILON;
112 #endif
113 const PetscReal bounds[28] = {0.0, 1.0, 0.0, 1.0, 0.0, PETSC_PI / 2., 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, PETSC_PI / 2., 0.0, PETSC_PI / 2., 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0};
114 const PetscReal analytic[14] = {0.250000000000000, 0.210657251225806988108092302182988001695680805674, 1.905238690482675827736517833351916563195085437332, 0.514041895890070761397629739576882871630921844127, -.444444444444444444444444444444444444444444444444, 0.785398163397448309615660845819875721049292349843, 1.198140234735592207439922492280323878227212663216, 2.000000000000000000000000000000000000000000000000, -1.08879304515180106525034444911880697366929185018, 2.221441469079183123507940495030346849307310844687, 1.570796326794896619231321691639751442098584699687, 1.772453850905516027298167483341145182797549456122, 1.253314137315500251207882642405522626503493370304, 0.500000000000000000000000000000000000000000000000};
115 void (*funcs[14])(const PetscReal[], void *, PetscReal *) = {func1, func2, func3, func4, func5, func6, func7, func8, func9, func10, func11, func12, func13, func14};
116 PetscInt f;
117
118 PetscFunctionBeginUser;
119 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
120 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Test Options", "none");
121 PetscCall(PetscOptionsBoundedInt("-digits", "The number of significant digits for the integral", "ex3.c", digits, &digits, NULL, 1));
122 PetscOptionsEnd();
123
124 /* Integrate each function */
125 for (f = 0; f < 14; ++f) {
126 PetscReal integral;
127
128 /* These can only be integrated accuractely using MPFR */
129 if ((f == 6) || (f == 7) || (f == 9) || (f == 11)) continue;
130 #if defined(PETSC_USE_REAL_SINGLE)
131 if (f == 8) continue;
132 #endif
133 PetscCall(PetscDTTanhSinhIntegrate(funcs[f], bounds[f * 2 + 0], bounds[f * 2 + 1], digits, NULL, &integral));
134 if (PetscAbsReal(integral - analytic[f]) > PetscMax(epsilon, PetscPowRealInt(10.0, -digits)) || PetscIsInfOrNanScalar(integral - analytic[f])) {
135 PetscCall(PetscPrintf(PETSC_COMM_SELF, "The integral of func%2" PetscInt_FMT " is wrong: %g (%g)\n", f + 1, (double)integral, (double)PetscAbsReal(integral - analytic[f])));
136 }
137 }
138 #if defined(PETSC_HAVE_MPFR)
139 for (f = 0; f < 14; ++f) {
140 PetscReal integral;
141
142 PetscCall(PetscDTTanhSinhIntegrateMPFR(funcs[f], bounds[f * 2 + 0], bounds[f * 2 + 1], digits, NULL, &integral));
143 if (PetscAbsReal(integral - analytic[f]) > PetscPowRealInt(10.0, -digits)) PetscCall(PetscPrintf(PETSC_COMM_SELF, "The integral of func%2d is wrong: %g (%g)\n", f + 1, (double)integral, (double)PetscAbsReal(integral - analytic[f])));
144 }
145 #endif
146 PetscCall(PetscFinalize());
147 return 0;
148 }
149
150 /*TEST
151 build:
152 requires: !complex
153
154 test:
155 suffix: 0
156 output_file: output/empty.out
157 TEST*/
158