Re: Linear Algebra package for Ada ?

[Laurent GASSER <lgasser@freesurf.fr>]

Beware that Fortran is one of the rare language to store
matrix entries column-wise. Most of the others do it
row-wise. So you may have to work with transposed matrices
somewhere...

An older book has elementary solvers:

Jean-Etienne Rombaldi
Algorithmique  num\'erique et Ada
Masson (1994)
ISBN 2-225-84384-8

This is a book entierly dedicated to numerical analysis in Ada 83.
The book goes through the topics of solving linear systems, of
eigenvalues, of approximations and interpolations, of numerical
integration, of differential equations (ODE) and finite differences
(with an example of finite elements). The approach is rather about
introducing the theory and the usage of the code on disc. The book is
in French, with a code not so independent of the compiler OpenAda for
PC-DOS. A very complete library, and I know that the author (now in
Corse) works at making a Window interface to it.

I played a little with conjugate gradient algorithms in Ada. With
the proper declarations of operators, you end up writing exactly
the formulas found in math books. I wish I had the time to make
a true package out of it. 

Surely sometimes in the future, when I have spare time... ;-)

I also found in my old mails the following:

Title: Ada-95 for numerics?
Date: April 2, 1996
From: rodemann@mathematik.uni-ulm.de
Newsgroup: com.lang.ada

Hello again!
Thanks for all your postings. First I want to make the following clear:
I am not criticizing the GNAT compiler --- sadly it is the only Ada-95 
compiler I am able to use by now. And it is a fine thing to use.
Since I am now working at a theoretical physics department and lots of
the software written here is still done in Fortran-77. Many of our problems
need quite a lot computing time. I guess for some applications, especially
matrices Fortran may be still the best choice. At least for there are already
many routines available (Numerical Recipes, Nag, EISPACK,...). But there
are other applications that may be done better in another language...perhaps
Ada-95 (complex and dynamic data structures, usage of generics and OO-
techniques). So I started coding some of _my_ problems in Ada-95 and check
for the speed. Yet, I am still a beginner in Ada-95, so maybe the coding style
is somewhat unlucky... 
Some of you wanted to look at my code, so I appended it to this article. 
Thanks and greetings
George
-----------------------------------------------------------------------------
--
--                                matrices.ads
--
--            A generic matrix package (only a few functions shown...)
--
-----------------------------------------------------------------------------
generic
   type Floating is digits <>;
package Matrices is
   type Matrix is array (integer range <>, integer range <>) of Floating;
   pragma Convention (Fortran, Matrix);
   function Unit (n : integer) return matrix;
   procedure Get (a : out Matrix);
   procedure Put (a : Matrix);
end Matrices;
-----------------------------------------------------------------------------
--
--                         eigensystems.ads
--
--      A generic library of routines taken from "Numerical Recipes in C" 
-- The code ist strongly oriented at the C version, only the generic character
-- has been added.
--
-----------------------------------------------------------------------------
with vectors;
with matrices;
generic
   type Floating is digits <>;
   with function sqrt (x : Floating) return Floating;
   with package v is new vectors  (Floating);
   with package m is new matrices (Floating);
   use v, m;
package eigensystems is
   procedure tred2 (a : in out Matrix; d, e : out vector);
end eigensystems;
-----------------------------------------------------------------------------
--
--                         vectors.ads
--
-----------------------------------------------------------------------------
generic
   type Floating is digits <>;
package Vectors is
   type Vector is array (integer range <>) of Floating;
   
   procedure Get (x : out Vector);
   procedure Put (x : Vector);
  
end Vectors;
-----------------------------------------------------------------------------
--
--                             matrices.adb
--
-- A generic matrix package (only a few functions shown here...)
--
-----------------------------------------------------------------------------
with Ada.Text_Io;
use Ada.Text_Io;
package body Matrices is
   function Unit (n : Integer) return Matrix is
      C : Matrix (1..N, 1..N);
   begin
      for i in C'range(1) loop
         for j in C'first(2)..i loop
            C(i,j) := 0.0;
            C(j,i) := 0.0;
         end loop;
         C(i,i) := 1.0;
      end loop;
      return (C);
   end Unit;
            
   package Floating_Io is new Float_io (Floating); use Floating_Io;
   procedure Get (a : out matrix) is
   begin
      for i in a'Range(1) loop
         for j in a'Range(2) loop
            Get (a(i,j));
         end loop;
      end loop;
   end Get;
   procedure Put (a : matrix) is
   begin
      for i in a'Range(1) loop
         for j in a'Range(2) loop
            Put (a(i,j)); New_Line;
         end loop;
         New_Line;
      end loop;
   end Put;
end Matrices;
-----------------------------------------------------------------------------
--
--                             mytest.adb
--
-- Main program: read in a symmetric matrix and brings it to tridiagonal form
--
-----------------------------------------------------------------------------
with Vectors;
with Matrices;
with Eigensystems;
with Ada.Text_Io;
with Ada.Numerics.Aux;
use Ada.Text_Io;
use Ada.Numerics.Aux;
procedure mytest is
   package MyVectors  is new Vectors  (Double); use MyVectors;
   package MyMatrices is new Matrices (Double); use MyMatrices;
   package MyEigensystems is new 
      Eigensystems (Double, Sqrt, MyVectors, MyMatrices);
   N : integer := 5;
   package Int_Io is new Integer_Io(integer); use Int_Io;
begin
   Get (N);
   declare
      a, b : Matrix (1..N,1..N);
      d, e : Vector (1..N);
   begin
      Get (b);
      for i in 1..100000 loop
         a := b;
         MyEigensystems.Tred2 (a, d, e);
      end loop;
      Put (d);
      Put (e);
   end;
end mytest;
-----------------------------------------------------------------------------
--
--                       eigensystems.adb
--
--           Contains the tred2 routine from "Numerical Recipes in C"
--
-----------------------------------------------------------------------------
package body eigensystems is
   pragma Optimize (Time);
   pragma Suppress (All_Checks);
   procedure tred2 (a : in out Matrix; d, e : out vector) is
      scale : Floating;
      hh    : Floating;
      h     : Floating;
      g     : Floating;
      f     : Floating;
      n     : constant integer := a'Last(1);
      l     : integer range 0..n;
      i, j, k : integer range 1..n;
   begin
      for i in reverse 2..n loop
         l := i-1;
         h := 0.0;
         scale := 0.0;      
         if l > 1 then
            for k in 1..l loop
               scale := scale + abs (a(i,k));
            end loop;
            if scale = 0.0 then 
              e(i) := a(i,l);
            else
               for k in 1..l loop
                  a(i,k) := a(i,k) / scale;
                  h := h + a(i,k) * a(i,k);
               end loop;
               f := a(i,l);
               if f >= 0.0 then
                  g := -sqrt(h);
               else
                  g := sqrt(h);
               end if;
               e(i) := scale * g;
               h := h - f*g;
               a(i,l) := f - g;
               f := 0.0;
               for j in 1..l loop
                  a(j,i) := a(i,j) / h;
                  g := 0.0;
                  for k in 1..j loop
                     g := g + a(j,k) * a(i,k);
                  end loop;
                  for k in j+1..l loop
                     g := g + a(k,j)*a(i,k);
                  end loop;
                  e(j) := g / h;
                  f := f + e(j) * a(i,j);
               end loop;
               hh := f / (h+h);
               for j in 1..l loop
                  f := a(i,j);
                  e(j) := e(j) - hh*f;
                  g := e(j);
                  for k in 1..j loop
                     a(j,k) := a(j,k) - (f*e(k) + g*a(i,k));
                  end loop;
               end loop;
            end if;
         else
            e(i) := a(i,l);   
         end if;
         d(i) := h;
      end loop;
      d(1) := 0.0;
      e(1) := 0.0;
      for i in 1..n loop
         l := i - 1;
         if d(i) /= 0.0 then
            for j in 1..l loop
               g := 0.0;
               for k in 1..l loop
                  g := g + a(i,k)*a(k,j);
               end loop;
               for k in 1..l loop
                  a(k,j) := a(k,j) - g*a(k,i);
               end loop;
            end loop;
         end if;
         d(i) := a(i,i);
         a(i,i) := 1.0;
         for j in 1..l loop
            a(j,i) := 0.0;
            a(i,j) := 0.0;
         end loop; 
      end loop;
   end tred2;
end eigensystems;
-----------------------------------------------------------------------------
--
--                         vectors.adb
--
-- A generic Vector package (This is just a stripped version
--
-----------------------------------------------------------------------------
with Ada.Text_Io;
use Ada.Text_Io;
package body vectors is
   
   package Floating_Io is new Float_Io (Floating); use Floating_Io;
   procedure Get (x : out Vector) is
   begin
      for i in x'Range loop
         Get (x(i));
      end loop;
   end Get;
   procedure Put (x : Vector) is
   begin
      for i in x'Range loop
         Put (x(i)); New_Line;
      end loop;
   end Put;
end vectors;

Robert Dewar wrote:
> 
> In article <8dv75n$qi1$1@news.uit.no>,
>   Reinert.Korsnes@npolar.no wrote:
> > Hi,
> >
> > Could anybody recommend me an Ada95 package for linear algebra
> > (including manipulation of vectors/matrices ?)
> 
> Why not use the interface to Linpack (is that the right spelling
> -- I refer to the standard Fortran linear algebra package)
> 
> Sent via Deja.com http://www.deja.com/
> Before you buy.

-- 
Laurent GASSER (gasser@hirondelle-lumineuse.ch)
http://www.hirondelle-lumineuse.ch
Computers do not solve problems, they execute solutions.