Re: [eigen] Skyline matrix

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


Hi,

thanks for the patch !

I'd be ok to have it in Eigen. In a first step I would probably put it in the unsupported/ set of modules because it sems quite preliminary for now.
However before we can apply the patch some cleaning is required: at the very least the copyrights have to be updated. It would also be nice to remove all the commented stuff.

Then we can discuss about the usefulness of having a SkylineMatrixBase. Indeed, such a base class can be needed for two reasons:

1 - To implement the _expression_ template paradigm for basics operations (+, -, scaling, coefficients wise operations like abs, sin, cos etc.). However, here I'm not sure that makes a lot of sense. It seems to me that a such a specialized matrix is only usefull for solving problems, so supporting such operators is not really needed, right ?

2 - Another reason would be to have various SkylineMatrix storage variants which would be sufficiently different to require different classes, but close enough to share the same (solver) algorithms. Again, I'm not sure that's the case here.

Nevertheless, it's still fine to keep it for now, and wait a bit for more use cases to see what's the best strategy.


gael.


On Wed, Nov 4, 2009 at 3:24 PM, guillaume saupin <guillaume.saupin@xxxxxx> wrote:
You'll find enclosed the code as a patch.


guillaume saupin wrote:
Hello !

I've finally (after some performances tests) implemented a basic skyline matrix using the following storage :

- an array stores the diagonal elements
- an array stores the upper elements in a column oriented way
- an array stores the lower elements in a row oriented way

I've implemented a LU decomposition for this format. I use the following algorithm (with minor modification to handle the sparseness of the matrix):

for each row
-- compute the lower part
  for each col < row
     LU(row, col) =  (LU(row, col)  - sum(k = 0 ; k < col) [LU(row, k)  * LU(k, col)] ) / LU(row, row) -- elements are access contiguously / Can be easily vectorize
-- compute the upper part
for each rrow < row
  LU(rrow, row) -= sum(k = 0 ; k < rrow) [LU(rrow, k)  * LU(k, row)] -- elements are access contiguously / Can be easily vectorize
-- compute the diag part
LU(row, row) -= sum(k = 0 ; k < row) [LU(rrow, k)  * LU(k, row)] -- elements are access contiguously / Can be easily vectorize

This algorithm is well suited for our storage format. After a few tests, it seems even more efficient to store the the upper and lower elements in the same array, with the elements from one row followed by the elements of the cols of same index, and so on ...

I've also implemented a product for the Skyline *  Dense operation..

Do you think that it would be possible to insert this code into eigen ?

guillaume saupin wrote:
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 is more like an extension of the tridiagonal storage.

I think that our storage is actually quite similar to the one proposed there for non symetric skyline (referred as NSK format by Saad) matrix. As in these format, we store the upper element in a column-oriented  way, and the lower element in a row-oriented way.
The only difference is that we store the diagonal elements in a separate array instead of keeping an index array to these elements. This allows contiguous access to diagonal elements.


However, I don't see how your storage can yield to a more efficient implementation of the LU decomposition than the standard skyline storage. For instance here is a basic algorithm without pivoting for the 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 2
 lu.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) since lu.col(k)..end(rrows) is just a small dense vector.

Line 2 is an outer product which again is trivially/automatically vectorized 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 sequentially stored in memory.

But perhaps there exists a special algorithm which perfectly fit your storage ? Is it possible to see your code somewhere ? Finally, if it is really more efficient then that would make sense to have it in Eigen.

The benefit of using the NSK, i.e. of storing the lower matrix in a column-oriented way, and the upper matrix in a row-oriented way is that in line 1 we access contigous elements (the column k of the lower matrix), and in line 2 we also access contiguous element (the row elements of the upper matrix). I don't know if it's efficient for vectorization, but without, it is quite efficient, as the locality is very good.



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


















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