Re: [eigen] Setting diagonal of a sparse matrix

[Peter <list@xxxxxxxxxxxxxxxxx>]

Dear Brad,



Am 22.09.2016 um 16:24 schrieb Brad Bell:
> 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


Here's a version of your code that compiles and runs.
You may want to use a different solver.


Best regards,
Peter

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

int main()
{
	typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix;
	typedef Eigen::Matrix<double, Eigen::Dynamic, 1> vector;
	typedef Eigen::TriangularView<matrix, Eigen::Lower> lower_view;
	typedef Eigen::SparseMatrix<double> sparse_matrix;
	typedef Eigen::TriangularView<sparse_matrix, Eigen::Lower> sparse_lower_view;

	vector     b   = vector::Ones(3);

	matrix     A   = matrix::Identity(3, 3);

	A.diagonal()     = A.diagonal() / 2.0;
	lower_view   L   = A.triangularView<Eigen::Lower>();
	vector     x   = L.solve<Eigen::OnTheLeft>(b);
	std::cout << "x =\n" << x << "\n";
	//
	sparse_matrix sA(3, 3);
	sA.setIdentity();

	sA.diagonal() = Eigen::Matrix<double, 3, 1>::Constant(0.5);

	std::cout << "sA: " << sA << std::endl;

	sparse_lower_view sL(sA);
	Eigen::ConjugateGradient<sparse_matrix , Eigen::Lower> solver;

	x   = solver.compute(sA).solve(b);
	std::cout << "x =\n" << x << "\n";

}