Re: [eigen] Skyline matrix |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Skyline matrix*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Sat, 14 Nov 2009 21:06:29 -0500*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=BYtKFWbW0EA/BtSmq1eoEit01ON3Qrkf6opfQnXVlm8=; b=FvK9mIQ5F3PrZsEoiasime90/KiWnB2XaIJruxXnVFORESFTzjTRRb21mML3LpRIUb qGasyAyJ9tqDpPHM8RS/FkPJ+gAguA7zN+rkZQ8R34vbDflbDogA0GLDEPBF8oLiMYYa hR0pWFz8VosCyiVWlR8Zt7WLv9EiEr9eyH9zU=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=u6KG5q/hVFuqBYeGdvg/YUTXt9q8nVi62/3gQd9D+0674RG5TmtiZ3ZZtzOg6szRDB xbiHox9esX2FWzugcLAMG6RC66J2s4vB10GsXfmHjOd5RpsiDHzeqCyDXYHEd9nlLTpt q3oCKI9q9cA/qA6eniJe1SVbUUmLefK7DSi9M=

Hey! kayhman contact intra cea fr, it's you? I just realized your email address was appearing unobfuscated in the credits. The script was assuming FirstName LastName <address> which is a reasonable assumption... anyway, if we don't have your name, how can we credit you? So please edit your config as shown, http://eigen.tuxfamily.org/index.php?title=Developer%27s_Corner#Configuring_Mercurial Benoit 2009/11/4 guillaume saupin <guillaume.saupin@xxxxxx>: > 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 >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>> >>> >>> >> >> >> >> > >

**Follow-Ups**:**Re: [eigen] Skyline matrix***From:*guillaume saupin

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

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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**[eigen] Segfault in test matrixExponential_2** - Next by Date:
**Re: [eigen] Segfault in test matrixExponential_2** - Previous by thread:
**Re: [eigen] Skyline matrix** - Next by thread:
**Re: [eigen] Skyline matrix**

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