Tim - thanks for the reference. The paper will come in handy. This is a longstanding issue, that we just haven’t got around to addressing yet, but perhaps now is a good time.
https://github.com/JuliaLang/julia/issues/3295 We have a very simplistic factorize() for sparse matrices that must have been implemented as a stopgap. This is what it currently does and that explains everything. # placing factorize here for now. Maybe add a new file function factorize(A::SparseMatrixCSC) m, n = size(A) if m == n Ac = cholfact(A) Ac.c.minor == m && ishermitian(A) && return Ac end return lufact(A) end -viral > On 06-Jan-2015, at 1:57 am, Tim Davis <[email protected]> wrote: > > That does sound like a glitch in the "\" algorithm, rather than in UMFPACK. > The OpenBLAS is pretty good. > > This is very nice in Julia: > > F = lufact (d["M"]) ; F \ d > > That's a great idea to have a factorization object like that. I have a > MATLAB toolbox that does > the same thing, but it's not a built-in function inside MATLAB. It's written > in M, so it can be slow for > small matrices. With it, however, I can do: > > F = factorize (A) ; % does an LU, Cholesky, QR, SVD, or whatever. Uses my > polyalgorithm for "\". > x = F\b ; > > I can do S = inverse(A); which returns a factorization, not an inverse, but > with a flag > set so that S*b does A\b (and yes, S\b would do A*b, since S keeps a copy of > A inside it, as well). > > You can also specify the factorization, such as > > F=factorize(A, 'lu') > F=factorize(A,'svd') ; etc. > > It's in SuiteSparse/MATLAB_tools/Factorize, if you're interested. I've > suggested the same > feature to The MathWorks. > > My factorize function includes a backslash polyalgorithm, if you're > interested in taking a look. > > Algorithm 930: FACTORIZE: an object-oriented linear system solver for MATLAB > T. A. Davis, ACM Transactions on Mathematical Software, Vol 39, Issue 4, pp. > 28:1 - 28:18, 2013. > http://faculty.cse.tamu.edu/davis/publications_files/Factorize_an_object_oriented_linear_system_solver_for_MATLAB.pdf > > On Mon, Jan 5, 2015 at 2:11 PM, Viral Shah <[email protected]> wrote: > The BLAS will certainly make a difference, but OpenBLAS is reasonably good. > > I also wonder what is happening in our \ polyalgorithm. The profile suggests > the code is trying Cholesky decomposition, but it really shouldn't since the > matrix is not symmetric. If I just do the lufact(), which essentially calls > Umfpack, I can match Matlab timing: > > @time F = lufact(d["M"]); F \ d["RHS"]; > > -viral > > > On Tuesday, January 6, 2015 12:31:34 AM UTC+5:30, Tim Davis wrote: > The difference could be the BLAS. MATLAB comes with its own BLAS library, > and the performance > of the BLAS has a huge impact on the performance of UMFPACK, particularly for > 3D discretizations. > > On Mon, Jan 5, 2015 at 6:21 AM, Ehsan Eftekhari <[email protected]> wrote: > I'm solving diffusion equation in Matlab on a 3D uniform grid (31x32x33) and > Julia. I use the "\" to solve the linear system of equations. Here is the > performance of the linear solver in Julia: > elapsed time: 2.743971424 seconds (35236720 bytes allocated) > > and Matlab (I used spparms('spumoni',1) to see what "\" does in Matlab): > sp\: bandwidth = 1056+1+1056. > sp\: is A diagonal? no. > sp\: is band density (0.00) > bandden (0.50) to try banded solver? no. > sp\: is A triangular? no. > sp\: is A morally triangular? no. > sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? no. > sp\: use Unsymmetric MultiFrontal PACKage with automatic reordering. > sp\: UMFPACK's factorization was successful. > sp\: UMFPACK's solve was successful. > Elapsed time is 0.819120 seconds. > > I have uploaded the sparse matrix (M) and the right-hand side (RHS) vectors > in a mat file here: > https://drive.google.com/open?id=0B8OOfC6oWXEPV2xYTWFMZTljU00&authuser=0 > > I read in the documents that Julia uses Umfpack for sparse matrices. My > question is why umfpack is faster when it is called from matlab? > > The matlab and julia codes are here: > https://drive.google.com/open?id=0B8OOfC6oWXEPbXFnYlh2TFBKV1k&authuser=0 > https://drive.google.com/open?id=0B8OOfC6oWXEPdlNfOEFKbnV5MlE&authuser=0 > > and the FVM codes are here: > https://github.com/simulkade/FVTool > https://github.com/simulkade/JFVM > > Thanks a lot in advance, > > Ehsan > > > On Wednesday, June 5, 2013 8:39:15 AM UTC+2, Viral Shah wrote: > I guess it is the last 20 years of sparse solver work packed into one > function. Not many fields can boast of providing this level of usability out > of their work. :-) > > There are a class of people who believe that things like \ encourage blackbox > usage, with people doing stuff they do not understand, and there are others > who believe in standing on the shoulders of giants. > > I find that we have taken a good approach in Julia, where we have \ and it > will have the perfect polyalgorithm at some point. But, you also have the > option of digging deeper with interfaces such as lufact(), cholfact(), > qrfact(), and finally, even if that does not work out for you, call the > LAPACK and SuiteSparse functions directly. > > -viral > > On Wednesday, June 5, 2013 9:42:12 AM UTC+5:30, Stefan Karpinski wrote: > Goodness. This is why there needs to be a polyalgorithm – no mortal user > could know all of this stuff! > > > On Tue, Jun 4, 2013 at 11:11 PM, Viral Shah <[email protected]> wrote: > Doug, > > Ideally, the backslash needs to look for diagonal matrices, triangular > matrices and permutations thereof, banded matrices and the least squares > problems (non-square). In case it is square, symmetric and hermitian, with a > heavy diagonal(?), cholesky can be attempted, with a fallback to LU. I > believe we do some of this in the dense \ polyalgorithm, but I am not sure if > we look for the banded cases yet. > > This is what Octave does: > http://www.gnu.org/software/octave/doc/interpreter/Sparse-Linear-Algebra.html#Sparse-Linear-Algebra > > This is Tim's Factorize for solving linear and least squares systems: > http://www.cise.ufl.edu/research/sparse/SuiteSparse/current/SuiteSparse/MATLAB_Tools/Factorize/Doc/factorize_demo.html > > -viral > > > On Tuesday, June 4, 2013 8:18:39 PM UTC+5:30, Douglas Bates wrote: > On Thursday, May 30, 2013 10:10:59 PM UTC-5, Mingming Wang wrote: > Hi, > > I am trying to port my MATLAB program to Julia. The for loop is about 25% > faster. But the backslash is about 10 times slower. It seems in MATLAB, the > backslash is parallelized automatically. Is there any plan in Julia to do > this? BTW, the matrix I am solving is sparse and symmetric. > > For a sparse symmetric matrix try > > cholfact(A)\b > > The simple > > A\b > > call will always use an LU decomposition from UMFPACK. > > > >
