[Scilab-Dev] eigs woes

[michael.baudin at contrib.scilab.org]

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