#include <flens/flens.h>

namespace flens {


struct MySyMatrix
    : public SymmetricMatrix<MySyMatrix>
{
    typedef double ElementType;
    typedef int    IndexType;
};

} // namespace flens


namespace flens { namespace blas {

template <typename ALPHA, typename VX, typename BETA, typename VY>
void
mv(const ALPHA &alpha, const MySyMatrix &A, const DenseVector<VX> &x,
   const BETA &beta, DenseVector<VY> &y)
{
    typedef typename DenseVector<VX>::IndexType  IndexType;
    const IndexType n = x.length();

    // we only resize empty left hand sides
    if (y.length()==0) {
        y.resize(n);
    }
    ASSERT(y.length()==n);


    for (IndexType i=1; i<=n; ++i) {
        y(i) = beta*y(i) + alpha*2*x(i);
        if (i>1) {
            y(i) -= alpha*x(i-1);
        }
        if (i<n) {
            y(i) -= alpha*x(i+1);
        }
    }
}

} } // namespace blas, flens