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 { 7 const PetscReal x = a[0]; 8 *val = x*PetscLogReal(1+x); 9 } 10 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 } 16 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 } 22 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 } 29 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 } 36 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 } 42 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 } 49 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 } 56 57 static void func9(const PetscReal x[], void *dummy, PetscReal *val) 58 { 59 *val = PetscLogReal(PetscCosReal(x[0])); 60 } 61 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 } 69 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 } 75 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 } 83 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 } 91 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 } 99 100 int main(int argc, char **argv) 101 { 102 #if PETSC_SCALAR_SIZE == 32 103 PetscInt digits = 7; 104 #elif PETSC_SCALAR_SIZE == 64 105 PetscInt digits = 14; 106 #else 107 PetscInt digits = 14; 108 #endif 109 /* for some reason in __float128 precision it cannot get more accuracy for some of the integrals */ 110 #if defined(PETSC_USE_REAL___FLOAT128) 111 const PetscReal epsilon = 2.2204460492503131e-16; 112 #else 113 const PetscReal epsilon = 2500.*PETSC_MACHINE_EPSILON; 114 #endif 115 const PetscReal bounds[28] = 116 { 117 0.0, 1.0, 118 0.0, 1.0, 119 0.0, PETSC_PI/2., 120 0.0, 1.0, 121 0.0, 1.0, 122 0.0, 1.0, 123 0.0, 1.0, 124 0.0, 1.0, 125 0.0, PETSC_PI/2., 126 0.0, PETSC_PI/2., 127 0.0, 1.0, 128 0.0, 1.0, 129 0.0, 1.0, 130 0.0, 1.0 131 }; 132 const PetscReal analytic[14] = 133 { 134 0.250000000000000, 135 0.210657251225806988108092302182988001695680805674, 136 1.905238690482675827736517833351916563195085437332, 137 0.514041895890070761397629739576882871630921844127, 138 -.444444444444444444444444444444444444444444444444, 139 0.785398163397448309615660845819875721049292349843, 140 1.198140234735592207439922492280323878227212663216, 141 2.000000000000000000000000000000000000000000000000, 142 -1.08879304515180106525034444911880697366929185018, 143 2.221441469079183123507940495030346849307310844687, 144 1.570796326794896619231321691639751442098584699687, 145 1.772453850905516027298167483341145182797549456122, 146 1.253314137315500251207882642405522626503493370304, 147 0.500000000000000000000000000000000000000000000000 148 }; 149 void (*funcs[14])(const PetscReal[], void *, PetscReal *) = {func1, func2, func3, func4, func5, func6, func7, func8, func9, func10, func11, func12, func13, func14}; 150 PetscInt f; 151 PetscErrorCode ierr; 152 153 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 154 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Test Options","none");CHKERRQ(ierr); 155 ierr = PetscOptionsBoundedInt("-digits", "The number of significant digits for the integral","ex3.c",digits,&digits,NULL,1);CHKERRQ(ierr); 156 ierr = PetscOptionsEnd();CHKERRQ(ierr); 157 158 /* Integrate each function */ 159 for (f = 0; f < 14; ++f) { 160 PetscReal integral; 161 162 /* These can only be integrated accuractely using MPFR */ 163 if ((f == 6) || (f == 7) || (f == 9) || (f == 11)) continue; 164 #if defined(PETSC_USE_REAL_SINGLE) 165 if (f == 8) continue; 166 #endif 167 ierr = PetscDTTanhSinhIntegrate(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, NULL, &integral);CHKERRQ(ierr); 168 if (PetscAbsReal(integral - analytic[f]) > PetscMax(epsilon, PetscPowRealInt(10.0, -digits)) || PetscIsInfOrNanScalar(integral - analytic[f])) { 169 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); 170 } 171 } 172 #if defined(PETSC_HAVE_MPFR) 173 for (f = 0; f < 14; ++f) { 174 PetscReal integral; 175 176 ierr = PetscDTTanhSinhIntegrateMPFR(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, NULL, &integral);CHKERRQ(ierr); 177 if (PetscAbsReal(integral - analytic[f]) > PetscPowRealInt(10.0, -digits)) { 178 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); 179 } 180 } 181 #endif 182 ierr = PetscFinalize(); 183 return ierr; 184 } 185 186 /*TEST 187 test: 188 suffix: 0 189 TEST*/ 190