Re: [eigen] Skyline matrix |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Skyline matrix*From*: guillaume saupin <guillaume.saupin@xxxxxx>*Date*: Wed, 28 Oct 2009 09:09:41 +0100

Gael Guennebaud a écrit :

On Tue, Oct 27, 2009 at 10:40 AM, guillaume saupin<guillaume.saupin@xxxxxx <mailto:guillaume.saupin@xxxxxx>> wrote:Gael Guennebaud a écrit : this is not an easy question because we are still unsure how to manage efficiently and with minimal code all these kind of "special" matrices. Of course an easy solution would be to make it inherit SparseMatrixBase and let the generic algorithms based on sparse iterators do the job. However, this approach won't be optimal because the iterators won't take advantage of the specificity of a SkylineMatrix which, e.g., should allow vectorization. Actually skyline matrices are very similar to banded matrices, and in particular it seems to me that they are similar enough to share the same algorithms. So for instance these two kind of matrices could inherit the same matrix base and the algorithms could be implemented in a generic way via the concept of "range inner vector" which is a column (or a row) vector with an associated start index... We also need efficient basic operators dealing with such vector. E.g. the addition of two such vector can be efficiently done by first copying the non overlapping part (start and end), and then sum the overlapping part using Eigen's core module feature. The way we implemented skyline matrix is not really similar to what is usually done for banded matrix, as we store the diagonal, the upper, and the lower elements in separate arrays. This way, we can achieve quite efficient LU decomposition (and solving) due to coherent memory accesses.ok I see, then it is not really what people call skyline storage(http://www.netlib.org/linalg/html_templates/node96.html), but it ismore like an extension of the tridiagonal storage.

However, I don't see how your storage can yield to a more efficientimplementation of the LU decomposition than the standard skylinestorage. For instance here is a basic algorithm without pivoting forthe standard skyline storage:for(int k = 0; k+1 < rows; ++k) { int rrows = rows-k-1; int rsize = size-k-1; // line 1 lu.col(k).end(rrows) /= lu.coeff(k,k); // line 2lu.corner(BottomRight,rrows,rsize).noalias() -= lu.col(k).end(rrows)* lu.row(k).end(rsize);}where lu is the working matrix, lu.col(k) is assumed to return a"range vector" as I described in my previous email.Here line 1 would be trivially optimized (i.e., vectorized) sincelu.col(k).end(rrows) is just a small dense vector.Line 2 is an outer product which again is trivially/automaticallyvectorized sicne it is impelmented as a sequence of:"col_range_vector_i -= scalar * col_range_vector".Here the locality is pretty good because the vector"lu.col(k).end(rrows)" which is reused multiple times is sequentiallystored in memory.But perhaps there exists a special algorithm which perfectly fit yourstorage ? Is it possible to see your code somewhere ? Finally, if itis really more efficient then that would make sense to have it in Eigen.

gael. So I don't think that our algorithms share common features with band matrix algorithms. They seem closer to what can be done with diagonal matrices. Anyway, I've used the SparseMatrix code as a base for my SkylineMatrix code. Currently it only supports very basic operations. I'll implement the LU decomposition soon. I will also try to use vectorization. guillaume These are just some initial thoughts and the discussion is very open! gael On Mon, Oct 19, 2009 at 12:55 PM, guillaume saupin <guillaume.saupin@xxxxxx <mailto:guillaume.saupin@xxxxxx> <mailto:guillaume.saupin@xxxxxx <mailto:guillaume.saupin@xxxxxx>>> wrote: Hello, We are planning to use your library in our projects, but we need a skyline matrix. Therefore I'd like to implement one for eigen, but I don't now where to start. Is there a specific class that can be a good starting point / skeleton to start with ? The SparseMatrixBase might be a good choice. Should this SkylineMatrix inherit from SparseMatrixBase, or be a separate class ? Thanks, guillaume

**References**:**[eigen] Skyline matrix***From:*guillaume saupin

**Re: [eigen] Skyline matrix***From:*Gael Guennebaud

**Re: [eigen] Skyline matrix***From:*guillaume saupin

**Re: [eigen] Skyline matrix***From:*Gael Guennebaud

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] rvalue refs - std::swap** - Next by Date:
**Re: [eigen] Mapping array of scalars into quaternions** - Previous by thread:
**Re: [eigen] Skyline matrix** - Next by thread:
**[eigen] Eigen is slow with complex matrices**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |