Hi guys,

I need a linear solver (Ax=b). Someone proposed Eigen for ITK (itk.org), and since my project is written in C++, it is a logical choice for me too.

My matrix is somewhat sparse and has dimensions few hundred squared. I will very rarely change the matrix A, but vector b will be different in every iteration. I want to accomplish this with minimum amount of coding (and subsequent maintenance), and I am willing to sacrifice some performance for it.

I have a few questions though:
1. Which version should I use, latest stable release (2.0.15), V3 beta2 (expected soon), or the development version (from Hg repo)?
2. How long does it take (approximately) to solve such a matrix (lets say 300x300)? Is it milliseconds or seconds?
3. Since my matrices (a few of them) are static, I can even precalculate some stuff before compilation and store it in a file or such. Can I get any speed advantage out of that? And what is the way?

I am just starting (to read tutorials and docs), and for that I need to settle on a version I am going to use. Forgive my ignorance if expressed in this post

Regards,
Dženan

Somewhat OT:
ITK is coupled with VTK, isn't it?
Benoit, are there some news about "merriage" of VTK and Eigen?
--Frank

VTK is made by the same company which maintains ITK. However, I currently use OpenSceneGraph, not VTK (so they are not "married" ).

Any suggestions about the version at least?

I just realized that the answer to my 3rd question is to calculate the inverse matrix, and then use that repeatedly.

If someone cares: I decided to use the latest (dev) version.

dzenanz wrote:I just realized that the answer to my 3rd question is to calculate the inverse matrix, and then use that repeatedly.

Hi,

As stated somewhere in the docs, you may check first that obtaining a decomposition (LDLT, LU, SVD, etc) then repeatedly apply solve() or solveInPlace() is not faster and more accurate than obtaining the inverse of your matrix. I'm quoting the doc by memory: "Inverses are popular in the field of mathematics, but less in numerical analysis." The doc also states that the above is not true for small matrices. I have the feeling that 300x300 is not that small, so there is a place for experimentation.

Cheers,
Romain

I recommend you to start with Eigen 3, and if your matrices are not very very sparse, then for such small problems (yes 300x300 is very small when we talk about sparse problems), I recommend you to use dense decompositions. For instance, to solve:

A X = B

with both A, B and X 300x300 of double, it takes ~0.017sec with LU, and 0.012sec if A is selfadjoint and using Cholesky. Here is my bench example:

#include <iostream>
#include <Eigen/LU>
#include <Eigen/Cholesky>
#include <bench/BenchTimer.h>
using namespace Eigen;

int main()
{

MatrixXd mat(300,300), x(300,300);
mat.setRandom();

BenchTimer t;
BENCH(t,4,10,x = mat.lu().solve(x));
std::cout << "LU : " << t.best()/10 << "\n";

t.reset();
BENCH(t,4,10,x = mat.selfadjointView<Lower>().llt().solve(x));
std::cout << "LLT : " << t.best()/10 << "\n";


return 0;
}

Thanks for stability recommendation (I might try it) and for times.