From mboxrd@z Thu Jan 1 00:00:00 1970 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on polar.synack.me X-Spam-Level: X-Spam-Status: No, score=-0.3 required=5.0 tests=BAYES_00,FREEMAIL_FROM, REPLYTO_WITHOUT_TO_CC autolearn=no autolearn_force=no version=3.4.4 X-Google-Language: ENGLISH,ASCII-7-bit X-Google-Thread: 103376,85917c9414146b3f X-Google-Attributes: gid103376,public From: Laurent GASSER Subject: Re: Linear Algebra package for Ada ? Date: 2000/04/24 Message-ID: <3904ACD0.3FE0EEDC@freesurf.fr> X-Deja-AN: 615108692 Content-Transfer-Encoding: 7bit References: <8dv75n$qi1$1@news.uit.no> <8e05ir$lc$1@nnrp1.deja.com> Content-Type: text/plain; charset=us-ascii X-Complaints-To: abuse@fr.clara.net X-Trace: nnrp3.clara.net 956607400 212.43.225.152 (Mon, 24 Apr 2000 21:16:40 BST) Organization: Hirondelle Lumineuse MIME-Version: 1.0 Reply-To: gasser@hirondelle-lumineuse.ch NNTP-Posting-Date: Mon, 24 Apr 2000 21:16:40 BST Newsgroups: comp.lang.ada Date: 2000-04-24T00:00:00+00:00 List-Id: 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.