polyfilt.h

Go to the documentation of this file.

#ifndef POLYFILT_H
#define POLYFILT_H

#include "matkitdef.h"

#ifdef USE_NAMESPACE
using namespace MATKIT;
namespace MATKIT {
#endif

class Vector;
class Matrix;
class SymmetricMatrix;
class SparseMatrix;

#ifdef USE_NAMESPACE
}  // end of namespace MATKIT
#endif


//    Newton - Hermite Polynomial Interpolation


// build P(z) by Newton's divided differences in the form
//     P(z) = a(1) + a(2)*(z-x(1)) + a(3)*(z-x(1))*(z-x(2)) + ... + a(n)*(z-x(1))*...*(z-x(n-1)),
// such that P(x(i)) = y(i) for i=1,...,n, where
//     x,y are input vectors of length n, and a is the output vector of length n
// if x(i)==x(j) for some i!=j, then it is assumed that the derivative of P(z) is to be zero at x(i),
//     and the Hermite polynomial interpolation is applied
// in general, if there are k x(i)'s with the same value x0, then
//     the j-th order derivative of P(z) is zero at z=x0 for j=1,...,k-1
Vector NewtonPolynomial(const Vector &x, const Vector &y);


// return the evaluated P(z0), i.e. the value of P(z) at z=z0, where P(z) is a Newton polynomial defined by
//    P(z) = a(1) + a(2)*(z-x(1)) + a(3)*(z-x(1))*(z-x(2)) + ... + a(n)*(z-x(1))*...*(z-x(n-1))
// this routine also works for evluating the function value of a Hermite interpolating polynomial,
//    which is in the same form as a Newton polynomial
Real NewtonPolynomialEvaluation(const Vector &a, const Vector &x, const Real z0);
  // end of group_poly_intp



//    Base Filter

// begin of group_base_filter

// compute a base filter P(z) which is a continuous, piecewise polynomial P(z) expanded
// in a basis of `translated' (i.e. scale-and-shift) Chebyshev polynomials in each interval
//
// The base filter P(z) equals P_j(z) for z in the j-th interval [intv(j), intv(j+1)), where
// P_j(z) a Hermite interpolating polynomial
//
// input:
// intv is a vector which defines the intervals; the j-th interval is [intv(j), intv(j+1))
// HiLowFlags determines the shape of the base filter P(z)
// Consider the j-th interval [intv(j), intv(j+1)]
// HighLowFlag[j-1]==1,  P(z)==1 for z in [intv(j), intv(j+1)]
//                 ==0,  P(z)==0 for z in [intv(j), intv(j+1)]
//                 ==-1, [intv(j), intv(j+1)] is a transition interval;
//                       P(intv(j)) and P(intv(j+1)) are defined such that P(z) is continuous
// baseDeg is the degree of smoothness of the Hermite (piecewise polynomial) interpolation
// to be precise, the i-th derivative of P(z) is zero, i.e. d^{i}P(z)/dz^i==0, at all interval end points z=intv(j) for i=1,...,baseDeg
//
// output:
// P(z) expanded in a basis of `translated' (scale-and-shift) Chebyshev polynomials
// to be preicse, for z in the j-th interval [intv(j),intv(j+1)), P(z) equals
//     P_j(z) = pp(1,j)*S_0(z) + pp(2,j)*S_1(z) + ... + pp(n,j)*S_{n-1}(z),
// where S_i(z) is the `translated' Chebyshev polynomial in that interval,
//     S_i((z-c)/h) = T_i(z),  c = (intv(j)+intv(j+1))) / 2,  h = (intv(j+1)-intv(j)) / 2,
// with T_i(z) the Chebyshev polynomial of the first kind,
//     T_0(z) = 1, T_1(z) = z, and T_i(z) = 2*z*T_{i-1}(z) - T_{i-2}(z) for i>=2
// the return matrix is the matrix of Chebyshev coefficients pp just described
//
// note that the degree of P(z) in each interval is (at most) 2*baseDeg+1, with 2*baseDeg+2 coefficients
// let n be the length of intv; then there are n-1 intervals
// therefore the return matrix pp is of size (2*baseDeg+2)-by-(n-1)
Matrix HermiteBaseFilterInChebyshevBasis(const Vector &intv, const int *HiLowFlags, mkIndex baseDeg);

// an instance of this class, taken by GetIntervals() and PolynomialFilterInterface::setFilter(), is a collection of options to determine the intervals which decides the base filter
// it is used as an input argument of GetIntervals(), and also
struct IntervalOptions {

    // interval weights, default 100,1,1,1,100
    Vector intervalWeights;

    // the following parameter is for high-pass filters (and low-pass filters with conversion) only
    // the (relative) length of transition interval (default 0.6)
    Real transitionIntervalRatio;

    // the following parameters are for mid-pass filters only
    // this parameter tells whether to reverse the interval or not; it is effective only for mid-pass filters (default false)
    bool reverseInterval;
    // initial length of `plateau' relative to the length of the interval of desired eigenvalues (default 0.1)
    Real initialPlateau;
    // the rate at which the plateau shrinks at each iteration (default 1.5)
    Real plateauShrinkRate;
    // initial shift step, relative to the length of the interval of desired eigenvalues (default 0.01)
    Real initialShiftStep;
    // the rate at which the shift step expands (default 1.5)
    Real shiftStepExpansionRate;
    // maximum number of inner iterations (default 30)
    mkIndex maxInnerIter;
    // a mid-pass filter P(x) should have P(a1)=P(b1), where [a1,b1] is the interval of desired eigenvalues
    // yLimitTol (default 0.0001) is the tolerance allowed for |P(a1)-P(b1)|
    Real yLimitTol;

    // the following parameters are for all filters
    // maximum number of outer iterations, default 50
    mkIndex maxOuterIter;

    // number of grid points, used to measure the maximum function value of the polynomial P(z) for z not in the interval which the eigenvalues are sought (default 200)
    // the value of numGridPoints will automatically be increased if it is too small
    mkIndex numGridPoints;
    // this is the bottom line (default 0.001) which the function value of the polynomial must be greater than for z is in the interval of desired eigenvalues
    Real yBottomLine;
    // the limit of height of ripples (not in the interval of desired eigenvalues) relative to the bottom of polynomial values for z in the interval of desired eigenvalues (default 0.8)
    Real yRippleLimit;

    // a default constructor to set default parameters
    IntervalOptions() {
        // set default interval weights
        intervalWeights.resize(5);
        intervalWeights(1) = 100.0;
        intervalWeights(2) = 1.0;
        intervalWeights(3) = 1.0;
        intervalWeights(4) = 1.0;
        intervalWeights(5) = 100.0;

        // for high-pass filters (and low-pass filters with conversion)
        transitionIntervalRatio = 0.6;

        // for mid-pass filters
        reverseInterval = false;
        initialPlateau = 0.1;
        plateauShrinkRate = 1.5;
        initialShiftStep = 0.01;
        shiftStepExpansionRate = 1.5;
        maxOuterIter = 50;
        yLimitTol = 0.0001;

        // for all filters
        maxInnerIter = 30;
        numGridPoints = 200;
        yBottomLine = 0.001;
        yRippleLimit = 0.8;
    }
};

// the routine GetIntervals() returns an instance of this class, which gives the information of a polynomial filter
struct PolynomialFilterInfo {
    // the following returned values are for all filters
    // The partition of the range of spectrum which decides the base filter and the polynomial filter.
    Vector intervals;
    // 1 means a high-pass filter (or low-pass filter with conversion); 2 means a mid-pass filter
    int filterType;
    // 0 means no acceptable is found; 1 means an OK filter is found; 2 means an optimal filter is found
    int filterOK;
    // between 0.0 and 1.0; the higher the better quality of the filter
    Real filterQualityIndex;
    // number of iterations to get the (transition) intervals
    mkIndex numIter;

    // total number of iterations performed; this may include a few extra iterations without improving the (transition) intervals, compared to numIter
    mkIndex totalNumIter;
    // the lowest polynomial value P(z) for z in the interval [a1,b1] of desired eigenvalues
    Real yLimit;
    // the height of (highest, if more than one) summit in the interval [a1,b1] of desired eigenvalues
    Real ySummit;
    // number of steps moving leftward
    mkIndex numLeftSteps;
    // the height of highest summit in the left-hand side of the interval of desired eigenvalues
    Real yLeftSummit;
    // the height of lowest bottom in the left-hand side of the interval of desired eigenvalues
    Real yLeftBottom;

    // the following returned values are for high-pass filters (or low-pass filters after conversion)
    // in general it is |p(a1)-p(b1)|, where [a1,b1] is the interval of desired eigenvalues
    Real yLimitGap;
    // number of steps moving rightward
    mkIndex numRightSteps;
    // the height of highest summit in the right-hand side of the interval of desired eigenvalues
    Real yRightSummit;
    // the height of lowest bottom in the right-hand side of the interval of desired eigenvalues
    Real yRightBottom;

    // a default constructor to initialize all zero
    PolynomialFilterInfo() {
        // zero initialization
        filterOK = 0;
        filterQualityIndex = 0.0;

        numIter = 0;
        totalNumIter = 0;
        yLimit = 0.0;
        yLimitGap = 0.0;
        ySummit = 0.0;

        numLeftSteps = 0;
        yLeftSummit = 0.0;
        yLeftBottom = 0.0;

        numRightSteps = 0;
        yRightSummit = 0.0;
        yRightBottom = 0.0;
    }
};


// this routine determines the intervals (including the transition one(s)) by an interative process
//
// frame is a vector consisting of 4 ordered elements:
//     [frame(1),frame(4)] is the interval which (tightly) contains all eigenvalues, and
//     [frame(2),frame(3)] is the interval in which the eigenvalues are sought
// baseDeg is the left-and-right degree of the base filter for each interval
// polyDeg is the (maximum possible) degree of s(z), with z*s(z) the polynomial filter
// intv is the output; the j-th interval is [intv(j),intv(j+1))
// opts is a collection of interval options
//
// the base filter P(z) is a piecewise polynomial from Hermite interpolation with degree baseDeg at each end point of intervals
//
// the polynomial filter Q(z) is in the form z*s(z), i.e. Q(0)==0, such that ||(1-P(z))-z*s(z)||_w is minimized with
// s(z) a polynomial of degree up to polyDeg
//
// the resulting polynomial filter Q(z) satisfies Q(x)>=Q(y) for x in [frame(2),frame(3)], and
// y in [frame(1),frame(4)] but not in [frame(2),frame(3)]
//
// the routine returns an instance of PolynomialFilterInfo which gives some information of the polynomial filter
PolynomialFilterInfo GetIntervals(Vector &intv, const Vector &frame, mkIndex polyDeg, mkIndex baseDeg, IntervalOptions &opts);

// same as GetIntervals() above, with default IntervalOptions
PolynomialFilterInfo GetIntervals(Vector &intervals, const Vector &frame, mkIndex polyDeg, mkIndex baseDeg);  // default interval options will be used

// end of group_base_filter



//    Chebyshev Polynomials


// translate the coefficients of a Newton polynomial to the coefficients
// in a basis of the `translated' Chebyshev polynomials
//
// input:
// a Newton polynomial defined by vectors a and x:
//     P(z) = a(1) + a(2)*(z-x(1)) + a(3)*(z-x(1))*(z-x(2)) + ... + a(n)*(z-x(1))*...*(z-x(n-1))
// the interval [aa,bb] defines the `translated' Chebyshev polynomials S_i(z) = T_i((z-c)/h),
//     where c=(aa+bb)/2 and h=(bb-aa)/2, and T_i is the Chebyshev polynomial of the first kind
// note that T_i is defined by T_0(z)=1, T_1(z)=z, and T_i(z)=2*z*T_{i-1}(z)+T_{i-2}(z) for i>=2
//
// output:
// a vector q containing the Chebyshev coefficients, such that
//     P(z) = q(1)*S_0(z) + q(2)*S_1(z) + ... + q(n)*S_{n-1}(z)
Vector ExpandNewtonPolynomialInChebyshevBasis(Real aa, Real bb, const Vector &a, const Vector &x);


// evaluate P(z) at z=z0, where P(z) is a polynomial expanded in a basis of
// the `translated' (i.e. scale-and-shift) Chebyshev polynomials
//
// input:
// c is a vector of Chebyshev coefficients which defines the polynomial
//     P(z) = c(1)*S_0(z) + c(2)*S_1(z) + ... + c(n)*S_{n-1}(z),
// where S_i is the `translated' Chebyshev polynomial S_i((z-c)/h) = T_i(z), with
//     c = (intv(j)+intv(j+1)) / 2,  h = (intv(j+1)-intv(j)) / 2
// note that T_i(z) is the Chebyshev polynomial of the first kind,
//     T_0(z) = 1, T_1(z) = z, and T_i(z) = 2*z*T_{i-1}(z) - T_{i-2}(z) for i>=2
//
// output:
// the evaluated value of P(z) at z=z0
Real PolynomialEvaluationInChebyshevBasis(const Vector &c, Real z0, Real aa=-1.0, Real bb=1.0);


// evaluate P(z) at z=z0, where P(z) is a piecewise polynomial expanded
// in a basis of the (optionally translated, i.e. scale-and-shift) Chebyshev polynomials for each interval
//
// input:
// intv is a vector which defines the intervals; the j-th interval is [intv(j), intv(j+1))
// pp is a matrix of Chebyshev coefficients which defines a piecewise polynomial P(z)
// in a basis of the `translated' Chebyshev polynomials in each interval
// the polynomial P_j(z) in the j-th interval, i.e. when z is in [intv(j), intv(j+1)), is defined by the j-th column of pp:
//     if basisTranslated == false, then
//         P_j(z) = pp(1,j)*T_0(z) + pp(2,j)*T_1(z) + ... + pp(n,j)*T_{n-1}(z),
//     where T_i(z) is the Chebyshev polynomial of the first kind,
//         T_0(z) = 1, T_1(z) = z, and T_i(z) = 2*z*T_{i-1}(z) - T_{i-2}(z) for i>=2
//     if basisTranslated == true, then
//         P_j(z) = pp(1,j)*S_0(z) + pp(2,j)*S_1(z) + ... + pp(n,j)*S_{n-1}(z),
//     where S_i is the `translated' Chebyshev polynomial S_i((z-c)/h) = T_i(z), with
//         c = (intv(j)+intv(j+1)) / 2,  h = (intv(j+1)-intv(j)) / 2
//
// output:
// the evaluated value of P(z) at z=z0
//
// note that if z0 falls below the first interval, then the polynomial in the first interval will be used for evaluting P(z0)
//           if z0 flies over  the last  interval, then the polynomial in the last  interval will be used for evaluting P(z0)
Real PiecewisePolynomialEvaluationInChebyshevBasis(const Matrix &pp, const Vector &intv, Real z0, bool basisTranslated=true);


// compute the weighted inner product of two piecewise polynomials expanded
// in a basis of `translated' (i.e. scale-and-shift) Chebyshev polynomials for each interval
//
// pp and qq are two matrices of Chebyshev coefficients which define the piecewise polynomials P(z) and Q(z), respectively
// For z in the j-th interval, P(z) equals
//     P_j(z) = pp(1,j)*S_0(z) + pp(2,j)*S_1(z) + ... + pp(n,j)*S_{n-1}(z),
// and Q(z) equals
//     Q_j(z) = qq(1,j)*S_0(z) + qq(2,j)*S_1(z) + ... + qq(n,j)*S_{n-1}(z),
// where S_i(z) is the `translated' Chebyshev polynomial in that interval,
//     S_i((z-c)/h) = T_i(z),  c = (aa+bb)) / 2,  h = (bb-aa) / 2,
// with T_i(z) the Chebyshev polynomial of the first kind,
//     T_0(z) = 1, T_1(z) = z, and T_i(z) = 2*z*T_{i-1}(z) - T_{i-2}(z) for i>=2
//
// the (scaled) j-th interval inner product is defined by
//     <P_j,Q_j> = (Pi/2)*( pp(1,j)*qq(1,j) + sum_{k} pp(k,j)*qq(k,j) ),
// which comes from the property
//     <T_0,T_0>=pi, <T_i,T_i>=pi/2 for i>=1, and <T_i,T_j>=0 for i!=j
//
// the weighted inner product is <P,Q> = sum_{j} intervalWeights(j)*<P_j,Q_j>,
// which is the return value
//
// note that for unit weights, pass an empty vector of intervalWeights (i.e. of length 0)
Real PiecewisePolynomialInnerProductInChebyshevBasis(const Matrix &pp, const Matrix &qq, const Vector &intervalWeights);


// compute Q(z) = z*P(z), where P(z) and Q(z) are piecewise polynomials expanded
// in a basis of `translated' (i.e. scale-and-shift) Chebyshev polynomials for each interval
//
// P(z) and Q(z) are stored as matrices of Chebyshev coefficients pp and qq, respectively
//
// for z in the j-th interval, P(z) equals
//     P_j(z) = pp(1,j)*S_0(z) + pp(2,j)*S_1(z) + ... + pp(n,j)*S_{n-1}(z),
// and Q(z) equals
//     Q_j(z) = qq(1,j)*S_0(z) + qq(2,j)*S_1(z) + ... + qq(n,j)*S_{n-1}(z),
// where S_i(z) is the `translated' Chebyshev polynomial in that interval,
//     S_i((z-c)/h) = T_i(z),  c = (intv(j)+intv(j+1))) / 2,  h = (intv(j+1)-intv(j)) / 2,
// with T_i(z) the Chebyshev polynomial of the first kind,
//     T_0(z) = 1, T_1(z) = z, and T_i(z) = 2*z*T_{i-1}(z) - T_{i-2}(z) for i>=2
//
// the returned matrix is qq which represents Q(z) = z*P(z)
Matrix PiecewisePolynomialInChebyshevBasisMultiplyX(const Matrix &pp, const Vector &intv);
  // end of group_chebyshev



//    Conjugate Residual Method in the Polynomial Space


// this routine employs a conjugate-residual-type algorithm in polynomial space to minimize ||P(z)-Q(z)||_w,
// where P(z), the base filter, is the input piecewise polynomial, and
//       Q(z) is the output polynomial satisfying Q(0)==1, i.e. the constant term of Q(z) is 1
// niter is the number of conjugate-residual iterations; therefore, the degree of Q(z) is up to niter+1
// both P(z) and Q(z) are expanded in the `translated' (scale-and-shift) Chebyshev basis for each interval,
// and presented as matrices of Chebyshev coefficients, denoted by pp and qq, respectively
//
// input:
// intv is a vector which defines the intervals; the j-th interval is [intv(j),intv(j+1))
// w is a vector of Chebyshev weights; the weight of j-th interval is w(j)
//     the interval weights define the inner product of two continuous functions and then
//     the derived w-norm ||P(z)-Q(z)||_w
// pp is a matrix of Chebyshev coefficients which defines the piecewise polynomial P(z)
// to be specific, for z in [intv(j), intv(j+1)), P(z) equals
//     P_j(z) = pp(1,j)*S_0(z) + pp(2,j)*S_1(z) + ... + pp(niter+2,j)*S_{niter+1}(z),
// where S_i(z) is the `translated' Chebyshev polynomial in that interval,
//     S_i((z-c)/h) = T_i(z),  c = (intv(j)+intv(j+1))) / 2,  h = (intv(j+1)-intv(j)) / 2,
// with T_i(z) the Chebyshev polynomial of the first kind,
//     T_0(z) = 1, T_1(z) = z, and T_i(z) = 2*z*T_{i-1}(z) - T_{i-2}(z) for i>=2
//
// output:
// the return matrix, denoted by qq, represents a polynomial Q(z) with degree up to 1+niter
// and satisfying Q(0)==1, such that ||P(z))-Q(z)||_w is minimized
// this polynomial Q(z) is expanded in the `translated' Chebyshev basis for each interval
// to be precise, considering z in [intv(j), intv(j+1)), Q(z) equals
//     Q_j(z) = qq(1,j)*S_0(z) + qq(2,j)*S_1(z) + ... + qq(niter+2,j)*S_{niter+1}(z)
//
// note:
// 1. since Q(0)==1, P(0)==1 is expected; if P(0)!=1, one can translate P(z)
//    for example, if P(0)==0, one can use 1-P(z) as input instead of P(z)
// 2. typically, the base filter, defined by pp and intv, is from Hermite interpolation
//    in intervals [intv(j),intv(j+1)) for j=1,...,nintv, with nintv the number of intervals
Matrix FilteredConjugateResidualPolynomial(const Matrix &pp, const Vector &intv,
                                           const Vector &w, mkIndex niter);


// this routine employs a conjugate-residual-type algorithm in polynomial space to compute
// x = x0 + s(A)*r0 with r0 = b - A*x0, such that ||(1-P(z))-z*s(z)||_w is minimized, where
// P(z) is a given piecewise polynomial, called the base filter,
// s(z) is a polynomial of degree up to niter, the number of conjugate-residual iterations,
// and b and x0 are given vectors
//
// note that P(z) is expanded in the `translated' (scale-and-shift) Chebyshev basis for each interval,
// and presented as a matrix of Chebyshev coefficients pp
//
// input:
// A is a sparse matrix
// x0, b are vectors
// niter is the number of conjugate-residual iterations
// intv is a vector which defines the intervals; the j-th interval is [intv(j),intv(j+1))
// w is a vector of Chebyshev weights; the weight of j-th interval is w(j)
//     the interval weights define the inner product of two continuous functions and then
//     the derived w-norm ||P(z)-Q(z)||_w
// pp is a matrix of Chebyshev coefficients which defines the piecewise polynomial P(z)
// to be specific, for z in [intv(j), intv(j+1)), P(z) equals
//     P_j(z) = pp(1,j)*S_0(z) + pp(2,j)*S_1(z) + ... + pp(niter+2,j)*S_{niter+1}(z),
// where S_i(z) is the `translated' Chebyshev polynomial in that interval,
//     S_i((z-c)/h) = T_i(z),  c = (intv(j)+intv(j+1))) / 2,  h = (intv(j+1)-intv(j)) / 2,
// with T_i(z) the Chebyshev polynomial of the first kind,
//     T_0(z) = 1, T_1(z) = z, and T_i(z) = 2*z*T_{i-1}(z) - T_{i-2}(z) for i>=2
// tol is the tolerance; if the residual polynomial in z-norm is dropped by a factor lower
//     than tol, then stop the conjugate-residual iteration
//
// output:
// the return vector is x = x0 + s(A)*r0 with r0 = b - A*x0, such that ||(1-P(z))-z*s(z)||_w is minimized,
// subject to that s(z) is a polynomial of degree up to niter, where P(z) is the base filter
// in short, z*s(z) approximates 1-P(z)
//
// note:
// 1. since z*s(z) approximates 1-P(z), P(0)==1 is expected; if P(0)!=1, one can translate P(z)
//    for example, if P(0)==0, one can use 1-P(z) as input instead of P(z)
// 2. typically, the base filter, defined by pp and intv, is from Hermite interpolation
//    in intervals [intv(j),intv(j+1)) for j=1,...,nintv, with nintv the number of intervals
// 3. a common application is to compute R(A)*b, where R(z) approximates 1-P(z)
//    in this case, one can set x0 = 0 and then the return vector is x = s(A)*b, where
//    z*s(z) approximates 1-P(z); therefore, A*x is the wanted R(A)*b
// see also PolynomialFilterInterface::filteredSparseMatrixPolynomialVectorProduct()
Vector FilteredConjugateResidualMatrixPolynomialVectorProduct(const SparseMatrix &A, const Vector &x0, const Vector &b,
                                           const Matrix &pp, const Vector &intv, const Vector &w,
                                           mkIndex niter, Real tol=0.0);


// an interface to define R(S), i.e. set R and S, and compute R(S)*v, where
// S is a sparse symmetric matrix,
// R(z) is a polynomial which approximates a piecewise polynomial 1-P(z), with P(z) the base filter, and
// v is the input vector
//
// in short, use setFilter() to define R(S) once, and
// use filteredSparseMatrixPolynomialVectorProduct() to compute R(S)*v for every v
class PolynomialFilterInterface {
private:
    // works for SparseMatrix only
    // a pointer to the input sparse symmetric matrix via setFilter()
    static const SparseMatrix *Sptr;
    // the sparse symmetric matrix S in computing R(S)*v; S is translated from the input sparse symmetric matrix S0 in setFilter()
    static SparseMatrix S;

    // parameters for polynomial filtering

    // the j-th interval is [intervals(j),intervals(j+1))
    // the intervals are decided by GetIntervals() invoked by setFilter() which takes parameters frame, baseDeg, and polyDeg:
    // 1. frame is a vector consisting of 4 ordered elements:
    //     [frame(1),frame(4)] is the interval which (tightly) contains all eigenvalues of S, and
    //     [frame(2),frame(3)] is the interval in which the eigenvalues are sought
    // 2. baseDeg is the left-and-right degree of the base filter for each interval
    // 3. polyDeg is the (maximum possible) degree of z*s(z), with s(z) the polynomial filter
    static Vector intervals;
    // intervalWeights(j) is the weight of j-th interval [intervals(j), intervals(j+1))
    static Vector intervalWeights;
    // maximum possible degree of s(z), with z*s(z) the polynomial filter
    static mkIndex polyDegree;
    // left-and-right degree of the base filter in each interval
    static mkIndex baseDegree;

    // a base filter, which is a piecewise polynomial typically from Hermite interpolation
    // for j-th interval, the polynomial is expanded in the `translated' (scale-and-shift) Chebyshev basis,
    // with the Chebyshev coefficients stored in baseFilter(:,j)
    static Matrix baseFilter;
    // a zero vector of length the number of rows/columns of S
    static Vector zn;

public:
    // the information of computing the polynomial filter; see class PolynomialFilterInfo in polyfilt.h for more information
    static PolynomialFilterInfo filterInfo;

    // routines for the polynomial filtered Lanczos eigensolver

    // this member function sets the sparse matrix S translated from S0, and also the base filter P(z)
    // S0 is the sparse symmetric matrix of concern
    // frame is a vector of 4 ordered elements
    //     [frame(1),frame(4)] is the interval which (tightly) contains all eigenvalues of S0, and
    //     [frame(2),frame(3)] is the interval in which the eigenvalues of S0 are sought
    // polyDeg is the (maximum possible) degree of s(z), with z*s(z) the polynomial filter
    // baseDeg is the left-and-right degree of the base filter for each interval
    // opts is a collection of options to determine the intervals
    static Vector setFilter(const SparseMatrix &S0, const Vector &frame, mkIndex polyDeg, mkIndex baseDeg, IntervalOptions &opts);

    // this routine computes R(S)*v, where it is assumed that S and R(z) have been defined by setFilter()
    // R(z) is in the form z*s(z) which minimizes ||(1-P(z))-z*s(z)||_w by a conjugate-residual-type algorithm,
    // where P(z) is a piecewise polynomial as the base filter, and s(z) is a polynomial of degree up to polyDegree,
    // and w-norm is defined by data members intervals and intervalWeights
    //
    // note that P(z) is defined by the data members baseFilter and intervals
    // it is expanded in the `translated' (scale-and-shift) Chebyshev basis for each interval,
    // and presented as a matrix of Chebyshev coefficients baseFilter
    //
    // see also FilteredConjugateResidualMatrixPolynomialVectorProduct()
    static Vector filteredSparseMatrixPolynomialVectorProduct(const Vector &v) {
        return (*Sptr)*FilteredConjugateResidualMatrixPolynomialVectorProduct(*Sptr, zn, v, baseFilter, intervals, intervalWeights, polyDegree);
    }
};
  // end of group_cr_poly


#endif