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.3 required=5.0 tests=BAYES_00,INVALID_MSGID
	autolearn=no autolearn_force=no version=3.4.4
X-Google-Language: ENGLISH,ASCII-7-bit
X-Google-Thread: 103376,95f6e17d8964ca3c
X-Google-Attributes: gid103376,public
From: Gautier <gautier.demontmollin@maths.unine.ch>
Subject: Re: BLAS binding
Date: 2000/05/11
Message-ID: <391ACCA0.7EA82894@maths.unine.ch>#1/1
X-Deja-AN: 622138715
Content-Transfer-Encoding: 7bit
References: <8fdrdf$pag$1@upsn21.u-psud.fr>
X-Accept-Language: en
Content-Type: text/plain; charset=us-ascii
X-Trace: 11 May 2000 17:07:15 +0100, mac13-32.unine.ch
Organization: Maths - Uni =?iso-8859-1?Q?Neuch=E2tel?=
MIME-Version: 1.0
Newsgroups: comp.lang.ada
Date: 2000-05-11T00:00:00+00:00
List-Id: <comp.lang.ada>

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;