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 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;

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();

124:   /* Integrate each function */
125:   for (f = 0; f < 14; ++f) {
126:     PetscReal integral;

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;

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: }

150: /*TEST
151:   build:
152:     requires: !complex

154:   test:
155:     suffix: 0
156: TEST*/