Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.


Gautier
Posts:
4
Registered:
12/15/04


Re: ADA?
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_p11;  for reuse; must be in [ 1..rows_max_p11 ] nnz: index:= nnz_max;  for reuse; must be in [ 1..nnz_max ] end record;
type p_CRS_matrix is access CRS_matrix;
type t_precond is ( none, diag, tridiag );
  Matrixvector 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 floatingpoint 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 compiletime, by instanciating the generic package. E.g. a Sparse matrix package for doubleprecision is created by this line of code package LFSparse is new Sparse(long_float, integer, Matrices.vector);
with hardcoded doubleprecision 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



