Le 02/06/2013 19:28, Volker Braun a écrit :
Sparse Matrices are a big problem for people like me working in the field of numerical solution of PDEs (elliptic and parabolic ones); and nowadays, there are some interesting researches.On Sunday, June 2, 2013 9:01:01 AM UTC+1, Charles Bouillaguet wrote:There is a presumably standard sparse-blas API : http://math.nist.gov/spblas/ Yes, though it doesn't seem to mandate any matrix storage format. So apparently you can't let it run on a given chunk of memory but you need to trust its internal (implementation-dependent) storage. Does anybody have any opinion on the performance of sparse blas and the reference implementation?
-> Whatever the data structure you use (and CSR or CSR like structures are used by everybody), the main problem is that the algorithms show a poor numerical intensity; let me explain: linear solvers are generally iterative krylov type methods (-> evaluate a polynomial of the matrix times a vector): the arithmetic intensity is the ratio:
q=bytes transfered between ram and cache /number of floating points performed
q is very low with sparse matrices and the performances are typically bounded by memory bandwidth to much less than 1 gigaflop. To improve this, the best is to modify the algorithms (a bit too long to explain here!) to try to make more operations by byte transfered: I think nobody knows a miraculous data structure, and so we keep the CSR.
-> Ok, this is the "scientific computing" perspective, where you want to solve say Poisson equation in 3d with 10^9 unknowns. But I am not sure it is what we want to do in Sage, at least nowadays; you can solve "moderate size" systems with sparse matrices (say a 2d Poisson equation -- there are 5 non zero terms by line-- with 10^6 unknowns using SuperLU which is a *direct* method (LU adapted to sparse matrices), which uses CSR structures, and it will work very well.
-> In that case, I recall that the main problem will be actually to *create* the CSR structure: you want to enter non zero coefficients (i,j)-> a_ij in any order. The mechanism used by by scipy (and thus sage) is very very slow.... much too slow (may be it changed... this must be tested).
-> Making CSR structures for other sets of coefficients is not difficult (at least in C++), if the coefficients have the same "sizeof".
t.d. -- You received this message because you are subscribed to the Google Groups "sage-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+unsubscr...@googlegroups.com. To post to this group, send email to sage-devel@googlegroups.com. Visit this group at http://groups.google.com/group/sage-devel?hl=en. For more options, visit https://groups.google.com/groups/opt_out.
<<attachment: tdumont.vcf>>
