| Re: [eigen] Skyline matrix |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
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