#ifndef EIGEN_LAURIEGAUTSCHI_H
#define EIGEN_LAURIEGAUTSCHI_H
namespace Eigen
{
/**
* \ingroup NumericalIntegration_Module
*
* \class LaurieGautschi
*
* \tparam Scalar floating point type
*
* \brief This class computes Kronrod abscissae & weights for arbitrary precision.
*
* Based on work of Dirk Laurie and Walter Gautschi.
* D. P. Laurie (1997). Calculation of Gauss-Kronrod Quadrature Rules.
* Mathematics of Computation, 66(219), 1133-1145
* Created by Pavel Holoborodko, November 7, 2011.
* Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder, and Matt Beall September 2014
*
* \TODO Ensure only appropriates types are used for Scalar, e.g. prohibit integers.
*/
template <typename Scalar>
class LaurieGautschi
{
public:
/**
* \brief Recurrence coefficients for monic Jacobi polynomials.
*
* This method generates the first N recurrence
* coefficients for monic Jacobi polynomials with parameters
* alpha and beta. These are orthogonal on [-1,1] relative to the
* weight function w(t)=(1-t)^a(1+t)^b. The N alpha-coefficients
* are stored in \a alphaOut, the n beta-coefficients in \a betaOut.
* http://en.wikipedia.org/wiki/Jacobi_polynomials
*
* Created by Dirk Laurie, 6-22-1998; edited by Walter Gautschi, 4-4-2002.
* Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder and Matt Beall, September 2014
*
* \param[in] N Number of recurrence coefficients
* \param[in] alpha Alpha parameter of the Jacobi-polynomials
* \param[in] beta Beta parameter of the Jacobi-polynomials
* \param[in,out] alphaOut N alpha-coefficients
* \param[in,out] betaOut N beta-coefficients
*/
static void r_jacobi(const Index N,
const Scalar alpha,
const Scalar beta,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& alphaOut,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& betaOut)
{
using std::pow;
using std::tgamma;
//TODO : make use the eigen assert facilities
assert(alpha > Scalar(-1));
assert(beta > Scalar(-1));
assert(alphaOut.rows() == betaOut.rows());
assert(alphaOut.rows() > 0);
assert(N <= alphaOut.rows());
alphaOut(0) = (beta-alpha)/(alpha+beta+Scalar(2.));
betaOut(0) = pow(Scalar(2.),(alpha+beta+Scalar(1.))) * tgamma(alpha+Scalar(1.)) *
tgamma(beta+Scalar(1.)) / tgamma(alpha+beta+Scalar(2.));
for (Index n = 1; n < N; ++n)
{
Scalar nAlphaBeta = Scalar(2.) * n + alpha + beta;
alphaOut(n) = (beta*beta - alpha*alpha) / (nAlphaBeta * (nAlphaBeta+Scalar(2.)));
betaOut(n) = Scalar(4.) * (n+alpha)*(n+beta)*n*(n+alpha+beta) /
(nAlphaBeta*nAlphaBeta*(nAlphaBeta+Scalar(1.)) * (nAlphaBeta-Scalar(1.)));
}
}
/**
* \brief Recurrence coefficients for monic Jacobi polynomials on [0,1].
*
* This method generates the first N recurrence
* coefficients for monic Jacobi polynomials on [0,1] with
* parameters alpha and beta. These are orthogonal on [0,1] relative
* to the weight function w(t)=(1-t)^alpha t^beta. The N alpha-
* coefficients are stored in \a alphaOut, the N beta-
* coefficients in \a betaOut.
* http://en.wikipedia.org/wiki/Jacobi_polynomials
*
* Created by Dirk Laurie, 6-22-1998; edited by Walter Gautschi, 4-4-2002.
* Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder and Matt Beall, September 2014
*
* \param[in] N Number of recurrence coefficients
* \param[in] alpha Alpha parameter of the Jacobi-polynomials
* \param[in] beta Beta parameter of the Jacobi-polynomials
* \param[in,out] alphaOut N alpha-coefficients
* \param[in,out] betaOut N beta-coefficients
*/
static void r_jacobi_01(const Index N,
const Scalar alpha,
const Scalar beta,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& alphaOut,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& betaOut)
{
using std::pow;
//TODO : make use the eigen assert facilities
assert(alpha > Scalar(-1));
assert(beta > Scalar(-1));
assert(alphaOut.rows() == betaOut.rows());
assert(alphaOut.rows() > 0);
assert(N <= alphaOut.rows());
r_jacobi(N, alpha, beta, alphaOut, betaOut);
for (Index n = 0; n < N; ++n)
{
alphaOut(n) = (Scalar(1.)+alphaOut(n)) / Scalar(2.);
}
betaOut(0) = betaOut(0) / pow(Scalar(2), alpha + beta + Scalar(1.));
for (Index n = 1; n < N; ++n)
{
betaOut(n) = betaOut(n) / Scalar(4.);
}
}
/**
* \brief Jacobi-Kronrod matrix.
*
* This method produces the alpha- and beta-elements in
* the Jacobi-Kronrod matrix of order 2N+1 for the weight
* function (or measure) w. The input data for the weight
* function w are the recurrence coefficients of the associated
* orthogonal polynomials, which are stored in \a alphaIn and \a betaIn .
* At least ceil(3*N/2)+1 coefficients should be provided.
* The 2N+1 alpha- and beta-elements are returned in \a alpha and \a beta
* respectively.
*
* Created by Dirk Laurie, 6-22.1998
* Edited by Pavel Holoborodko, November 7, 2011
* Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder and Matt Beall, September 2014
*
* \param[in] N Number of nodes
* \param[in,out] alphaIn The recurrence coefficients of the associated orthogonal polynomials
* \param[in,out] betaIn The recurrence coefficients of the associated orthogonal polynomials
* \param[in,out] alpha Alpha-elements in the Jacobi-Kronrod matrix of order 2N+1
* \param[in,out] beta Beta-elements in the Jacobi-Kronrod matrix of order 2N+1
*
*/
static void r_kronrod(const Index N,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& alphaIn,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& betaIn,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& alpha,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& beta)
{
using std::ceil;
using std::floor;
//TODO : make use the eigen assert facilities
assert(alphaIn.rows() == betaIn.rows());
assert(alphaIn.rows() >= ceil(3*N/2) + 1 );
assert(alphaIn.rows() > 0);
assert(alpha.rows() == 2*N+1);
assert(beta.rows() == 2*N+1);
alpha = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N+1);
beta = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N+1);
for (Index k = 0; k <= (Index)floor(3*N/2) + 1; ++k)
{
alpha(k) = alphaIn(k);
}
for (Index k=0; k <= (Index)ceil(3*N/2) + 1; ++k)
{
beta(k) = betaIn(k);
}
Eigen::Array<Scalar, Eigen::Dynamic, 1> sigma = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero((Index)floor(N/2) + 2);
Eigen::Array<Scalar, Eigen::Dynamic, 1> tempVector = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero((Index)floor(N/2) + 2);
tempVector(1) = beta(N+1);
for (Index m = 0; m < N-2+1; ++m)
{
Scalar u = 0;
for(Index k = (Index)floor((m+1) / 2); k >= 0;--k)
{
Index l = m-k;
u = u + (alpha(k+N+1)-alpha(l)) * tempVector(k+1) + beta(k+N+1)*sigma(k) - beta(l)*sigma(k+1);
sigma(k+1) = u;
}
Eigen::Array<Scalar, Eigen::Dynamic, 1> swap = sigma;
sigma = tempVector;
tempVector = swap;
}
for (Index j = (Index)floor(N/2); j>=0; --j)
{
sigma(j+1) = sigma(j);
}
for (Index m = N-1; m < 2*N-3+1; ++m)
{
Index k = m+1-N;
Index j = 0;
Scalar u = 0;
for (k = m+1-N; k < (Index)floor((m-1) / 2) + 1; ++k)
{
Index l = m-k;
j = N-1-l;
u = u - (alpha(k+N+1)-alpha(l))*tempVector(j+1) - beta(k+N+1)*sigma(j+1) + beta(l)*sigma(j+2);
sigma(j+1) = u;
}
k = (Index)floor((m+1) / 2);
if (m % 2 == 0)
{
alpha(k+N+1) = alpha(k) + (sigma(j+1)-beta(k+N+1)*sigma(j+2)) / tempVector(j+2);
}
else
{
beta(k+N+1) = sigma(j+1) / sigma(j+2);
}
Eigen::Array<Scalar, Eigen::Dynamic, 1> swap = sigma;
sigma = tempVector;
tempVector = swap;
}
alpha(2*N) = alpha(N-1)-beta(2*N)*sigma(1) / tempVector(1);
}
/**
* \brief Gauss-Kronrod quadrature formula.
*
* This method generates the (2N+1)-point Gauss-Kronrod
* quadrature rule for the weight function w encoded by the
* recurrence matrix (alpha,beta) of order [ceil(3*n/2)+1]x2 containing
* in its first and second column respectively the alpha- and
* beta-coefficients in the three-term recurrence relation
* for w. The 2N+1 nodes, in increasing order, are output
* into \a nodes, the corresponding weights into \a weights.
*
* Created by Dirk Laurie, June 22, 1998.
* Edited by Pavel Holoborodko, November 7, 2011:
* Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder
* and Matt Beall, September 2014
*
* \param[in] N Number of nodes
* \param[in,out] alpha 2N alpha coefficients (input)
* \param[in,out] beta 2N beta coefficients (input)
* \param[in,out] nodes 2N+1 nodes
* \param[in,out] weights 2N+1 weights corresponding to \a nodes
*/
static void kronrod(const Index N,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& alpha,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& beta,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& nodes,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& weights)
{
using std::abs;
using std::sqrt;
//TODO : make use the eigen assert facilities
assert(N>0);
assert(alpha.rows() == 2*N);
assert(alpha.rows() == beta.rows());
assert(nodes.rows() == 2*N+1);
assert(nodes.rows() == weights.rows());
Eigen::Array<Scalar, Eigen::Dynamic, 1> alpha0 = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N+1);
Eigen::Array<Scalar, Eigen::Dynamic, 1> beta0 = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N+1);
r_kronrod(N, alpha, beta, alpha0, beta0);
// \TODO : CHECK NEEDED LIKE THE ONE ON LINE 21 IN KRONROD.M
// Do we have an approximately equal function in Eigen?
assert(abs(beta0.sum() - (Scalar) (2*N+1)) > 1e-5);
Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> J = Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic>::Zero(2*N+1, 2*N+1);
for (Index k = 0; k < 2*N; ++k)
{
J(k,k) = alpha0(k);
J(k,k+1) = sqrt(beta0(k+1));
J(k+1,k) = J(k,k+1);
}
J(2*N,2*N) = alpha0(2*N);
//TODO : Is this assumption of positive definiteness correct?
SelfAdjointEigenSolver< Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> > es(J);
//TODO : make use the eigen assert facilities
assert(es.info() == Eigen::Success);
nodes = es.eigenvalues();
Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> V = es.eigenvectors();
weights = beta0(0) * (V.row(0).array() * V.row(0).array()).matrix();
}
/**
* \brief Gauss quadrature rule.
*
* Given a weight function w encoded by (alpha,beta) of the
* first N recurrence coefficients for the associated orthogonal
* polynomials, the first column of (alpha,beta) containing the N alpha-
* coefficients and the second column the N beta-coefficients,
* the method generates the nodes and weights of
* the N-point Gauss quadrature rule for the weight function.
* The Nodes, in increasing order, are stored in \a nodes ,
* the N corresponding weights are stored in \a weights .
*
* \param[in] N Number of nodes
* \param[in,out] alpha 2N alpha coefficients (input)
* \param[in,out] beta 2N beta coefficients (input)
* \param[in,out] nodes 2N+1 nodes
* \param[in,out] weights 2N+1 weights corresponding to \a nodes
*/
static void gauss(const Index N,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& alpha,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& beta,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& nodes,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& weights)
{
using std::sqrt;
//TODO : make use the eigen assert facilities
assert(N > 0);
assert(alpha.rows() == 2*N);
assert(alpha.rows() == beta.rows());
assert(nodes.rows() == N);
assert(nodes.rows() == weights.rows());
Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> J = Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic>::Zero(N,N);
J(0,0) = alpha(0);
for (Index n=1; n<N; ++n)
{
J(n,n) = alpha(n);
J(n,n-1) = sqrt(beta(n));
J(n-1,n) = J(n,n-1);
}
//TODO : Is this assumption of positive definiteness correct?
SelfAdjointEigenSolver< Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> > es(J);
//TODO : make use the eigen assert facilities
assert(es.info() == Eigen::Success);
nodes = es.eigenvalues();
Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> V = es.eigenvectors();
weights = beta(0)*(V.row(0).array()*V.row(0).array()).matrix();
}
/**
* \brief Arbitrary precision Kronrod abscissae & weights.
*
* This method computes Kronrod points for (-1,1) with any required precision.
* The result is a vector 2N+1 nodes (N points on either side of zero and zero)
* and the corresponding weights.
*
* \param[in] N Number of nodes
* \param[in,out] nodes Returns a vector of 2N+1 nodes
* \param[in,out] w Returns a vector of weights corresponding to \a nodes
*/
static void mpkronrod(const Index N,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& nodes,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& weights)
{
//TODO : make use the eigen assert facilities
assert(nodes.rows() == 2*N+1);
assert(weights.rows() == 2*N+1);
assert(N>0);
Eigen::Array<Scalar, Eigen::Dynamic, 1> alpha = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N);
Eigen::Array<Scalar, Eigen::Dynamic, 1> beta = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N);
r_jacobi_01( 2*N, Scalar(0), Scalar(0), alpha, beta);
kronrod(N,alpha,beta,nodes,weights);
for (Index i=0; i<nodes.rows(); ++i)
{
nodes(i) = Scalar(2.)*nodes(i) - Scalar(1.);
weights(i) = Scalar(2.)*weights(i);
}
}
/**
* \brief Arbitrary precision Gauss abscissae & weights.
*
* This method computes Kronrod points for (-1,1) with any required precision.
* The result is a vector of N nodes the corresponding weights.
*
* \param[in] N Number of nodes
* \param[in,out] nodes Returns a vector of 2N+1 nodes
* \param[in,out] weights Returns a vector of weights corresponding to \a nodes
*/
static void mpgauss(const Index N,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& nodes,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& weights)
{
//TODO : make use the eigen assert facilities
assert(nodes.rows() == N);
assert(weights.rows() == N);
assert(N > 0);
Eigen::Array<Scalar, Eigen::Dynamic, 1> alpha = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N);
Eigen::Array<Scalar, Eigen::Dynamic, 1> beta = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N);
r_jacobi_01(2*N, Scalar(0), Scalar(0), alpha, beta);
gauss(N, alpha, beta, nodes, weights);
for (Index i = 0;i<nodes.rows();++i)
{
nodes(i) = Scalar(2.)*nodes(i) - Scalar(1.);
weights(i) = Scalar(2.)*weights(i);
}
}
/**
* \brief A method ofr computing Gauss/Kronrod abscissae and weights
* \param[in] nNodes Gauss-Legendre degree
* \param[out] abscGaussKronrod Gauss/Kronrod abscissae
* \param[out] weightGaussKronrod Gauss/Kronrod weights
* \param[out] abscGauss Gauss abscissae
* \param[out] weightGauss Gauss weights
*/
static void computeAbscissaeAndWeights(unsigned int nNodes,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& abscGaussKronrod,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& weightGaussKronrod,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& abscGauss,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& weightGauss)
{
using std::abs;
Eigen::Array<Scalar, Eigen::Dynamic, 1> xGK = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*nNodes+1);
Eigen::Array<Scalar, Eigen::Dynamic, 1> wGK = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*nNodes+1);
Eigen::Array<Scalar, Eigen::Dynamic, 1> xG = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(nNodes);
Eigen::Array<Scalar, Eigen::Dynamic, 1> wG = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(nNodes);
LaurieGautschi::mpkronrod(nNodes,xGK,wGK);
LaurieGautschi::mpgauss(nNodes,xG,wG);
unsigned int arraySize = nNodes + 1;
abscGaussKronrod = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize);
weightGaussKronrod = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize);
abscGauss = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize/2);
weightGauss = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize/2);
for (unsigned int i=0; i<arraySize; ++i)
{
abscGaussKronrod(i) = abs(xGK(i));
weightGaussKronrod(i) = wGK(i);
}
abscGaussKronrod(arraySize-1) = Scalar(0);
for (unsigned int i=0; i<arraySize/2; ++i)
{
abscGauss(i) = abs(xG(i));
weightGauss(i) = wG(i);
}
}
};
}//namespace Eigen
#endif //EIGEN_LAURIEGAUTSCHI_H