From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 40224 invoked by alias); 7 Feb 2016 00:03:11 -0000 Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org Received: (qmail 38691 invoked by uid 89); 7 Feb 2016 00:03:09 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-0.0 required=5.0 tests=BAYES_20,RCVD_IN_DNSWL_NONE,SPF_HELO_PASS,SPF_PASS autolearn=ham version=3.3.2 spammy=Storage, coded, manipulation, HX-Microsoft-Antispam-PRVS:sk:BY1PR03 X-HELO: na01-by2-obe.outbound.protection.outlook.com Received: from mail-by2on0125.outbound.protection.outlook.com (HELO na01-by2-obe.outbound.protection.outlook.com) (207.46.100.125) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with (AES256-SHA256 encrypted) ESMTPS; Sun, 07 Feb 2016 00:03:06 +0000 Authentication-Results: sourceware.org; dkim=none (message not signed) header.d=none;sourceware.org; dmarc=none action=none header.from=colorado.edu; Received: from [192.168.0.75] (75.166.52.243) by BY1PR0301MB0903.namprd03.prod.outlook.com (10.160.195.142) with Microsoft SMTP Server (TLS) id 15.1.403.16; Sun, 7 Feb 2016 00:03:03 +0000 Subject: Re: Sparse matrix extension To: Alexis Tantet References: <569E6C33.1090505@colorado.edu> <569EA1A9.2080101@colorado.edu> CC: "gsl-discuss@sourceware.org" From: Patrick Alken X-Enigmail-Draft-Status: N1110 Message-ID: <56B689B1.5090005@colorado.edu> Date: Sun, 07 Feb 2016 00:03:00 -0000 User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:38.0) Gecko/20100101 Thunderbird/38.4.0 MIME-Version: 1.0 In-Reply-To: Content-Type: text/plain; charset="utf-8" Content-Transfer-Encoding: 7bit X-ClientProxiedBy: DM2PR10CA0060.namprd10.prod.outlook.com (10.141.241.28) To BY1PR0301MB0903.namprd03.prod.outlook.com (25.160.195.142) X-MS-Office365-Filtering-Correlation-Id: 81b25709-e0c1-4bb1-7025-08d32f520cf9 X-Microsoft-Exchange-Diagnostics: 1;BY1PR0301MB0903;2:2wUerVfxUuyIS58v5vb2p7kzXBy0tqK9e0RtE8+furPszwaI/JE6v+axHur9Bmdsqkb78Xp6oRS9q9TdKlDwTOHHsjSBSsoU6UuV2UAE2vybFpTjrMiAPKBVLX5slm0CVJ78uhvWrDOhxrBZUPKl9lUJ7fyKmawzuHMAt1FJ4bCWYvxkDbuyMJ5ybpQw09Xn;3:fCKy+H1QAPJmfVpWDrW+fIz0H6hd20qXe4D43kPTF77YEtwQyMbEHquXZfhhSks07UGwnOuqIR7qBTgEp3CO/3MrcCK1gnNZ78/nfA3U1IUfpGQqJygq0xgGmy7LLi/2;25:+HyYXNh0ntAAHOVAOfJROGEByY7GhmM8ODudq/4cn3TdZTUEWtX1+XbaP56MsJUAyZP4qHEmEALtVAjzjbF24HaNrWSyjoN6HbYXya7En2zoXzXvhqXO9Zrtz3eAM2El3w0Fj3DsQNCBEd4XVjATZCudnd9Vr9NtSSyMe0h0F3szhUdhZXAhFxExRkoPe4TSUMBelvGBoKU2ZE8u6VKD0xiTJjIN9QVXL0bQCENdqd4bY96qVOiVW/NtC04gXqJlE1IlD8IVAycuhOyQ7ufWW1wpf1lsYteWZXC8TamcLBdeEO1MsxqT/CNEvXLCrDnt X-Microsoft-Antispam: UriScan:;BCL:0;PCL:0;RULEID:;SRVR:BY1PR0301MB0903; X-Microsoft-Antispam-PRVS: X-Exchange-Antispam-Report-Test: UriScan:; X-Exchange-Antispam-Report-CFA-Test: BCL:0;PCL:0;RULEID:(601004)(2401047)(5005006)(8121501046)(10201501046)(3002001);SRVR:BY1PR0301MB0903;BCL:0;PCL:0;RULEID:;SRVR:BY1PR0301MB0903; X-Microsoft-Exchange-Diagnostics: 1;BY1PR0301MB0903;4:jyyuV9MXPRh4daLWzk/gqP7lsQ6YQxVDlQ/x+ueu/d4RUZ4F9upV/zxUkdBJmwG2+0qVLHwChjJhKeRzYiSf5AQpqayvC9q6zgvS88GO5KA0smkdbJ/ZtamMme831vOPa0Pbi7ghsl12XVfK3cRYWRj0BpFPQV5smGvgUmzbNBJe/8IS9OHnEygXD4Z5zO9Idi6nY7bmrgEvwdVq7xoA9Z4/Ws6asIpdIS2VDm1ctqNfexg+BaHRXMcZYxlC1BQRFQw7grVBhtiR0OCGvl1kUNobB8kQcxnNeYzh4AIx7HBsQrw9bPkEMZ54qJ+ttjBBv8BWxhjBm9PXOHVbM4b8GpQo/czYvT65rACunvX19Wb3qMCftIlHZUYLu6k/mMxF X-Forefront-PRVS: 08457955C4 X-Forefront-Antispam-Report: SFV:NSPM;SFS:(10019020)(6049001)(6009001)(485384002)(52314003)(43784003)(377454003)(24454002)(479174004)(76104003)(88552002)(4001350100001)(5004730100002)(90282001)(66066001)(92566002)(5001960100002)(65956001)(4326007)(189998001)(87976001)(117156001)(110136002)(2906002)(89122001)(3480700003)(42186005)(33656002)(3846002)(36756003)(75432002)(76176999)(47776003)(2950100001)(77096005)(23676002)(93886004)(65816999)(83506001)(5008740100001)(54356999)(50986999)(19580395003)(80316001)(40100003)(87266999)(65806001)(15975445007)(64126003)(122386002)(1096002)(50466002)(230700001)(19580405001)(586003);DIR:OUT;SFP:1102;SCL:1;SRVR:BY1PR0301MB0903;H:[192.168.0.75];FPR:;SPF:None;MLV:sfv;LANG:en; X-Microsoft-Exchange-Diagnostics: =?utf-8?B?MTtCWTFQUjAzMDFNQjA5MDM7MjM6dTVNUWo0b0hPcmdsNGlDdjZOOGdMRXVT?= =?utf-8?B?bkhINlpnY2djZEc0RjRLaVBKaGVKTGJEVmlSNjV2Rnk0ZElzbWdWcUVmU1lU?= =?utf-8?B?a3MrUXUwVjhJemxuUW1zWFB2NFdKVXRhREIzSHZ5R3lQbUxOaUpZclBaa0Zt?= =?utf-8?B?RHEwQkJIN3AzNnhhL2hRT3NpaWpCaERRSTR4ZTdZWmFMN2ZxN0dQc01lNzRK?= =?utf-8?B?VEZJQkpEVXE3Vm9NRDhmb0JrYzBSdzE4N0puMXZXc2pvNmwzUll3ZG13dTFm?= =?utf-8?B?NFBKdnNDM3BzWmNzM2VwdWJjOGRrbkEzcUZoYUl5bmVsYmZoZklrUjVCSWlQ?= =?utf-8?B?L0VMeTl0R1dZVDhBOHZSRFp6ellzeHZUeWE3YU1DOFlwT1V2SlVNUW5kUUxy?= =?utf-8?B?bnBsMVZxWmpvRHNLdTRobjQybE1HWVZZVnFYTGFicThZclFhL3VHMTBIdXlP?= =?utf-8?B?bTRYanVqRTZVcm1xWmJLb3FjRHFWUUJqN0hXUW5yWklRcmY2aGxCM3FkU24x?= =?utf-8?B?QTRQZXpmemtNZ2lXelVZTXVYZmlueWhmM3pLTm9ncUMvd25IOG1hS0h4OS9O?= =?utf-8?B?NU1FdFllK2dMdUtZSmtBWDJFcERlRXNXcVdHVDNLaTdKRlZ1UDVnY1BPTkZZ?= =?utf-8?B?RUlLOEpuNEdiSlplaHNHVk1vUlIzSURjRVNRc0Flc08yaS9HeDZ3alY1MWd4?= =?utf-8?B?Q1YydUV1eXBlSDNBQk9Ma1ArTHVqRk9oT2FDeGptU0xYYllPbWpsK3pSRXJL?= =?utf-8?B?dkg4blg2Ni82T2JrNmd2UjFaYkp2K1dmYkswN1JFZ3R5Nmo0TEpiQzJNMzJB?= =?utf-8?B?US91MWxnU0xWbGVhTXYrSkt3QUt5MFV3S3pQcks0VU5wMVVBb3FWUjlGMFYz?= =?utf-8?B?UVdnN2Jucmg3Qm92Z3pvdHJ6WHVLaWdqU1MvSEVtZUovNDgreDhXQ2UvOFdS?= =?utf-8?B?UGw0TmpzSnFxL2FaNFlYRlBMVHpsYjRJckZKdnA5VWY2SHhSclZCRlNEL3lK?= =?utf-8?B?bm5pWURxWHBCK1YrTGF2N1FLbHlFOVVEWHRuTUZ1UlNLbjNzOGh3SEFjRVNF?= =?utf-8?B?OE9tdWs2c2VmSHdDL3FwanZsRTh4L0lyTWJhcDZ3U2FFeitCQ2srbXZtNTc5?= =?utf-8?B?ZzdPbVppMk0xNU15NGU0cERORVRUcERzcjU4a2J0Vk80NjlHOVk4UUJtQmsx?= =?utf-8?B?c2NIcGU1M2ZVV0MwZ3d4clM3UEdDaU0zNnF0WHU5QWxSMmx6Q08xR2YrMGx2?= =?utf-8?B?aXZ6RXFXTXRpOE1Pc2dLdjVQWTVuM3FtdnJMbEU4eG41R0N2Y1Q0NXh1WlN0?= =?utf-8?B?QTArSTZnN1pVQ0Q3UjhienlQb3BQemNPRXJ3MDhIY0czNnh4QXE2VVhRWFFu?= =?utf-8?B?OWVyS1NRV0JDUXkyRW9iUVVlQ2ZHRmh2RzNWR0lqOHBFRDB2UkFZcnRwTU1h?= =?utf-8?B?bVB0UTgyOGN5ZVJlL3U2SHlDT0Y3c3NLckVXOGNMcUpzdVhBc0hrVDRobW1m?= =?utf-8?B?M0VVLzQ2SHljN2grRllPczNwa2ZTUTNkdzRra0VFVmU5cUV3b2dacnRqVEtx?= =?utf-8?B?YmIwVmVXZFh6cjI1RG1mb2RJd083WXJjT0NkaVphM2hyNEZCRjczbzlQUGRK?= =?utf-8?B?eDJrNFA1SnU3dzNzZzcwdGdsVGJMbnV4Rkx5ek5LZ2Fpdk5MVlZYaGJSK0gv?= =?utf-8?B?UkFEbGFPdEoxbXIxdlgwYlIxRnNnTU9aMUx1Yk9PN0FwM1doNkdxQXR1L3Jw?= =?utf-8?B?TmRISGIwaVlMd21YVUpPT21yOVNzV0tQWElzRUlLdkVZNitoQkIrRlZjWU1S?= =?utf-8?B?TkdWLzNNeGY3SXhMRzdiVTZrQ0p5OVpKR1FuamJIZm5SZTFkcXQ0bEtnRGU2?= =?utf-8?Q?Tlfa/jXBLB4mk=3D?= X-Microsoft-Exchange-Diagnostics: 1;BY1PR0301MB0903;5:3MHeBtd+5Rz4Fqkyf2IMFUiObgT1pRIl+6Sz1H1ZkZX2yzDz4nj9DsyFIphh1+vyxWPuJ0Zbap0m54W4sjgLuz/45ghtYPaYlm3IDFqrUs3cyMZMT9PSp3O/yTib0jtyjHJXG+Ig+PBjnb7ZCtStrA==;24:IDL1KZtsxlrpUU3W1rNVjSOdVvar51x8yVZBI4MC3M4RPLxc52XisgOou3GE8qLoxrRfVSVdfvR3EHXQBk7lAlp+9jAP2KM3WRZhh4vNvq4= SpamDiagnosticOutput: 1:23 SpamDiagnosticMetadata: NSPM X-OriginatorOrg: colorado.edu X-MS-Exchange-CrossTenant-OriginalArrivalTime: 07 Feb 2016 00:03:03.5556 (UTC) X-MS-Exchange-CrossTenant-FromEntityHeader: Hosted X-MS-Exchange-Transport-CrossTenantHeadersStamped: BY1PR0301MB0903 X-SW-Source: 2016-q1/txt/msg00005.txt.bz2 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 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 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. >>>> >>> > >