opsRootsFromJacobiMatrix.hh

Go to the documentation of this file.

#ifndef NPSTAT_OPSROOTSFROMJACOBIMATRIX_HH_
#define NPSTAT_OPSROOTSFROMJACOBIMATRIX_HH_
 
/*!
// \file opsRootsFromJacobiMatrix.hh
//
// \brief Roots of orthogonal polynomials calculated from corresponding
//        Jacobi matrices
//
// Author: I. Volobouev
//
// February 2023
*/
 
#include <vector>
#include <cassert>
#include <algorithm>
#include <stdexcept>
 
#include "npstat/nm/Matrix.hh"
 
namespace npstat {
    /** The roots are returned in the increasing order */
    template<class OPS>
    std::vector<double> opsRootsFromJacobiMatrix(
        const OPS& ops, const unsigned degree)
    {
        std::vector<double> roots(degree);
        if (degree)
        {
            if (degree > ops.maxDegree()) throw std::invalid_argument(
                "In npstat::opsRootsFromJacobiMatrix: "
                "degree argument is out of range");
            const Matrix<double> jacobiMat = ops.jacobiMatrix(degree);
            assert(jacobiMat.nRows() == degree);
            assert(jacobiMat.nColumns() == degree);
            if (degree == 1U)
                roots[0] = jacobiMat[0][0];
            else
                jacobiMat.tdSymEigen(&roots[0], degree);
            const double s = ops.scale();
            const double shift = ops.location();
            if (!(shift == 0.0 && s == 1.0))
                for (unsigned i=0; i<degree; ++i)
                    roots[i] = roots[i]*s + shift;
            std::sort(roots.begin(), roots.end());
        }
        return roots;
    }
}
 
#endif // NPSTAT_OPSROOTSFROMJACOBIMATRIX_HH_