[Scilab-Dev] eigs woes

[Guillaume Horel]

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, <michael.baudin at contrib.scilab.org> 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
>> dev at lists.scilab.org
>> http://lists.scilab.org/**mailman/listinfo/dev<http://lists.scilab.org/mailman/listinfo/dev>
>>
>
> ______________________________**_________________
> dev mailing list
> dev at lists.scilab.org
> http://lists.scilab.org/**mailman/listinfo/dev<http://lists.scilab.org/mailman/listinfo/dev>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/dev/attachments/20121127/c0cb4522/attachment.htm>