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 <[email protected]> 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 <[email protected]> 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. >>>> >>> > >
