1 static char help[] = "Tests quadrature.\n\n"; 2 3 #include <petscdt.h> 4 5 static void func1(PetscReal x, PetscReal *val) 6 { 7 *val = x*PetscLogReal(1+x); 8 } 9 10 static void func2(PetscReal x, PetscReal *val) 11 { 12 *val = x*x*PetscAtanReal(x); 13 } 14 15 static void func3(PetscReal x, PetscReal *val) 16 { 17 *val = PetscExpReal(x)*PetscCosReal(x); 18 } 19 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 } 25 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 } 31 32 static void func6(PetscReal x, PetscReal *val) 33 { 34 *val = PetscSqrtReal(1-x*x); 35 } 36 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 } 42 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 } 48 49 static void func9(PetscReal x, PetscReal *val) 50 { 51 *val = PetscLogReal(PetscCosReal(x)); 52 } 53 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 } 60 61 static void func11(PetscReal x, PetscReal *val) 62 { 63 *val = 1/(1-2*x+2*x*x); 64 } 65 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 } 72 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 } 79 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 } 86 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; 139 140 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 141 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Test Options","none");CHKERRQ(ierr); 142 ierr = PetscOptionsBoundedInt("-digits", "The number of significant digits for the integral","ex3.c",digits,&digits,NULL,1);CHKERRQ(ierr); 143 ierr = PetscOptionsEnd();CHKERRQ(ierr); 144 145 /* Integrate each function */ 146 for (f = 0; f < 14; ++f) { 147 PetscReal integral; 148 149 /* These can only be integrated accuractely using MPFR */ 150 if ((f == 6) || (f == 7) || (f == 9) || (f == 11)) continue; 151 #if defined(PETSC_USE_REAL_SINGLE) 152 if (f == 8) continue; 153 #endif 154 ierr = PetscDTTanhSinhIntegrate(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, &integral);CHKERRQ(ierr); 155 if (PetscAbsReal(integral - analytic[f]) > PetscMax(epsilon, PetscPowRealInt(10.0, -digits)) || PetscIsInfOrNanScalar(integral - analytic[f])) { 156 ierr = PetscPrintf(PETSC_COMM_SELF, "The integral of func%2d is wrong: %g (%g)\n", f+1, (double)integral, (double) PetscAbsReal(integral - analytic[f]));CHKERRQ(ierr); 157 } 158 } 159 #if defined(PETSC_HAVE_MPFR) 160 for (f = 0; f < 14; ++f) { 161 PetscReal integral; 162 163 ierr = PetscDTTanhSinhIntegrateMPFR(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, &integral);CHKERRQ(ierr); 164 if (PetscAbsReal(integral - analytic[f]) > PetscPowRealInt(10.0, -digits)) { 165 ierr = PetscPrintf(PETSC_COMM_SELF, "The integral of func%2d is wrong: %g (%g)\n", f+1, (double)integral, (double)PetscAbsReal(integral - analytic[f]));CHKERRQ(ierr); 166 } 167 } 168 #endif 169 ierr = PetscFinalize(); 170 return ierr; 171 } 172 173 /*TEST 174 test: 175 suffix: 0 176 TEST*/ 177