Mercurial > octave

comparison libinterp/corefcn/quadcc.cc @ 15995:906ae1705522

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
quadcc.cc: Use heap instead of stack for large temporary variables. * libinterp/corefcn/quadcc.cc: Use heap instead of stack for large temporary variables. Use Octave coding conventions. Place debug code under #if.
author Rik <rik@octave.org>
date Mon, 04 Feb 2013 11:10:14 -0800
parents 352ef2fb3ad1
children 1aff8c38c53c
comparison
equal deleted inserted replaced
15994:352ef2fb3ad1 15995:906ae1705522
31 #include "defun.h" 31 #include "defun.h"
32 #include "error.h" 32 #include "error.h"
33 #include "oct-obj.h" 33 #include "oct-obj.h"
34 #include "utils.h" 34 #include "utils.h"
35 35
36 //#include "oct.h" 36
37 //#include "defun.h" 37 /* Extended debugging */
38 38 #define DEBUG_QUADCC 0
39 /* Define the size of the interval heap. */ 39
40 #define cquad_heapsize 50 40 /* Define the minimum size of the interval heap. */
41 #define min_cquad_heapsize 200
41 42
42 43
43 /* Data of a single interval */ 44 /* Data of a single interval */
44 typedef struct 45 typedef struct
45 { 46 {
1550 static const int idx[4] = { 0, 5, 14, 31 }; 1551 static const int idx[4] = { 0, 5, 14, 31 };
1551 static const double w = M_SQRT2 / 2; 1552 static const double w = M_SQRT2 / 2;
1552 static const int ndiv_max = 20; 1553 static const int ndiv_max = 20;
1553 1554
1554 /* The interval heap. */ 1555 /* The interval heap. */
1555 cquad_ival ivals[cquad_heapsize]; 1556 cquad_ival *ivals;
1556 int heap[cquad_heapsize]; 1557 int *heap;
1557 1558
1558 /* Arguments left and right */ 1559 /* Arguments left and right */
1559 int nargin = args.length (); 1560 int nargin = args.length ();
1560 octave_function *fcn; 1561 octave_function *fcn;
1561 double a, b, tol, iivals[cquad_heapsize], *sing; 1562 double a, b, tol, *iivals, *sing;
1562 1563
1563 /* Variables needed for transforming the integrand. */ 1564 /* Variables needed for transforming the integrand. */
1564 bool wrap = false; 1565 bool wrap = false;
1565 double xw; 1566 double xw;
1566 1567
1586 1587
1587 if (args(0).is_function_handle () || args(0).is_inline_function ()) 1588 if (args(0).is_function_handle () || args(0).is_inline_function ())
1588 fcn = args(0).function_value (); 1589 fcn = args(0).function_value ();
1589 else 1590 else
1590 { 1591 {
1591 std::string fcn_name = unique_symbol_name ("__quadcc_fcn_"); 1592 std::string fcn_name = unique_symbol_name ("__quadcc_fcn__");
1592 std::string fname = "function y = "; 1593 std::string fname = "function y = ";
1593 fname.append (fcn_name); 1594 fname.append (fcn_name);
1594 fname.append ("(x) y = "); 1595 fname.append ("(x) y = ");
1595 fcn = extract_function (args(0), "quadcc", fcn_name, fname, 1596 fcn = extract_function (args(0), "quadcc", fcn_name, fname,
1596 "; endfunction"); 1597 "; endfunction");
1597 } 1598 }
1598 1599
1599 if (!args(1).is_real_scalar ()) 1600 if (! args(1).is_real_scalar ())
1600 { 1601 {
1601 error ("quadcc: lower limit of integration (A) must be a single real scalar"); 1602 error ("quadcc: lower limit of integration (A) must be a single real scalar");
1602 return retval; 1603 return retval;
1603 } 1604 }
1604 else 1605 else
1605 a = args(1).double_value (); 1606 a = args(1).double_value ();
1606 1607
1607 if (!args(2).is_real_scalar ()) 1608 if (! args(2).is_real_scalar ())
1608 { 1609 {
1609 error ("quadcc: upper limit of integration (B) must be a single real scalar"); 1610 error ("quadcc: upper limit of integration (B) must be a single real scalar");
1610 return retval; 1611 return retval;
1611 } 1612 }
1612 else 1613 else
1613 b = args(2).double_value (); 1614 b = args(2).double_value ();
1614 1615
1615 if (nargin < 4 || args(3).is_empty ()) 1616 if (nargin < 4 || args(3).is_empty ())
1616 tol = 1.0e-6; 1617 tol = 1.0e-6;
1617 else if (!args(3).is_real_scalar () || args(3).double_value () <= 0) 1618 else if (! args(3).is_real_scalar () || args(3).double_value () <= 0)
1618 { 1619 {
1619 error ("quadcc: tolerance (TOL) must be a single real scalar > 0"); 1620 error ("quadcc: tolerance (TOL) must be a single real scalar > 0");
1620 return retval; 1621 return retval;
1621 } 1622 }
1622 else 1623 else
1623 tol = args(3).double_value (); 1624 tol = args(3).double_value ();
1624 1625
1625 if (nargin < 5) 1626 if (nargin < 5)
1626 { 1627 {
1627 nivals = 1; 1628 nivals = 1;
1628 iivals[0] = a;
1629 iivals[1] = b;
1630 } 1629 }
1631 else if (!(args(4).is_real_scalar () || args(4).is_real_matrix ())) 1630 else if (!(args(4).is_real_scalar () || args(4).is_real_matrix ()))
1632 { 1631 {
1633 error ("quadcc: list of singularities (SING) must be a vector of real values"); 1632 error ("quadcc: list of singularities (SING) must be a vector of real values");
1634 return retval; 1633 return retval;
1635 } 1634 }
1636 else 1635 else
1637 { 1636 {
1638 nivals = 1 + args(4).length (); 1637 nivals = 1 + args(4).numel ();
1639 if (nivals > cquad_heapsize) 1638 }
1640 { 1639
1641 error ("quadcc: maximum number of singular points is limited to %i", 1640 int cquad_heapsize = (nivals >= min_cquad_heapsize ? nivals + 1
1642 cquad_heapsize-1); 1641 : min_cquad_heapsize);
1643 return retval; 1642 heap = new int [cquad_heapsize];
1644 } 1643 ivals = new cquad_ival [cquad_heapsize];
1644 iivals = new double [cquad_heapsize];
1645
1646 if (nivals == 1)
1647 {
1648 iivals[0] = a;
1649 iivals[1] = b;
1650 }
1651 else
1652 {
1653 // Intervals around singularities
1645 sing = args(4).array_value ().fortran_vec (); 1654 sing = args(4).array_value ().fortran_vec ();
1646 iivals[0] = a; 1655 iivals[0] = a;
1647 for (i = 0; i < nivals - 2; i++) 1656 for (i = 0; i < nivals - 1; i++)
1648 iivals[i + 1] = sing[i]; 1657 iivals[i + 1] = sing[i];
1649 iivals[nivals] = b; 1658 iivals[nivals] = b;
1650 } 1659 }
1651 1660
1652 /* If a or b are +/-Inf, transform the integral. */ 1661 /* If a or b are +/-Inf, transform the integral. */
1653 if (xisinf (a) || xisinf (b)) 1662 if (xisinf (a) || xisinf (b))
1654 { 1663 {
1655 wrap = true; 1664 wrap = true;
1656 for (i = 0; i <= nivals; i++) 1665 for (i = 0; i < nivals + 1; i++)
1657 if (xisinf (iivals[i])) 1666 if (xisinf (iivals[i]))
1658 iivals[i] = gnulib::copysign (1.0, iivals[i]); 1667 iivals[i] = gnulib::copysign (1.0, iivals[i]);
1659 else 1668 else
1660 iivals[i] = 2.0 * atan (iivals[i]) / M_PI; 1669 iivals[i] = 2.0 * atan (iivals[i]) / M_PI;
1661 } 1670 }
1662 1671
1663 1672
1664 /* Initialize the heaps. */ 1673 /* Initialize the heaps. */
1665 for (i = 0; i < cquad_heapsize; i++) 1674 for (i = 0; i < cquad_heapsize; i++)
1666 heap[i] = i; 1675 heap[i] = i;
1667
1668 1676
1669 /* Create the first interval(s). */ 1677 /* Create the first interval(s). */
1670 igral = 0.0; 1678 igral = 0.0;
1671 err = 0.0; 1679 err = 0.0;
1672 for (j = 0; j < nivals; j++) 1680 for (j = 0; j < nivals; j++)
1679 nnans = 0; 1687 nnans = 0;
1680 ColumnVector ex (33); 1688 ColumnVector ex (33);
1681 if (wrap) 1689 if (wrap)
1682 { 1690 {
1683 for (i = 0; i <= n[3]; i++) 1691 for (i = 0; i <= n[3]; i++)
1684 ex (i) = tan (M_PI / 2 * (m + xi[i] * h)); 1692 ex(i) = tan (M_PI / 2 * (m + xi[i] * h));
1685 } 1693 }
1686 else 1694 else
1687 { 1695 {
1688 for (i = 0; i <= n[3]; i++) 1696 for (i = 0; i <= n[3]; i++)
1689 ex (i) = m + xi[i] * h; 1697 ex(i) = m + xi[i] * h;
1690 } 1698 }
1691 fargs(0) = ex; 1699 fargs(0) = ex;
1692 fvals = feval (fcn, fargs, 1); 1700 fvals = feval (fcn, fargs, 1);
1693 if (fvals.length () != 1 || !fvals(0).is_real_matrix ()) 1701 if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1694 { 1702 {
1695 error ("quadcc: integrand F must return a single, real-valued vector"); 1703 error ("quadcc: integrand F must return a single, real-valued vector");
1696 return retval; 1704 return retval;
1697 } 1705 }
1698 Matrix effex = fvals(0).matrix_value (); 1706 Matrix effex = fvals(0).matrix_value ();
1701 error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input"); 1709 error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1702 return retval; 1710 return retval;
1703 } 1711 }
1704 for (i = 0; i <= n[3]; i++) 1712 for (i = 0; i <= n[3]; i++)
1705 { 1713 {
1706 iv->fx[i] = effex (i); 1714 iv->fx[i] = effex(i);
1707 if (wrap) 1715 if (wrap)
1708 { 1716 {
1709 xw = ex(i); 1717 xw = ex(i);
1710 iv->fx[i] *= (1.0 + xw * xw) * M_PI / 2; 1718 iv->fx[i] *= (1.0 + xw * xw) * M_PI / 2;
1711 } 1719 }
1712 neval++; 1720 neval++;
1713 if (!xfinite (iv->fx[i])) 1721 if (! xfinite (iv->fx[i]))
1714 { 1722 {
1715 nans[nnans++] = i; 1723 nans[nnans++] = i;
1716 iv->fx[i] = 0.0; 1724 iv->fx[i] = 0.0;
1717 } 1725 }
1718 } 1726 }
1780 /* Put our finger on the interval with the largest error. */ 1788 /* Put our finger on the interval with the largest error. */
1781 iv = &(ivals[heap[0]]); 1789 iv = &(ivals[heap[0]]);
1782 m = (iv->a + iv->b) / 2; 1790 m = (iv->a + iv->b) / 2;
1783 h = (iv->b - iv->a) / 2; 1791 h = (iv->b - iv->a) / 2;
1784 1792
1785 /* printf 1793 #if (DEBUG_QUADCC)
1794 printf
1786 ("quadcc: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n", 1795 ("quadcc: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
1787 heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth); 1796 heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
1788 */ 1797 #endif
1798
1789 /* Should we try to increase the degree? */ 1799 /* Should we try to increase the degree? */
1790 if (iv->depth < 3) 1800 if (iv->depth < 3)
1791 { 1801 {
1792 1802
1793 /* Keep tabs on some variables. */ 1803 /* Keep tabs on some variables. */
1797 { 1807 {
1798 ColumnVector ex (n[d] / 2); 1808 ColumnVector ex (n[d] / 2);
1799 if (wrap) 1809 if (wrap)
1800 { 1810 {
1801 for (i = 0; i < n[d] / 2; i++) 1811 for (i = 0; i < n[d] / 2; i++)
1802 ex (i) = 1812 ex(i) = tan (M_PI / 2 * (m + xi[(2 * i + 1) * skip[d]] * h));
1803 tan (M_PI / 2 * (m + xi[(2 * i + 1) * skip[d]] * h));
1804 } 1813 }
1805 else 1814 else
1806 { 1815 {
1807 for (i = 0; i < n[d] / 2; i++) 1816 for (i = 0; i < n[d] / 2; i++)
1808 ex (i) = m + xi[(2 * i + 1) * skip[d]] * h; 1817 ex(i) = m + xi[(2 * i + 1) * skip[d]] * h;
1809 } 1818 }
1810 fargs(0) = ex; 1819 fargs(0) = ex;
1811 fvals = feval (fcn, fargs, 1); 1820 fvals = feval (fcn, fargs, 1);
1812 if (fvals.length () != 1 || !fvals(0).is_real_matrix ()) 1821 if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1813 { 1822 {
1814 error ("quadcc: integrand F must return a single, real-valued vector"); 1823 error ("quadcc: integrand F must return a single, real-valued vector");
1815 return retval; 1824 return retval;
1816 } 1825 }
1817 Matrix effex = fvals(0).matrix_value (); 1826 Matrix effex = fvals(0).matrix_value ();
1822 } 1831 }
1823 neval += effex.length (); 1832 neval += effex.length ();
1824 for (i = 0; i < n[d] / 2; i++) 1833 for (i = 0; i < n[d] / 2; i++)
1825 { 1834 {
1826 j = (2 * i + 1) * skip[d]; 1835 j = (2 * i + 1) * skip[d];
1827 iv->fx[j] = effex (i); 1836 iv->fx[j] = effex(i);
1828 if (wrap) 1837 if (wrap)
1829 { 1838 {
1830 xw = ex(i); 1839 xw = ex(i);
1831 iv->fx[j] *= (1.0 + xw * xw) * M_PI / 2; 1840 iv->fx[j] *= (1.0 + xw * xw) * M_PI / 2;
1832 } 1841 }
1833 } 1842 }
1834 } 1843 }
1835 nnans = 0; 1844 nnans = 0;
1836 for (i = 0; i <= 32; i += skip[d]) 1845 for (i = 0; i <= 32; i += skip[d])
1837 { 1846 {
1838 if (!xfinite (iv->fx[i])) 1847 if (! xfinite (iv->fx[i]))
1839 { 1848 {
1840 nans[nnans++] = i; 1849 nans[nnans++] = i;
1841 iv->fx[i] = 0.0; 1850 iv->fx[i] = 0.0;
1842 } 1851 }
1843 } 1852 }
1886 if ((m + h * xi[0]) >= (m + h * xi[1]) 1895 if ((m + h * xi[0]) >= (m + h * xi[1])
1887 || (m + h * xi[31]) >= (m + h * xi[32]) 1896 || (m + h * xi[31]) >= (m + h * xi[32])
1888 || iv->err < fabs (iv->igral) * std::numeric_limits<double>::epsilon () * 10) 1897 || iv->err < fabs (iv->igral) * std::numeric_limits<double>::epsilon () * 10)
1889 { 1898 {
1890 1899
1891 /* printf 1900 #if (DEBUG_QUADCC)
1901 printf
1892 ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n", 1902 ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
1893 heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, 1903 heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
1894 iv->depth); 1904 iv->depth);
1895 */ 1905 #endif
1906
1896 /* Keep this interval's contribution */ 1907 /* Keep this interval's contribution */
1897 err_final += iv->err; 1908 err_final += iv->err;
1898 igral_final += iv->igral; 1909 igral_final += iv->igral;
1899 /* Swap with the last element on the heap */ 1910 /* Swap with the last element on the heap */
1900 t = heap[nivals - 1]; 1911 t = heap[nivals - 1];
1903 nivals--; 1914 nivals--;
1904 /* Fix up the heap */ 1915 /* Fix up the heap */
1905 i = 0; 1916 i = 0;
1906 while (2 * i + 1 < nivals) 1917 while (2 * i + 1 < nivals)
1907 { 1918 {
1908
1909 /* Get the kids */ 1919 /* Get the kids */
1910 j = 2 * i + 1; 1920 j = 2 * i + 1;
1911 /* If the j+1st entry exists and is larger than the jth, 1921 /* If the j+1st entry exists and is larger than the jth,
1912 use it instead. */ 1922 use it instead. */
1913 if (j + 1 < nivals 1923 if (j + 1 < nivals
1946 { 1956 {
1947 ColumnVector ex (n[0] - 1); 1957 ColumnVector ex (n[0] - 1);
1948 if (wrap) 1958 if (wrap)
1949 { 1959 {
1950 for (i = 0; i < n[0] - 1; i++) 1960 for (i = 0; i < n[0] - 1; i++)
1951 ex (i) = tan (M_PI / 2 * (ml + xi[(i + 1) * skip[0]] * hl)); 1961 ex(i) = tan (M_PI / 2 * (ml + xi[(i + 1) * skip[0]] * hl));
1952 } 1962 }
1953 else 1963 else
1954 { 1964 {
1955 for (i = 0; i < n[0] - 1; i++) 1965 for (i = 0; i < n[0] - 1; i++)
1956 ex (i) = ml + xi[(i + 1) * skip[0]] * hl; 1966 ex(i) = ml + xi[(i + 1) * skip[0]] * hl;
1957 } 1967 }
1958 fargs(0) = ex; 1968 fargs(0) = ex;
1959 fvals = feval (fcn, fargs, 1); 1969 fvals = feval (fcn, fargs, 1);
1960 if (fvals.length () != 1 || !fvals(0).is_real_matrix ()) 1970 if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1961 { 1971 {
1962 error ("quadcc: integrand F must return a single, real-valued vector"); 1972 error ("quadcc: integrand F must return a single, real-valued vector");
1963 return retval; 1973 return retval;
1964 } 1974 }
1965 Matrix effex = fvals(0).matrix_value (); 1975 Matrix effex = fvals(0).matrix_value ();
1970 } 1980 }
1971 neval += effex.length (); 1981 neval += effex.length ();
1972 for (i = 0; i < n[0] - 1; i++) 1982 for (i = 0; i < n[0] - 1; i++)
1973 { 1983 {
1974 j = (i + 1) * skip[0]; 1984 j = (i + 1) * skip[0];
1975 ivl->fx[j] = effex (i); 1985 ivl->fx[j] = effex(i);
1976 if (wrap) 1986 if (wrap)
1977 { 1987 {
1978 xw = ex(i); 1988 xw = ex(i);
1979 ivl->fx[j] *= (1.0 + xw * xw) * M_PI / 2; 1989 ivl->fx[j] *= (1.0 + xw * xw) * M_PI / 2;
1980 } 1990 }
1981 } 1991 }
1982 } 1992 }
1983 nnans = 0; 1993 nnans = 0;
1984 for (i = 0; i <= 32; i += skip[0]) 1994 for (i = 0; i <= 32; i += skip[0])
1985 { 1995 {
1986 if (!xfinite (ivl->fx[i])) 1996 if (! xfinite (ivl->fx[i]))
1987 { 1997 {
1988 nans[nnans++] = i; 1998 nans[nnans++] = i;
1989 ivl->fx[i] = 0.0; 1999 ivl->fx[i] = 0.0;
1990 } 2000 }
1991 } 2001 }
2042 { 2052 {
2043 ColumnVector ex (n[0] - 1); 2053 ColumnVector ex (n[0] - 1);
2044 if (wrap) 2054 if (wrap)
2045 { 2055 {
2046 for (i = 0; i < n[0] - 1; i++) 2056 for (i = 0; i < n[0] - 1; i++)
2047 ex (i) = tan (M_PI / 2 * (mr + xi[(i + 1) * skip[0]] * hr)); 2057 ex(i) = tan (M_PI / 2 * (mr + xi[(i + 1) * skip[0]] * hr));
2048 } 2058 }
2049 else 2059 else
2050 { 2060 {
2051 for (i = 0; i < n[0] - 1; i++) 2061 for (i = 0; i < n[0] - 1; i++)
2052 ex (i) = mr + xi[(i + 1) * skip[0]] * hr; 2062 ex(i) = mr + xi[(i + 1) * skip[0]] * hr;
2053 } 2063 }
2054 fargs(0) = ex; 2064 fargs(0) = ex;
2055 fvals = feval (fcn, fargs, 1); 2065 fvals = feval (fcn, fargs, 1);
2056 if (fvals.length () != 1 || !fvals(0).is_real_matrix ()) 2066 if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
2057 { 2067 {
2058 error ("quadcc: integrand F must return a single, real-valued vector"); 2068 error ("quadcc: integrand F must return a single, real-valued vector");
2059 return retval; 2069 return retval;
2060 } 2070 }
2061 Matrix effex = fvals(0).matrix_value (); 2071 Matrix effex = fvals(0).matrix_value ();
2066 } 2076 }
2067 neval += effex.length (); 2077 neval += effex.length ();
2068 for (i = 0; i < n[0] - 1; i++) 2078 for (i = 0; i < n[0] - 1; i++)
2069 { 2079 {
2070 j = (i + 1) * skip[0]; 2080 j = (i + 1) * skip[0];
2071 ivr->fx[j] = effex (i); 2081 ivr->fx[j] = effex(i);
2072 if (wrap) 2082 if (wrap)
2073 { 2083 {
2074 xw = ex(i); 2084 xw = ex(i);
2075 ivr->fx[j] *= (1.0 + xw * xw) * M_PI / 2; 2085 ivr->fx[j] *= (1.0 + xw * xw) * M_PI / 2;
2076 } 2086 }
2077 } 2087 }
2078 } 2088 }
2079 nnans = 0; 2089 nnans = 0;
2080 for (i = 0; i <= 32; i += skip[0]) 2090 for (i = 0; i <= 32; i += skip[0])
2081 { 2091 {
2082 if (!xfinite (ivr->fx[i])) 2092 if (! xfinite (ivr->fx[i]))
2083 { 2093 {
2084 nans[nnans++] = i; 2094 nans[nnans++] = i;
2085 ivr->fx[i] = 0.0; 2095 ivr->fx[i] = 0.0;
2086 } 2096 }
2087 } 2097 }
2193 } 2203 }
2194 } 2204 }
2195 2205
2196 } 2206 }
2197 2207
2198 /* If the heap is about to overflow, remove the last two 2208 /* If the heap is about to overflow, remove the last two intervals. */
2199 intervals. */
2200 while (nivals > cquad_heapsize - 2) 2209 while (nivals > cquad_heapsize - 2)
2201 { 2210 {
2202 iv = &(ivals[heap[nivals - 1]]); 2211 iv = &(ivals[heap[nivals - 1]]);
2203 /* printf 2212 #if (DEBUG_QUADCC)
2213 printf
2204 ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n", 2214 ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
2205 heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, 2215 heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
2206 iv->depth); 2216 iv->depth);
2207 */ 2217 #endif
2208 err_final += iv->err; 2218 err_final += iv->err;
2209 igral_final += iv->igral; 2219 igral_final += iv->igral;
2210 nivals--; 2220 nivals--;
2211 } 2221 }
2212 2222
2220 } 2230 }
2221 2231
2222 } 2232 }
2223 2233
2224 /* Dump the contents of the heap. */ 2234 /* Dump the contents of the heap. */
2225 /* for (i = 0; i < nivals; i++) 2235 #if (DEBUG_QUADCC)
2236 for (i = 0; i < nivals; i++)
2226 { 2237 {
2227 iv = &(ivals[heap[i]]); 2238 iv = &(ivals[heap[i]]);
2228 printf 2239 printf
2229 ("quadcc: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i, ndiv=%i\n", 2240 ("quadcc: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i, ndiv=%i\n",
2230 i, heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth, 2241 i, heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth,
2231 iv->rdepth, iv->ndiv); 2242 iv->rdepth, iv->ndiv);
2232 } 2243 }
2233 */ 2244 #endif
2245
2246 /* Clean up heap memory */
2247 delete [] heap;
2248 delete [] ivals;
2249 delete [] iivals;
2250
2234 /* Clean up and present the results. */ 2251 /* Clean up and present the results. */
2235 if (nargout > 2) 2252 if (nargout > 2)
2236 retval(2) = neval; 2253 retval(2) = neval;
2237 if (nargout > 1) 2254 if (nargout > 1)
2238 retval(1) = err; 2255 retval(1) = err;
2263 %! y(idx < 1e-10) = NaN; 2280 %! y(idx < 1e-10) = NaN;
2264 %!endfunction 2281 %!endfunction
2265 2282
2266 %!test 2283 %!test
2267 %! [q, err, npoints] = quadcc ("__nansin", -pi, pi); 2284 %! [q, err, npoints] = quadcc ("__nansin", -pi, pi);
2268 %! assert (q, 0, eps); 2285 %! assert (q, 0, 1e-6);
2269 %! assert (err, 0, 15*eps); 2286 %! assert (err, 0, 15*eps);
2270 2287
2271 %% Test input validation 2288 %% Test input validation
2272 %!error (quadcc ()) 2289 %!error (quadcc ())
2273 %!error (quadcc (@sin)) 2290 %!error (quadcc (@sin))