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=-1.9 required=5.0 tests=BAYES_00 autolearn=ham autolearn_force=no version=3.4.4 X-Google-Language: ENGLISH,ASCII-7-bit X-Google-Thread: 103376,35ed85a9d92a025 X-Google-Attributes: gid103376,public From: rodemann@mathematik.uni-ulm.de (Joerg Rodemann) Subject: Ada-95 for numerics? Date: 1996/04/02 Message-ID: <4jqofv$jtp@rigel.rz.uni-ulm.de> X-Deja-AN: 145387693 organization: University of Ulm, SAI, Germany newsgroups: comp.lang.ada Date: 1996-04-02T00:00:00+00:00 List-Id: 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;