Actual source code: ex3.c
petsc-3.7.7 2017-09-25
1: static char help[] = "Tests quadrature.\n\n";
3: #include <petscdt.h>
7: static void func1(PetscReal x, PetscReal *val)
8: {
9: *val = x*PetscLogReal(1+x);
10: }
14: static void func2(PetscReal x, PetscReal *val)
15: {
16: *val = x*x*PetscAtanReal(x);
17: }
21: static void func3(PetscReal x, PetscReal *val)
22: {
23: *val = PetscExpReal(x)*PetscCosReal(x);
24: }
28: static void func4(PetscReal x, PetscReal *val)
29: {
30: const PetscReal u = PetscSqrtReal(2.0 + x*x);
31: *val = PetscAtanReal(u)/((1.0 + x*x)*u);
32: }
36: static void func5(PetscReal x, PetscReal *val)
37: {
38: if (x == 0.0) *val = 0.0;
39: else *val = PetscSqrtReal(x)*PetscLogReal(x);
40: }
44: static void func6(PetscReal x, PetscReal *val)
45: {
46: *val = PetscSqrtReal(1-x*x);
47: }
51: static void func7(PetscReal x, PetscReal *val)
52: {
53: if (x == 1.0) *val = PETSC_INFINITY;
54: else *val = PetscSqrtReal(x)/PetscSqrtReal(1-x*x);
55: }
59: static void func8(PetscReal x, PetscReal *val)
60: {
61: if (x == 0.0) *val = PETSC_INFINITY;
62: else *val = PetscLogReal(x)*PetscLogReal(x);
63: }
67: static void func9(PetscReal x, PetscReal *val)
68: {
69: *val = PetscLogReal(PetscCosReal(x));
70: }
74: static void func10(PetscReal x, PetscReal *val)
75: {
76: if (x == 0.0) *val = 0.0;
77: else if (x == 1.0) *val = PETSC_INFINITY;
78: *val = PetscSqrtReal(PetscTanReal(x));
79: }
82: static void func11(PetscReal x, PetscReal *val)
83: {
84: *val = 1/(1-2*x+2*x*x);
85: }
88: static void func12(PetscReal x, PetscReal *val)
89: {
90: if (x == 0.0) *val = 0.0;
91: else if (x == 1.0) *val = PETSC_INFINITY;
92: else *val = PetscExpReal(1-1/x)/PetscSqrtReal(x*x*x-x*x*x*x);
93: }
96: static void func13(PetscReal x, PetscReal *val)
97: {
98: if (x == 0.0) *val = 0.0;
99: else if (x == 1.0) *val = 1.0;
100: else *val = PetscExpReal(-(1/x-1)*(1/x-1)/2)/(x*x);
101: }
104: static void func14(PetscReal x, PetscReal *val)
105: {
106: if (x == 0.0) *val = 0.0;
107: else if (x == 1.0) *val = 1.0;
108: else *val = PetscExpReal(1-1/x)*PetscCosReal(1/x-1)/(x*x);
109: }
114: int main(int argc, char **argv)
115: {
116: #if PETSC_SCALAR_SIZE == 32
117: PetscInt digits = 7;
118: #elif PETSC_SCALAR_SIZE == 64
119: PetscInt digits = 14;
120: #else
121: PetscInt digits = 14;
122: #endif
123: /* for some reason in __float128 precision it cannot get more accuracy for some of the integrals */
124: #if defined(PETSC_USE_REAL___FLOAT128)
125: const PetscReal epsilon = 2.2204460492503131e-16;
126: #else
127: const PetscReal epsilon = 2500.*PETSC_MACHINE_EPSILON;
128: #endif
129: const PetscReal bounds[28] =
130: {
131: 0.0, 1.0,
132: 0.0, 1.0,
133: 0.0, PETSC_PI/2.,
134: 0.0, 1.0,
135: 0.0, 1.0,
136: 0.0, 1.0,
137: 0.0, 1.0,
138: 0.0, 1.0,
139: 0.0, PETSC_PI/2.,
140: 0.0, PETSC_PI/2.,
141: 0.0, 1.0,
142: 0.0, 1.0,
143: 0.0, 1.0,
144: 0.0, 1.0
145: };
146: const PetscReal analytic[14] =
147: {
148: 0.250000000000000,
149: 0.210657251225806988108092302182988001695680805674,
150: 1.905238690482675827736517833351916563195085437332,
151: 0.514041895890070761397629739576882871630921844127,
152: -.444444444444444444444444444444444444444444444444,
153: 0.785398163397448309615660845819875721049292349843,
154: 1.198140234735592207439922492280323878227212663216,
155: 2.000000000000000000000000000000000000000000000000,
156: -1.08879304515180106525034444911880697366929185018,
157: 2.221441469079183123507940495030346849307310844687,
158: 1.570796326794896619231321691639751442098584699687,
159: 1.772453850905516027298167483341145182797549456122,
160: 1.253314137315500251207882642405522626503493370304,
161: 0.500000000000000000000000000000000000000000000000
162: };
163: void (*funcs[14])(PetscReal, PetscReal *) = {func1, func2, func3, func4, func5, func6, func7, func8, func9, func10, func11, func12, func13, func14};
164: PetscInt f;
165: PetscErrorCode ierr;
167: PetscInitialize(&argc, &argv, PETSC_NULL, help);
169: PetscOptionsBegin(PETSC_COMM_WORLD,"","Test Options","none");
170: PetscOptionsInt("-digits", "The number of significant digits for the integral","ex3.c",digits,&digits,NULL);
171: PetscOptionsEnd();
173: /* Integrate each function */
174: for (f = 0; f < 14; ++f) {
175: PetscReal integral;
177: /* These can only be integrated accuractely using MPFR */
178: if ((f == 6) || (f == 7) || (f == 9) || (f == 11)) continue;
179: #ifdef PETSC_USE_REAL_SINGLE
180: if (f == 8) continue;
181: #endif
182: PetscDTTanhSinhIntegrate(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, &integral);
183: if (PetscAbsReal(integral - analytic[f]) > PetscMax(epsilon, PetscPowRealInt(10.0, -digits)) || PetscIsInfOrNanScalar(integral - analytic[f])) {
184: PetscPrintf(PETSC_COMM_SELF, "The integral of func%2d is wrong: %g (%g)\n", f+1, (double)integral, (double) PetscAbsReal(integral - analytic[f]));
185: }
186: }
187: #ifdef PETSC_HAVE_MPFR
188: for (f = 0; f < 14; ++f) {
189: PetscReal integral;
191: PetscDTTanhSinhIntegrateMPFR(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, &integral);
192: if (PetscAbsReal(integral - analytic[f]) > PetscPowRealInt(10.0, -digits)) {
193: PetscPrintf(PETSC_COMM_SELF, "The integral of func%2d is wrong: %g (%g)\n", f+1, (double)integral, (double)PetscAbsReal(integral - analytic[f]));
194: }
195: }
196: #endif
197: PetscFinalize();
198: return 0;
199: }