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