Re: [eigen] Setting diagonal of a sparse matrix

[Brad Bell <bradbell@xxxxxxxxxx>]

I was using eigen-3.2.7. Changing to using hg to clone the current version of the repository from
     https://bitbucket.org/eigen/eigen/
the following example now runs (Thanks Peter):

# include <Eigen/SparseCore>
# include <iostream>

int main()
{   using Eigen::Dynamic;
    using Eigen::Lower;
    typedef Eigen::Matrix<double, Dynamic, 1> vector;
    typedef Eigen::SparseMatrix<double> sparse_matrix;
    typedef Eigen::TriangularView<sparse_matrix, Lower> sparse_lower_view;
    //
    sparse_matrix sA(3, 3);
    sA.setIdentity();
    sA.diagonal() = sA.diagonal() / 2.0;
    vector b      = vector::Ones(3);
    //
    sparse_lower_view sL(sA);
    vector   x   = sL.solve(b);
    std::cout << "x =\n" << x << "\n";
}


   
   
On 9/22/2016 7:24 AM, Brad Bell wrote:

My actual problem is quite complex. I do some computation, end of with a symmetric matrix, and must divide the diagonal by two. On top of that, in some cases, I must take the lower triangular view and do a solve. The exact mathematics can be found at
    http://www.coin-or.org/CppAD/Doc/cholesky.htm

I have a dense implementation in Eigen
    http://www.coin-or.org/CppAD/Doc/atomic_eigen_cholesky.hpp.htm
and verified that the corresponding derivative computations are correct
    http://www.coin-or.org/CppAD/Doc/atomic_eigen_cholesky.cpp.htm
I also have not yet figured out how to get a lower triangular view for sparse matrices (to do a solve on the lower triangular view). To be specific, the following code does not compile when I change the
    # if 0
to
    # if 1

# include <Eigen/SparseCore>
# include <iostream>

int main()
{   using Eigen::Dynamic;
    using Eigen::Lower;
    typedef Eigen::Matrix<double, Dynamic, Dynamic> matrix;
    typedef Eigen::Matrix<double, Dynamic, 1> vector;
    typedef Eigen::TriangularView<matrix, Lower> lower_view;
    typedef Eigen::SparseMatrix<double> sparse_matrix;
    typedef Eigen::TriangularView<sparse_matrix, Lower> sparse_lower_view;
    //
    matrix     A   = matrix::Identity(3, 3);
    vector     b   = vector::Ones(3);
    A.diagonal()     = A.diagonal() / 2.0;
    lower_view   L   = A.triangularView<Lower>();
    vector     x   = L.solve<Eigen::OnTheLeft>(b);
    std::cout << "x =\n" << x << "\n";
    //
    sparse_matrix sA(3, 3);
    for(size_t i = 0; i < 3; i++)
        sA.insert(i, i) = 0.5;
# if 0
    sparse_lower_view sL = sA.triangularView<Lower>();
    x   = sL.solve<Eigen::OnTheLeft>(b);
    std::cout << "x =\n" << x << "\n";
# endif
}

... snip ...