File: converged.h

package info (click to toggle)
cubature 1.0.4%2Bds-1
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, sid
  • size: 268 kB
  • sloc: ansic: 1,513; makefile: 80; sh: 34
file content (73 lines) | stat: -rw-r--r-- 2,399 bytes parent folder | download | duplicates (6)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
/* Body of convergence test, shared between hcubature.c and
   pcubature.c.  We use an #include file because the two routines use
   somewhat different data structures, and define macros ERR(j) and
   VAL(j) to get the error and value estimates, respectively, for
   integrand j. */
{
     unsigned j;
#    define SQR(x) ((x) * (x))
     switch (norm) {
	 case ERROR_INDIVIDUAL:
	      for (j = 0; j < fdim; ++j)
		   if (ERR(j) > reqAbsError && ERR(j) > fabs(VAL(j))*reqRelError)
			return 0;
	      return 1;
	      
	 case ERROR_PAIRED:
	      for (j = 0; j+1 < fdim; j += 2) {
		   double maxerr, serr, err, maxval, sval, val;
		   /* scale to avoid overflow/underflow */
		   maxerr = ERR(j) > ERR(j+1) ? ERR(j) : ERR(j+1);
		   maxval = VAL(j) > VAL(j+1) ? VAL(j) : VAL(j+1);
		   serr = maxerr > 0 ? 1/maxerr : 1;
		   sval = maxval > 0 ? 1/maxval : 1;
		   err = sqrt(SQR(ERR(j)*serr) + SQR(ERR(j+1)*serr)) * maxerr;
		   val = sqrt(SQR(VAL(j)*sval) + SQR(VAL(j+1)*sval)) * maxval;
		   if (err > reqAbsError && err > val*reqRelError)
			return 0;
	      }
	      if (j < fdim) /* fdim is odd, do last dimension individually */
		   if (ERR(j) > reqAbsError && ERR(j) > fabs(VAL(j))*reqRelError)
			return 0;
	      return 1;

	 case ERROR_L1: {
	      double err = 0, val = 0;
	      for (j = 0; j < fdim; ++j) {
		   err += ERR(j);
		   val += fabs(VAL(j));
	      }
	      return err <= reqAbsError || err <= val*reqRelError;
	 }

	 case ERROR_LINF: {
	      double err = 0, val = 0;
	      for (j = 0; j < fdim; ++j) {
		   double absval = fabs(VAL(j));
		   if (ERR(j) > err) err = ERR(j);
		   if (absval > val) val = absval;
	      }
	      return err <= reqAbsError || err <= val*reqRelError;
	 }

	 case ERROR_L2: {
	      double maxerr = 0, maxval = 0, serr, sval, err = 0, val = 0;
	      /* scale values by 1/max to avoid overflow/underflow */
	      for (j = 0; j < fdim; ++j) {
		   double absval = fabs(VAL(j));
		   if (ERR(j) > maxerr) maxerr = ERR(j);
		   if (absval > maxval) maxval = absval;
	      }
	      serr = maxerr > 0 ? 1/maxerr : 1;
	      sval = maxval > 0 ? 1/maxval : 1;
	      for (j = 0; j < fdim; ++j) {
		   err += SQR(ERR(j) * serr);
		   val += SQR(fabs(VAL(j)) * sval);
	      }
	      err = sqrt(err) * maxerr;
	      val = sqrt(val) * maxval;
	      return err <= reqAbsError || err <= val*reqRelError;
	 }
     }
     return 1; /* unreachable */
}