Re: BLAS binding

[Gautier <gautier.demontmollin@maths.unine.ch>]

Marvelous! When I find some time I'll use your binding.
For now I'm using an excerpt, of course also with "digits <>"
genericity for float type and with hastly conception - see below.
Maybe the "BLAS_Precision" could be added as generic parametre
for equilibirum with the generic "digits <>" ? Another solution
would be to pass as generic:

 1) the type you want to use,            type Float_Type is digits <>;
 2) the type which you know is single,   type Single is digits <>;
 3) the type which you know is double    type Double is digits <>;

and automatically identify the type with the 'digits or another more
subtle criterion. Surely for Ada 95, the Interface.Fortran package
contains 2) and 3) so you can determine the "BLAS_Precision" only
with 1) as generic, but then you'd lose the Ada 83 compatibility.

Gautier
===========

-- Linear Algebra routines for Sparse package, mapped to BLAS library

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

package SparseB1 is
  procedure Copy( u: in vector; v: out vector );
  function "*"(u,v: vector) return real;
  procedure Add_scaled( factor: real; u: in vector; v: in out vector );
  procedure Scale( factor: real; u: in out vector );
end SparseB1;

package body SparseB1 is

  -- identify floating point type
  -- wishful thinking about float=single, long_float=double!

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

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


  procedure scopy(n: natural; x: in  vector; incx: integer;
                              y: out vector; incy: integer);
  procedure dcopy(n: natural; x: in vector; incx: integer;
                              y: out vector; incy: integer);
  pragma Interface(FORTRAN, scopy);
  pragma Interface(FORTRAN, dcopy);

  function sdot(n: natural; x: in vector; incx: integer;
                            y: in vector; incy: integer) return real;
  function ddot(n: natural; x: in vector; incx: integer;
                            y: in vector; incy: integer) return real;
  pragma Interface(FORTRAN, sdot); 
  pragma Interface(FORTRAN, ddot);

  procedure saxpy(n: natural; a: real; x: in vector; incx: integer;
                  y: in out vector; incy: integer);
  procedure daxpy(n: natural; a: real; x: in vector; incx: integer;
                  y: in out vector; incy: integer);
  pragma Interface(FORTRAN, daxpy);
  pragma Interface(FORTRAN, saxpy);

  procedure sscal(n: natural; a: real; x: in out vector; incx: integer);
  procedure dscal(n: natural; a: real; x: in out vector; incx: integer);
  pragma Interface(FORTRAN, sscal);
  pragma Interface(FORTRAN, dscal);


  procedure Copy( u: in vector; v: out vector ) is
     begin
       if is_single then
         scopy(u'length, u,1,v,1);
       elsif is_double then
         dcopy(u'length, u,1,v,1);
       else
         v:= u;
       end if;
     end;

   function "*"(u,v: vector) return real is
     begin
       if is_single then
         return sdot(u'length, u,1,v,1);
       elsif is_double then
         return ddot(u'length, u,1,v,1); 
       else
         declare uv: real:= 0.0;
         begin
           for i in u'range loop
             uv:= uv + u(i)*v(i);
           end loop;
           return uv;
         end;
       end if;
     end;
                                             
  procedure Add_scaled( factor: real; u: in vector; v: in out vector ) is
     begin
       if is_single then
         saxpy(u'length, factor, u,1,v,1);
       elsif is_double then
         daxpy(u'length, factor, u,1,v,1);
       else
         for i in u'range loop
           v(i):= v(i) + factor * u(i);
         end loop;
       end if;
     end;
 
  procedure Scale( factor: real; u: in out vector ) is
     begin
       if is_single then
         sscal(u'length, factor, u,1);
       elsif is_double then
         dscal(u'length, factor, u,1);
       else
         for i in u'range loop
           u(i):= factor * u(i);
         end loop;
       end if;
     end;

  pragma inline("*", Add_scaled, Scale);

end SparseB1;