Re: ANN: Ada 2005 Math Extensions, 20100901 release

[jonathan <johnscpg@googlemail.com>]

On Sep 7, 1:56 pm, sjw <simon.j.wri...@mac.com> wrote:
> On Sep 5, 1:32 am, jonathan <johns...@googlemail.com> wrote:
>
>
>
> > On Sep 4, 12:12 pm, Ada novice <ycalleecha...@gmx.com> wrote:
>
> > > I tried compiling test_extensions.adb but I get an error message:
>
> > >     20. with Ada.Numerics.Generic_Complex_Arrays;
> > >     21. with Ada.Numerics.Generic_Arrays;
> > >              |
> > >         >>> "Ada.Numerics.Generic_Arrays" is not a predefined library
> > > unit
>
> > > ...
> > > YC
>
> > Don't worry, I got the same error message, knew I would
> > not be smart enough to fix it, and just replaced all
> > the dots with underscores. Now my version is named
> > ada_numerics_generic_arrays.adb
>
> > I got around to attaching my suite of test matrices to
> > the new LAPACK binding.  A few rambling observations..
> > First the good news is that the binding seems to work
> > the way it's supposed to .. no problems to report there.
> > This is good news because it is much easier to fill in
> > missing LAPACK routines now that Simon has shown us how
> > to do it. In particular I seem to recall that you (YC)
> > requested a routine for generalized Eigen systems:
> >    A*v = lambda*B*v
> > for matrices A and B. That should be easy to add to
> > the present interface. The not-so-good news is that I
> > tested Lapack's eigendecomposition for real general
> > (non-symmetric) matrices, dgeev, and I was surprized
> > to see it fail catastrophically on 3
> > or 4 matrices in the test suite, and perform poorly
> > on another 7 or 8.  It worked fine on small matrices
> > (N x N = 16 x 16) but degraded rapidly for larger N,
> > so I don't think I am misinterpreting the test results.
>
> I see similar behaviour (can't go up to 787 on this Mac, or rather I
> haven't tried too hard) but it starts showing up at 64x64 and is
> obvious at 128x128.
>
> Your test says that Eig.Eigensystem overwrites the input matrix -- it
> doesn't (of course, at the cost of another on-stack copy :-)

My mistake. I copied the warning from Peters_Eigen
and from the original dgeev routine, both of which
do overwrite A. Other comments were also foolishly
retained from the original Peters_Eigen_tst_1, so
I have revised the file.

> Also, I now find myself confused about pragma Convention (Fortran) --
> my code expects (effectively) pragma Convention (Ada) for matrices,
> and commenting-out your pragma made no difference to the results!!!

The pragma Convention was also from Peters_Eigen_tst_1
where it did seem to matter.  I left it in
numerics_tst_1 so that an interested user can toggle
it on and off to see if it matters .. it's easy to
forget about this option, and in some programs it can
amount to a factor of 2 in running time.

> My plan, not being a mathematician and not having any immediate need
> for matrix manipulations myself, was to implement interfaces to
> existing libraries, specifically LAPACK 'because it's there', on the
> assumption that they would be reasonably good (both in performance and
> accuracy). If that's not so, a new horizon beckons ...

Lapack is the gold standard in this business.
The reason I was grateful that you wrote this for
us is that I knew I would never be interested
in writing my own versions of 95% of the routines
in Lapack. Of the matrices that dgeev.f failed
on, most are so pathological and badly scaled
that I doubt anything of the sort would occur
in a real application, and I doubt anyone would
care much if dgeev.f failed on them. There was one
interesting exception to this rule: Lower_Integers,
which is (nearly) lower triangular, with the
integers 1, 2, 3... on the diagonal. This
is (nearly) the easiest matrix in the test suite.
Don't know what went wrong with dgeev here, but
in this case the failure was appalling. I suspect that
making an eigensolver that never fails is considered
a holy grail problem in numerical linear algebra. (Come
to think of it, I wonder if anyone can break Peters_Eigen
with its Jacobi preconditioner enabled .. or disabled for
that matter.  Let them try;-)

Mainly I wrote numerics_tst_1.adb to show how to use
the new LAPACK binding in a significant program.  I thought
it might make the new software easier to use and understand.
Well it helped me anyway!

J.