Thanks. I'm surprised eigs shipped with code like that in the first place. I'm going to file a bug and submit a cleaned up version of the patch. It still doesn't explain why it's more than 5 times slower on my example than the octave version while it's just a wrapper around the arpack function, but that's for another day :)
Guillaume On Tue, Nov 27, 2012 at 3:10 AM, <[email protected]> wrote: > Hi, > > I think that you are discussing the code : > > if(~isreal(AMSB)) > Lup = umf_lufact(AMSB); > [L, U, p, q, R] = umf_luget(Lup); > R = diag(R); > P = zeros(nA, nA); > Q = zeros(nA, nA); > for i = 1:nA > P(i,p(i)) = 1; > Q(q(i),i) = 1; > end > umf_ludel(Lup); > else > [hand, rk] = lufact(AMSB); > [P, L, U, Q] = luget(hand); > ludel(hand); > end > > extracted from eigs. > > The lufact function is designed to be used with lusolve. > The luget function was designed essentially for debugging purposes, > or to communicate with other algorithms. > In general, it should not be used to compute the solution. > > The other problem is the for loop, which should be avoided, because > with large sparse matrices, this will fail for performance reasons. > I'm not sure, but this loop should be vectorisable quite easily. > > In general, the umf functions are much faster, as shown by B. Pinçon. > But there might another technical reason here that I do not see. > > Best regards, > > Michaël > > Le 2012-11-26 05:53, Guillaume Horel a écrit : > >> This is a rather long email about my fight with the eigs function in >> scilab. It might be better suited for a bug report, but I wanted to >> try out this list first. >> >> It's a boundary problem for the helmholtz equation: >> http://bpaste.net/show/60377/ [1] >> >> >> On scilab-5.4.0, the code fails with the following error message: >> >> eigs: Impossible to invert complex sparse matrix. >> at line 333 of function speigs called by : >> at line 112 of function eigs called by : >> [D V] = eigs(A, [], M2,'SM'); >> at line 54 of exec file called by : >> exec('/home/guillaume/test-**eigs.sce', -1) >> >> In the code of the eigs function, it turns out that there is a test >> to check if the factors of the LU decomposition of A-sigma I are >> complex (which is the majority of cases if you start from a complex >> matrix), and the code fails with this error. Is there any reason to >> have such a test in there? >> I also realized that the code was computing the inverse by actually >> calling inv(L) and inv(U), which completely defeats the purpose of >> doing a LU decomposition. Long story short, the following patch fixes >> the two aforementionned issues: http://bpaste.net/show/60375/ [2] >> >> >> I still have two questions: >> - I'm still unsure about why the code needs two different lu solvers: >> lufact and mf_lufact. Unless lufact is really faster than umf_lufact >> for real numbers, I think that just using umf_lufact should be enough >> and would further simplify the code. >> - after applying this patch scilab computes my eigenfunctions fine. >> However the performance is pretty disappointing. On my laptop, scilab >> takes around 15s to compute 100 eigenvalues, but on the same machine >> octave takes less than 3s. I checked the octave eigs function and it's >> all written in C. However I would think the computation time is >> dominated by the calls to the various arpack functions and lu solvers >> rather than all the plumbing, but maybe I'm wrong... Any thoughts? >> >> Thanks for your insights, >> Guillaume >> >> >> Links: >> ------ >> [1] http://bpaste.net/show/60377/ >> [2] http://bpaste.net/show/60375/ >> >> ______________________________**_________________ >> dev mailing list >> [email protected] >> http://lists.scilab.org/**mailman/listinfo/dev<http://lists.scilab.org/mailman/listinfo/dev> >> > > ______________________________**_________________ > dev mailing list > [email protected] > http://lists.scilab.org/**mailman/listinfo/dev<http://lists.scilab.org/mailman/listinfo/dev> >
_______________________________________________ dev mailing list [email protected] http://lists.scilab.org/mailman/listinfo/dev
