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

153:   PetscInitialize(&argc, &argv, NULL, help);
154:   PetscOptionsBegin(PETSC_COMM_WORLD,"","Test Options","none");
155:   PetscOptionsBoundedInt("-digits", "The number of significant digits for the integral","ex3.c",digits,&digits,NULL,1);
156:   PetscOptionsEnd();

158:   /* Integrate each function */
159:   for (f = 0; f < 14; ++f) {
160:     PetscReal integral;

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:     PetscDTTanhSinhIntegrate(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, NULL, &integral);
168:     if (PetscAbsReal(integral - analytic[f]) > PetscMax(epsilon, PetscPowRealInt(10.0, -digits)) || PetscIsInfOrNanScalar(integral - analytic[f])) {
169:       PetscPrintf(PETSC_COMM_SELF, "The integral of func%2d is wrong: %g (%g)\n", f+1, (double)integral, (double) PetscAbsReal(integral - analytic[f]));
170:     }
171:   }
172: #if defined(PETSC_HAVE_MPFR)
173:   for (f = 0; f < 14; ++f) {
174:     PetscReal integral;

176:     PetscDTTanhSinhIntegrateMPFR(funcs[f], bounds[f*2+0], bounds[f*2+1], digits, NULL, &integral);
177:     if (PetscAbsReal(integral - analytic[f]) > PetscPowRealInt(10.0, -digits)) {
178:       PetscPrintf(PETSC_COMM_SELF, "The integral of func%2d is wrong: %g (%g)\n", f+1, (double)integral, (double)PetscAbsReal(integral - analytic[f]));
179:     }
180:   }
181: #endif
182:   PetscFinalize();
183:   return 0;
184: }

186: /*TEST
187:   test:
188:     suffix: 0
189: TEST*/