This function applies the adaptive integration algorithm QAGS taking account of the user-supplied locations of singular points. The array pts of length npts should contain the endpoints of the integration ranges defined by the integration region and locations of the singularities. For example, to integrate over the region @math{(a,b)} with break-points at @math{x_1, x_2, x_3} (where @math{a < x_1 < x_2 < x_3 < b}) the following pts array should be used
pts[0] = a pts[1] = x_1 pts[2] = x_2 pts[3] = x_3 pts[4] = b
with npts = 5.
If you know the locations of the singular points in the integration
region then this routine will be faster than QAGS
.