NEC-LIST: Debugging a LAPACK version of nec2

From: Tim Molteno <tim_at_email.domain.hidden>
Date: Mon, 14 Jan 2008 09:44:45 -0500

Hi Everyone

I am the author of nec2++ (a translation/rewrite of nec2 into C++). I have
been struggling for a while with the original LU decomposition routines.

I have begun to get little suspicious of the FACTRS routine (called
lu_decomposition() in the C++ version).

In the nec2++ source (http://alioth.debian.org/projects/necpp/), there is a
file called atlas_test.cpp. It contains several versions of the routines. One
translated directly from the nec2c code, another from the fortran of burke
(see the intro comments in matrix_algebra.cpp). Basically for small test
arrays the original code does not do a correct LU decomposition (when
compared against octave)

For example, the following 4x4 matrix

    3 1 -4 2
    3 1 0 2
    2 13 -1 0
   -2 3 -1 4

produces

l =

   1.00000 0.00000 0.00000 0.00000
   1.00000 1.00000 0.00000 0.00000
   0.00000 -1.00000 1.00000 0.00000
   2.00000 -1.33333 0.00000 1.00000

u =

    3.00000 1.00000 0.66667 -0.66667
    0.00000 12.33330 0.00000 0.29730
    0.00000 0.00000 -4.00000 0.17568
    0.00000 0.00000 0.00000 5.72973

And these factors are incorrect (can not be multiplied together with a
permutation matrix to get the original matrix back)..

I have gone back to the original book referred to in the source (Ralston) but
I havent found the problem with the FACTRS code.

After substituting LAPACK routines (now in the CVS archive for nec2++), the
testbench gives slightly different results. Typically between 0.1 and 5%
different. It gives them a lot faster too! There might be examples that
produce dramatically different results.

If anobody has experience with this piece of code, or with the effect of using
LAPACK instead of the FACTRS code, I would love to know whether you ran into
this problem.

Kind Regards,

Tim Molteno
Department of Physics
University of Otago
Dunedin, New Zealand

-- 
The NEC-List mailing list
NEC-List_at_robomod.net
http://www.robomod.net/mailman/listinfo/nec-list
Received on Mon Jan 14 2008 - 14:45:03 EST

This archive was generated by hypermail 2.2.0 : Sat Oct 02 2010 - 00:10:46 EDT