Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Replies: 10   Last Post: Nov 24, 1999 4:55 PM

 Messages: [ Previous | Next ]
 Gautier Posts: 4 Registered: 12/15/04
Posted: Nov 24, 1999 4:55 PM

> Is there much numerical analysis done in ADA?

I'd say: not much.

> Could it be (at least in theory) an alternative to Fortran?

Yes, even in practice...
Just for supercomputers, there is still a lack of compilers.

> Are there any serios performance measurments?

I did some CPU time measurments.
For writing Ada like Fortran, no problem - with a decent compiler
like DEC/Compaq or GNAT.

For using Ada features, no problem if you are careful.
An example to illustrate it, with genericity and composite type
definition - see comment after excerpt...

~~~~~~~~~~~~~
generic

type real is digits <>;
type index is range <>;
type vector is array(index range <>) of real;

package Sparse is

type index_array is array(index range <>) of index;

--------------------------------------------------
-- Define a matrix with compressed rows storage --
--------------------------------------------------

type CRS_matrix( rows_max_p1, nnz_max: index ) is
record
val: vector(1..nnz_max);
col_ind: index_array(1..nnz_max);
row_ptr: index_array(1..rows_max_p1); -- p1 means +1 -> extra row ptr
symmetric: boolean;
rows: index:= rows_max_p1-1; -- for re-use; must be in [ 1..rows_max_p1-1 ]
nnz: index:= nnz_max; -- for re-use; must be in [ 1..nnz_max ]
end record;

type p_CRS_matrix is access CRS_matrix;

type t_precond is ( none, diag, tridiag );

----------------------------------
-- Matrix-vector multiplication --
----------------------------------
~~~~~~~~~~~~~

To obtain performances equal to Fortran, you can have to
avoid indirect addressing due to the record's fields
(DEC/Compaq Ada needs it; GNAT optimizer doesn't) :

~~~~~~~~~~~~~
rows: constant index:= A.rows;
val: vector renames A.val;
col_ind: index_array renames A.col_ind;
row_ptr: index_array renames A.row_ptr;

begin

if not A.symmetric then
-- *** la matrice est memorisee sous forme non symetrique

for i in 1..rows loop
deb := row_ptr(i);
fin := row_ptr(i + 1) - 1;
wi_sum := 0.0;
for j in deb .. fin loop
wi_sum := wi_sum + val(j) * u(col_ind(j));
end loop;
w(i) := wi_sum;
end loop;
~~~~~~~~~~~~~
Identification of the generic floating-point type allows you to
use the correct BLAS routines
-- identify floating point type

is_single: constant boolean:=
real'digits = float'digits;

is_double: constant boolean:=
real'digits = long_float'digits;

All that is solved at compile-time, by instanciating the generic
package. E.g. a Sparse matrix package for double-precision is
created by this line of code
package LFSparse is new Sparse(long_float, integer, Matrices.vector);

with hard-coded double-precision operations and calls to d* routines in BLAS!

* A good usage of procedure encapsulation saves a lot of parametre passings.
* Subtyping spares a lot of dimension passings (executable size really shrinks!)
e.g.
procedure FT_elem(gelem: t_element) is
nsd: constant positive:= spatial_dimension(gelem);
subtype sp_vector is vector(1..nsd); -- vecteur spatial

* Usage of exceptions saves a lot of error codes transmission.
* A last thing: use the (standard) pragma suppress(all_checks)...

That's all for a `short' overview...

Gautier

Date Subject Author
11/17/99 Helmut Zeisel
11/17/99 P.G.Hamer
11/17/99 Helmut Zeisel
11/17/99 P.G.Hamer
11/17/99 Tabon
11/17/99 Warner Bruns
11/17/99 Warner Bruns
11/17/99 Helmut Zeisel
11/17/99 Warner Bruns
11/17/99 Tabon
11/24/99 Gautier