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