Re: [eigen] Setting diagonal of a sparse matrix

[Brad Bell <bradbell@xxxxxxxxxx>]

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
}





On 9/22/2016 7:05 AM, Christoph Hertzberg wrote:
> I'm not sure if there is a better way to 'replace' the diagonal, if it 
> does exist partially already. But setting the entire Matrix S to the 
> diagonal of D goes like this:
>   S = D.diagonal().asDiagonal();  // clears all off-diagonals, of course
> If you know that there are no entries yet (or don't care), you can do 
> the same with += or -= instead.
> As the documentation says, S.diagonal()=... only works if all diagonal 
> entries already exist.
>
>
> Christoph
>
> On 2016-09-22 15:14, Brad Bell wrote:
> > How should one efficiently set the diagonal of an Eigen sparse matrix ?
> >
> > The following program does not compile if I set
> >     USE_SPARSE_DIAGONAL 1
> >
> > # include <Eigen/SparseCore>
> > # include <iostream>
> >
> > # define USE_SPARSE_DIAGONAL 0
> > int main()
> > {   using Eigen::Dynamic;
> >     Eigen::Matrix<double, 3, 3> D = Eigen::Matrix<double, 3,
> > 3>::Identity();
> >     Eigen::SparseMatrix<double> S(3, 3);
> >     D.diagonal() = D.diagonal() / 2.0;
> > # if USE_SPARSE_DIAGONAL
> >     S.diagonal() = D.diagonal();
> > # else
> >     for(size_t i = 0; i < 3; i++)
> >         S.coeffRef(i, i) = D(i, i);
> > # endif
> >     std::cout << "S =\n" << S << "\n";
> > }
> > ~