public inbox for gsl-discuss@sourceware.org
 help / color / mirror / Atom feed
From: Patrick Alken <alken@colorado.edu>
To: Alexis Tantet <alexis.tantet@gmail.com>
Cc: "gsl-discuss@sourceware.org" <gsl-discuss@sourceware.org>
Subject: Re: Sparse matrix extension
Date: Sun, 07 Feb 2016 00:03:00 -0000	[thread overview]
Message-ID: <56B689B1.5090005@colorado.edu> (raw)
In-Reply-To: <CAMWWPT2d=PnK1cexcZ+OLQAYJudRk9QLGac-Wcn33OaES5UjJA@mail.gmail.com>

Hello,

  I've been working my way through your patch. There is a new branch on
the git called 'sparse' with all the current changes. You've done a lot
of work on your patch so its taking me a while to get through all of it.

  The main functionality of CRS is in there, along with basic stuff like
memcpy, transpose_memcpy, transpose. Also dgemv works (and hence also
the GMRES linear solver). I did some initial work on the fprintf/fscanf
for the triplet case - for the compressed matrices, I am wondering if we
should stick to a standard format such as Harwell-Boeing, I'm not sure
at the moment. We might be able to avoid this issue by making a binary
fwrite/fread interface which just directly writes the 3 arrays to disk.

Things I still need to add for the CRS storage:

1. spmatrix_add
2. spblas_dgemm
3. fprintf/fscanf
4. transpose - for now I'll keep it the way you coded it, but I think it
would be better to keep the matrix in the same storage format during a
transpose, ie: the user might want to pass A^T to an external solver
library which requires CCS and they might not expect the transpose
routine to change it to CRS. I did a quick search for in-place transpose
algorithms in CCS/CRS and they do appear to exist though I haven't had
time to implement it yet.

I removed the inner/outer idx stuff from gsl_spmatrix - its not
necessary and also I'd like to avoid modifying the gsl_spmatrix struct
since any changes would break binary compatibility for the next release.
Also your modification to gsl_spmatrix_set - I see what you're trying to
do but I think its better to add a new function:

double *gsl_spmatrix_ptr(gsl_spmatrix *m, size_t i, size_t j)

which returns a pointer to the data array for element (i,j), and the
user can then add whatever value they want (this is how its done with
gsl_matrix and gsl_vector). I think it should be straightforward to add
this function, possibly even for the compressed storage too.

Anyway, just wanted to let you know things are progressing (though a
little slowly).

Patrick

On 01/20/2016 11:31 AM, Alexis Tantet wrote:
> Hi Patrick,
>
> I've cloned gsl.git, made my modifications to spmatrix/ and spblas/
> and pushed them to master at:
> https://github.com/atantet/gsl
>
> This time, I've tried to modify the code as little as possible. I've
> usually added comments labelled "AT" with each modification.
> I also have a few questions about the code which I've labelled "?AT".
> Finally, I've added the corresponding tests to both test.c (successful
> on my machine).
>
> I see that duplicates are handled for the triplet format by replacing.
> Thus, I removed the function taking care of sorting and summing
> duplicates that I had added to _compress (the trick was to transpose
> with copy to reorder first, transpose back and remove the duplicates
> then, taken from Eigen C++). Instead, I have added an option to sum
> the duplicates (useful e.g. when counting links of a graph or
> transitions of a Markov chain). One may need the third option of
> repeating the triplets, but then avl_insert should also be modified.
>
> However, if I've understood well, no sorting of the inner indices is
> performed right now. This could be problematic in the future (e.g. in
> functions such as _equal).
>
> For _spdgemm for CRS matrices, I've used the trick (A B)^t = B^t A^t
> to convert the problem to CSC. I've tested it successfully but you may
> not like the coding style.
>
> I hope that helps,
> Best,
> Alexis
>
> On Tue, Jan 19, 2016 at 9:50 PM, Patrick Alken <alken@colorado.edu> wrote:
>>
>> On 01/19/2016 12:55 PM, Alexis Tantet wrote:
>>> Hi Patrick,
>>>
>>> Thanks for the quick reply! I wanted to be sure that this contribution
>>> is useful before to spend time on the merging with the latest version.
>>> I will create the gsl.git repository and work on it during the week.
>>>
>>> I had already had a look at the documentation but did not know about
>>> the iterative solvers (a link between each modules would be useful).
>>> My contribution indeed fits in the sparse matrix module + the update
>>> of the dgemm and dgemv functions to support CRS (an update may also be
>>> needed for the solvers).
>>
>> The solvers only call dgemv (as far as I remember) so they shouldn't need an
>> update once dgemv is updated.
>>
>>> I have also developed a simple C++ object allowing to use gsl_spmatrix
>>> as a user-defined matrix in ARPACK++ (a maintained fork of ARPACK++
>>> can be found at https://github.com/m-reuter/arpackpp), allowing to
>>> avoid having to use other libraries such as superLU. It could be
>>> useful to others, maybe as an extension. Now that I think about it,
>>> the iterative solvers could also be used to support the shift and
>>> invert modes (see ARPACK++ documentation). What do you think (I could
>>> work on it)?
>> I've never used ARPACK, but if you want to make an extension to interface
>> GSL/ARPACK its fine with me - such an extension would never be added to the
>> main GSL code since GSL tries to be as standalone as possible.
>>
>>
>>> If you have major comments, the sooner the better, so that I can work
>>> on them while merging.
>>>
>>> Thank you for your interest,
>>> Alexis
>>>
>>> On Tue, Jan 19, 2016 at 6:02 PM, Patrick Alken <alken@colorado.edu> wrote:
>>>> Hi Alexis,
>>>>
>>>>    This looks like very good work! Adding compressed row storage has been
>>>> on
>>>> my todo list for a while. The 'gslsp' extension is unfortunately very out
>>>> of
>>>> date, and the current git contains newer code (including a GMRES
>>>> iterative
>>>> linear solver). I removed the gslsp extension from the web page a while
>>>> back
>>>> to reflect this. You can browse the latest manual to see the current
>>>> sparse
>>>> matrix capabilities (http://www.gnu.org/software/gsl/manual/gsl-ref.pdf)
>>>> -
>>>> there are 3 chapters: sparse matrices, sparse blas and sparse linear
>>>> algebra
>>>> - it looks like your contributions will fit into the sparse matrices
>>>> chapter.
>>>>
>>>>    Would you be able to verify that your changes are compatible with the
>>>> current gsl.git repository? This will make it much easier for me to merge
>>>> everything into the git when ready. It would be best if you made a new
>>>> branch of gsl.git, and add your changes so I can then pull them from
>>>> github
>>>> or somewhere. I will try to find some time in the next few days to look
>>>> over
>>>> your code.
>>>>
>>>> Thanks again,
>>>> Patrick
>>>>
>>>>
>>>> On 01/19/2016 09:43 AM, Alexis Tantet wrote:
>>>>> Dear GSLers,
>>>>>
>>>>> As a scientist rather than a developer, I have developed an extension
>>>>> of the sparse matrix module (CRS, I/O, manipulation, see below), which
>>>>> I have tested. These modifications conserve the structure of the
>>>>> original module and be useful for a large number of sparse matrices
>>>>> users.
>>>>>
>>>>> I'm not familiar with the contributing process here. My repository can
>>>>> be found there:
>>>>> https://github.com/atantet/gslsp
>>>>> Unfortunately, I did not know of the gsl.git repository and I forked it
>>>>> froml:
>>>>> https://github.com/drjerry/gslsp ,
>>>>> which seems to be a bit older than gsl.git.
>>>>>
>>>>> How can I push/merge to gsl.git ? Should it be as an update or another
>>>>> extension? Is it necessary to adapt to the newest version of the code
>>>>> ?
>>>>>
>>>>> Best regards,
>>>>> Alexis Tantet
>>>>>
>>>>> CHANGES.md:
>>>>>
>>>>> Extension of the sparse matrix module of GSL
>>>>>
>>>>> ===================================
>>>>>
>>>>> Introduction
>>>>> ------------
>>>>>
>>>>> Usages of sparse matrices are numerous in scientific computing.
>>>>> When numerical linear algebra problems become large, sparse
>>>>> matrices become necessary to avoid memory overload and unnecessary
>>>>> computations, at the cost of element access and matrix construction.
>>>>> As a result, most large scale linear solvers or eigen solvers perform
>>>>> on sparse matrices.
>>>>>
>>>>> Fortunately, a very useful sparse matrix module has recently been
>>>>> introduced to GSL.
>>>>> However, important features are still lacking, such has
>>>>> Compressed Row Storage (CRS) matrices, input/output functions and
>>>>> other matrix properties and manipulation functions.
>>>>> This new version attempts to address this, conserving the original
>>>>> structure of the module and conventions.
>>>>>
>>>>> Major changes
>>>>> -------------
>>>>>
>>>>> * Add CRS format and update functions manipulating compressed matrices :
>>>>>        - additional flag GSL_SPMATRIX_CRS and macro GSLSP_ISMATRIX (
>>>>> gsl_spmatrix.h )
>>>>>        - additional members innerSize and outerSize used to iterate
>>>>> matrix elements ( gsl_spmatrix.h )
>>>>>        - rename some variables for coherence ( gsl_spmatrix.h , *.c )
>>>>>        - update all functions on compressed matrices ( *.c )
>>>>> * Allow to sum duplicate elements when compressing ( spcompress.c ) :
>>>>>        - modify gsl_spmatrix_compress
>>>>>        - add    gsl_spmatrix_sum_duplicate
>>>>> * CCS <-> CRS and fast transpose inplace in spswap.c :
>>>>>        - add gsl_spmatrix_switch_major
>>>>>        - add gsl_spmatrix_transpose
>>>>> * Add printing and scanning functions in spio.c :
>>>>>        - add gsl_spmatrix_fprintf
>>>>>        - add gsl_spmatrix_fscanf
>>>>> * Add manipulation functions in spmanip.c (particularly useful for
>>>>> Markov chain transition matrices) :
>>>>>        - add gsl_spmatrix_get_rowsum : get vector of sum over row
>>>>> elements
>>>>>        - add gsl_spmatrix_get_colsum : get vector of sum over column
>>>>> elements
>>>>>        - add gsl_spmatrix_div_rows   : divide all elements of each row
>>>>> by a vector element
>>>>>        - add gsl_spmatrix_div_cols   : divide all elements of each
>>>>> column by a vector element
>>>>> * Add test functions in atprop.c :
>>>>>        - add gsl_spmatrix_gt_elements : greater than test for each
>>>>> matrix
>>>>> element
>>>>>        - add gsl_spmatrix_ge_elements : greater or equal than test for
>>>>> each matrix element
>>>>>        - add gsl_spmatrix_lt_elements : lower than test for each matrix
>>>>> element
>>>>>        - add gsl_spmatrix_le_elements : lower or equal than test for
>>>>> each matrix element
>>>>>        - add gsl_spmatrix_any         : test if any non-zero element in
>>>>> matrix
>>>>>
>>>>> Other minor changes have been made, such as error tests.
>>>>> test.c has also been updated to test new features.
>>>>
>>>
>
>

  reply	other threads:[~2016-02-07  0:03 UTC|newest]

Thread overview: 25+ messages / expand[flat|nested]  mbox.gz  Atom feed  top
     [not found] <CAMWWPT3uJj4Vrn7ut6+F18gY===zd6+1r1UJhz0hcCj--zwtdg@mail.gmail.com>
2016-01-19 16:43 ` Alexis Tantet
2016-01-19 17:02   ` Patrick Alken
2016-01-19 19:55     ` Alexis Tantet
2016-01-19 20:50       ` Patrick Alken
2016-01-20 18:31         ` Alexis Tantet
2016-02-07  0:03           ` Patrick Alken [this message]
2016-02-07  1:03             ` Alexis Tantet
2016-02-07 17:25               ` Patrick Alken
2016-02-07 19:32                 ` Alexis Tantet
2016-02-07 20:14                   ` Patrick Alken
2016-02-07 20:31                     ` Alexis Tantet
2016-02-07 21:34                       ` Patrick Alken
2016-02-08  0:59                         ` Alexis Tantet
2016-02-10 13:16                           ` Alexis Tantet
2016-02-10 13:48                             ` Alexis Tantet
2016-02-10 15:56                             ` Patrick Alken
2016-02-12 10:43                               ` Alexis Tantet
2016-02-13 19:42                                 ` Patrick Alken
2016-02-14 11:06                                   ` Alexis Tantet
2016-02-14 18:11                                     ` Patrick Alken
2016-02-14 18:25                                       ` Brian Gladman
2016-02-15  5:10                                         ` Patrick Alken
2016-02-15 11:09                                           ` Brian Gladman
     [not found]                                             ` <CAMWWPT0J9ENRZjJHLO=cxot4DGdSLer+n2HkBVnhFnO0oiVV8g@mail.gmail.com>
2016-02-15 13:56                                               ` Alexis Tantet
2016-02-15 17:17                                                 ` Patrick Alken

Reply instructions:

You may reply publicly to this message via plain-text email
using any one of the following methods:

* Save the following mbox file, import it into your mail client,
  and reply-to-all from there: mbox

  Avoid top-posting and favor interleaved quoting:
  https://en.wikipedia.org/wiki/Posting_style#Interleaved_style

* Reply using the --to, --cc, and --in-reply-to
  switches of git-send-email(1):

  git send-email \
    --in-reply-to=56B689B1.5090005@colorado.edu \
    --to=alken@colorado.edu \
    --cc=alexis.tantet@gmail.com \
    --cc=gsl-discuss@sourceware.org \
    /path/to/YOUR_REPLY

  https://kernel.org/pub/software/scm/git/docs/git-send-email.html

* If your mail client supports setting the In-Reply-To header
  via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).