Actual source code: ex3.c
1: static char help[] = "Tests quadrature.\n\n";
3: #include <petscdt.h>
5: static void func1(const PetscReal a[], void *dummy, PetscReal *val)
6: {
7: const PetscReal x = a[0];
8: *val = x*PetscLogReal(1+x);
9: }
11: static void func2(const PetscReal a[], void *dummy, PetscReal *val)
12: {
13: const PetscReal x = a[0];
14: *val = x*x*PetscAtanReal(x);
15: }
17: static void func3(const PetscReal a[], void *dummy, PetscReal *val)
18: {
19: const PetscReal x = a[0];
20: *val = PetscExpReal(x)*PetscCosReal(x);
21: }
23: static void func4(const PetscReal a[], void *dummy, 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: }
30: static void func5(const PetscReal a[], void *dummy, 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: }
37: static void func6(const PetscReal a[], void *dummy, PetscReal *val)
38: {
39: const PetscReal x = a[0];
40: *val = PetscSqrtReal(1-x*x);
41: }
43: static void func7(const PetscReal a[], void *dummy, 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: }
50: static void func8(const PetscReal a[], void *dummy, 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: }
57: static void func9(const PetscReal x[], void *dummy, PetscReal *val)
58: {
59: *val = PetscLogReal(PetscCosReal(x[0]));
60: }
62: static void func10(const PetscReal a[], void *dummy, 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: }
70: static void func11(const PetscReal a[], void *dummy, PetscReal *val)
71: {
72: const PetscReal x = a[0];
73: *val = 1/(1-2*x+2*x*x);
74: }
76: static void func12(const PetscReal a[], void *dummy, 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: }
84: static void func13(const PetscReal a[], void *dummy, 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: }
92: static void func14(const PetscReal a[], void *dummy, 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: }
100: int main(int argc, char **argv)
101: {
102: #if PETSC_SCALAR_SIZE == 32
103: PetscInt digits = 7;
104: #elif PETSC_SCALAR_SIZE == 64
105: PetscInt digits = 14;
106: #else
107: PetscInt digits = 14;
108: #endif
109: /* for some reason in __float128 precision it cannot get more accuracy for some of the integrals */
110: #if defined(PETSC_USE_REAL___FLOAT128)
111: const PetscReal epsilon = 2.2204460492503131e-16;
112: #else
113: const PetscReal epsilon = 2500.*PETSC_MACHINE_EPSILON;
114: #endif
115: const PetscReal bounds[28] =
116: {
117: 0.0, 1.0,
118: 0.0, 1.0,
119: 0.0, PETSC_PI/2.,
120: 0.0, 1.0,
121: 0.0, 1.0,
122: 0.0, 1.0,
123: 0.0, 1.0,
124: 0.0, 1.0,
125: 0.0, PETSC_PI/2.,
126: 0.0, PETSC_PI/2.,
127: 0.0, 1.0,
128: 0.0, 1.0,
129: 0.0, 1.0,
130: 0.0, 1.0
131: };
132: const PetscReal analytic[14] =
133: {
134: 0.250000000000000,
135: 0.210657251225806988108092302182988001695680805674,
136: 1.905238690482675827736517833351916563195085437332,
137: 0.514041895890070761397629739576882871630921844127,
138: -.444444444444444444444444444444444444444444444444,
139: 0.785398163397448309615660845819875721049292349843,
140: 1.198140234735592207439922492280323878227212663216,
141: 2.000000000000000000000000000000000000000000000000,
142: -1.08879304515180106525034444911880697366929185018,
143: 2.221441469079183123507940495030346849307310844687,
144: 1.570796326794896619231321691639751442098584699687,
145: 1.772453850905516027298167483341145182797549456122,
146: 1.253314137315500251207882642405522626503493370304,
147: 0.500000000000000000000000000000000000000000000000
148: };
149: void (*funcs[14])(const PetscReal[], void *, PetscReal *) = {func1, func2, func3, func4, func5, func6, func7, func8, func9, func10, func11, func12, func13, func14};
150: PetscInt f;
151: PetscErrorCode ierr;
153: PetscInitialize(&argc, &argv, NULL, help);
154: PetscOptionsBegin(PETSC_COMM_WORLD,"","Test Options","none");
155: PetscOptionsBoundedInt("-digits", "The number of significant digits for the integral","ex3.c",digits,&digits,NULL,1);
156: PetscOptionsEnd();
158: /* Integrate each function */
159: for (f = 0; f < 14; ++f) {
160: PetscReal integral;
162: /* These can only be integrated accuractely using MPFR */
163: if ((f == 6) || (f == 7) || (f == 9) || (f == 11)) continue;
164: #if defined(PETSC_USE_REAL_SINGLE)
165: if (f == 8) continue;
166: #endif
167: PetscDTTanhSinhIntegrate(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, NULL, &integral);
168: if (PetscAbsReal(integral - analytic[f]) > PetscMax(epsilon, PetscPowRealInt(10.0, -digits)) || PetscIsInfOrNanScalar(integral - analytic[f])) {
169: PetscPrintf(PETSC_COMM_SELF, "The integral of func%2d is wrong: %g (%g)\n", f+1, (double)integral, (double) PetscAbsReal(integral - analytic[f]));
170: }
171: }
172: #if defined(PETSC_HAVE_MPFR)
173: for (f = 0; f < 14; ++f) {
174: PetscReal integral;
176: PetscDTTanhSinhIntegrateMPFR(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, NULL, &integral);
177: if (PetscAbsReal(integral - analytic[f]) > PetscPowRealInt(10.0, -digits)) {
178: PetscPrintf(PETSC_COMM_SELF, "The integral of func%2d is wrong: %g (%g)\n", f+1, (double)integral, (double)PetscAbsReal(integral - analytic[f]));
179: }
180: }
181: #endif
182: PetscFinalize();
183: return 0;
184: }
186: /*TEST
187: test:
188: suffix: 0
189: TEST*/