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