* BLAS binding
@ 2000-05-11 0:00 Duncan Sands
2000-05-11 0:00 ` Gautier
0 siblings, 1 reply; 2+ messages in thread
From: Duncan Sands @ 2000-05-11 0:00 UTC (permalink / raw)
I'm interested in doing a medium-thin binding to the BLAS
(Basic Linear Algebra Subprograms) and to some of LAPACK.
I've written a (preliminary) binding to the level 1 BLAS
to show the kind of thing I have in mind (see below or
http://topo.math.u-psud.fr/~sands/BLAS/ ): you just
instantiate it with your favorite floating point type
and it automagically calls the BLAS routines of the
appropriate precision. I added checking of array bounds
(easy in Ada) and exploited overloading to give simplified
forms of the procedures.
Please let me know if you have any comments. Help is
welcome!
Duncan Sands.
with Interfaces.Fortran;
generic
type Float_Type is digits <>;
package Real_BLAS is
pragma Elaborate_Body (Real_BLAS);
type Precision is (Double, Single, Unsupported);
BLAS_Precision : constant Precision;
Argument_Error : exception;
Unsupported_Precision_Error : exception;
type Vector is array (Integer range <>) of Float_Type'Base;
pragma Convention (Fortran, Vector);
type Matrix is array (Integer range <>, Integer range <>)
of Float_Type'Base;
pragma Convention (Fortran, Matrix);
-- Level 1
-- ROTG: Generate plane rotation
procedure ROTG (
A, B : in out Float_Type'Base;
C, S : out Float_Type'Base
);
pragma Inline (ROTG);
type Modified_Givens is new Vector (1 .. 5);
-- ROTMG: Generate modified plane rotation
procedure ROTMG (
D1, D2 : in out Float_Type'Base;
A : in out Float_Type'Base;
B : in Float_Type'Base;
PARAM : out Modified_Givens
);
pragma Inline (ROTMG);
-- ROT: Apply plane rotation
procedure ROT (
N : in Natural;
X : in out Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer;
C, S : in Float_Type'Base
);
procedure ROT (
X, Y : in out Vector;
C, S : in Float_Type'Base
);
pragma Inline (ROT);
-- ROTM: Apply modified plane rotation
procedure ROTM (
N : in Natural;
X : in out Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer;
PARAM : in Modified_Givens
);
procedure ROTM (
X, Y : in out Vector;
PARAM : in Modified_Givens
);
pragma Inline (ROTM);
-- SWAP: x <-> y
procedure SWAP (
N : in Natural;
X : in out Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer
);
procedure SWAP (X, Y : in out Vector);
pragma Inline (SWAP);
-- SCAL: x <- alpha x
procedure SCAL (
N : in Natural;
ALPHA : in Float_Type'Base;
X : in out Vector;
INCX : in Integer
);
procedure SCAL (
ALPHA : in Float_Type'Base;
X : in out Vector
);
pragma Inline (SCAL);
-- COPY: y <- x
procedure COPY (
N : in Natural;
X : in Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer
);
procedure COPY (
X : in Vector;
Y : in out Vector
);
pragma Inline (COPY);
-- AXPY: y <- alpha x + y
procedure AXPY (
N : in Natural;
ALPHA : in Float_Type'Base;
X : in Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer
);
procedure AXPY (
ALPHA : in Float_Type'Base;
X : in Vector;
Y : in out Vector
);
pragma Inline (AXPY);
-- DOT: dot <- x^T y
function DOT (
N : Natural;
X : Vector;
INCX : Integer;
Y : Vector;
INCY : Integer
) return Float_Type'Base;
function DOT (X, Y : Vector) return Float_Type'Base;
pragma Inline (DOT);
-- DSDOT: dot <- alpha + x^T y (double precision accumulation)
function DSDOT (
N : Natural;
ALPHA : Float_Type'Base;
X : Vector;
INCX : Integer;
Y : Vector;
INCY : Integer
) return Float_Type'Base;
function DSDOT (
ALPHA : Float_Type'Base;
X, Y : Vector
) return Float_Type'Base;
pragma Inline (DSDOT);
-- NRM2: nrm2 <- ||x||_2
function NRM2 (
N : Natural;
X : Vector;
INCX : Integer
) return Float_Type'Base;
function NRM2 (X : Vector) return Float_Type'Base;
pragma Inline (NRM2);
-- ASUM: asum <- ||x||_1
function ASUM (
N : Natural;
X : Vector;
INCX : Integer
) return Float_Type'Base;
function ASUM (X : Vector) return Float_Type'Base;
pragma Inline (ASUM);
-- AMAX: amax <- 1st k such that |x_k| = max(|x_i|)
function AMAX (
N : Natural;
X : Vector;
INCX : Integer
) return Integer;
function AMAX (X : Vector) return Integer;
pragma Inline (AMAX);
private
-- Strings that will be appended/prepended to each BLAS subroutine name
-- to form the external name used when importing it
Name_Prepend : constant String := "";
Name_Append : constant String := "_";
BLAS_Precision : constant Precision := Precision'Val (
Precision'Pos (Unsupported) -
(
Boolean'Pos (
Float_Type'Base'Digits = Interfaces.Fortran.Real'Digits and
Float_Type'Base'Size = Interfaces.Fortran.Real'Size
) - Boolean'Pos (False)
) - 2 *
(
Boolean'Pos (
Float_Type'Base'Digits = Interfaces.Fortran.Double_Precision'Digits
and Float_Type'Base'Size =
Interfaces.Fortran.Double_Precision'Size
) - Boolean'Pos (False)
)
);
-- Used this method to encourage compile time evaluation
-- Will raise Constraint_Error if Double and Single precision are the same
end Real_BLAS;
package body Real_BLAS is
use Interfaces.Fortran;
----------
-- AMAX --
----------
function ISAMAX (
N : Fortran_Integer;
X : Vector;
INCX : Fortran_Integer
) return Fortran_Integer;
pragma Import (Fortran, ISAMAX, Name_Prepend & "isamax" & Name_Append);
function IDAMAX (
N : Fortran_Integer;
X : Vector;
INCX : Fortran_Integer
) return Fortran_Integer;
pragma Import (Fortran, IDAMAX, Name_Prepend & "idamax" & Name_Append);
function AMAX (
N : Natural;
X : Vector;
INCX : Integer
) return Integer is
begin
-- Use of INCX rather than abs INCX is deliberate
if N > 0 and (N - 1) * INCX >= X'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
return X'First - 1 + Integer (
ISAMAX (Fortran_Integer (N), X, Fortran_Integer (INCX))
);
when Double =>
return X'First - 1 + Integer (
IDAMAX (Fortran_Integer (N), X, Fortran_Integer (INCX))
);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end AMAX;
function AMAX (X : Vector) return Integer is
begin
case BLAS_Precision is
when Single =>
return X'First - 1 + Integer (
ISAMAX (X'Length, X, 1)
);
when Double =>
return X'First - 1 + Integer (
IDAMAX (X'Length, X, 1)
);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end AMAX;
----------
-- ASUM --
----------
function SASUM (
N : Fortran_Integer;
X : Vector;
INCX : Fortran_Integer
) return Float_Type'Base;
pragma Import (Fortran, SASUM, Name_Prepend & "sasum" & Name_Append);
function DASUM (
N : Fortran_Integer;
X : Vector;
INCX : Fortran_Integer
) return Float_Type'Base;
pragma Import (Fortran, DASUM, Name_Prepend & "dasum" & Name_Append);
function ASUM (
N : Natural;
X : Vector;
INCX : Integer
) return Float_Type'Base is
begin
-- Use of INCX rather than abs INCX is deliberate
if N > 0 and (N - 1) * INCX >= X'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
return SASUM (Fortran_Integer (N), X, Fortran_Integer (INCX));
when Double =>
return DASUM (Fortran_Integer (N), X, Fortran_Integer (INCX));
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end ASUM;
function ASUM (X : Vector) return Float_Type'Base is
begin
case BLAS_Precision is
when Single =>
return SASUM (X'Length, X, 1);
when Double =>
return DASUM (X'Length, X, 1);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end ASUM;
----------
-- AXPY --
----------
procedure SAXPY (
N : in Fortran_Integer;
ALPHA : in Float_Type'Base;
X : in Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer
);
pragma Import (Fortran, SAXPY, Name_Prepend & "saxpy" & Name_Append);
procedure DAXPY (
N : in Fortran_Integer;
ALPHA : in Float_Type'Base;
X : in Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer
);
pragma Import (Fortran, DAXPY, Name_Prepend & "daxpy" & Name_Append);
procedure AXPY (
N : in Natural;
ALPHA : in Float_Type'Base;
X : in Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer
) is
begin
if N > 0 and (
(N - 1) * abs INCX >= X'Length or (N - 1) * abs INCY >= Y'Length
) then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SAXPY (
Fortran_Integer (N),
ALPHA,
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Double =>
DAXPY (
Fortran_Integer (N),
ALPHA,
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end AXPY;
procedure AXPY (
ALPHA : in Float_Type'Base;
X : in Vector;
Y : in out Vector
) is
begin
if X'Length /= Y'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SAXPY (X'Length, ALPHA, X, 1, Y, 1);
when Double =>
DAXPY (X'Length, ALPHA, X, 1, Y, 1);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end AXPY;
----------
-- COPY --
----------
procedure SCOPY (
N : in Fortran_Integer;
X : in Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer
);
pragma Import (Fortran, SCOPY, Name_Prepend & "scopy" & Name_Append);
procedure DCOPY (
N : in Fortran_Integer;
X : in Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer
);
pragma Import (Fortran, DCOPY, Name_Prepend & "dcopy" & Name_Append);
procedure COPY (
N : in Natural;
X : in Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer
) is
begin
if N > 0 and (
(N - 1) * abs INCX >= X'Length or (N - 1) * abs INCY >= Y'Length
) then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SCOPY (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Double =>
DCOPY (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end COPY;
procedure COPY (
X : in Vector;
Y : in out Vector
) is
begin
if X'Length /= Y'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SCOPY (X'Length, X, 1, Y, 1);
when Double =>
DCOPY (X'Length, X, 1, Y, 1);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end COPY;
---------
-- DOT --
---------
function SDOT (
N : Fortran_Integer;
X : Vector;
INCX : Fortran_Integer;
Y : Vector;
INCY : Fortran_Integer
) return Float_Type'Base;
pragma Import (Fortran, SDOT, Name_Prepend & "sdot" & Name_Append);
function DDOT (
N : Fortran_Integer;
X : Vector;
INCX : Fortran_Integer;
Y : Vector;
INCY : Fortran_Integer
) return Float_Type'Base;
pragma Import (Fortran, DDOT, Name_Prepend & "ddot" & Name_Append);
function DOT (
N : Natural;
X : Vector;
INCX : Integer;
Y : Vector;
INCY : Integer
) return Float_Type'Base is
begin
if N > 0 and (
(N - 1) * abs INCX >= X'Length or (N - 1) * abs INCY >= Y'Length
) then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
return SDOT (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Double =>
return DDOT (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end DOT;
function DOT (X, Y : Vector) return Float_Type'Base is
begin
if X'Length /= Y'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
return SDOT (X'Length, X, 1, Y, 1);
when Double =>
return DDOT (X'Length, X, 1, Y, 1);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end DOT;
-----------
-- DSDOT --
-----------
function SDSDOT (
N : Fortran_Integer;
ALPHA : Float_Type'Base;
X : Vector;
INCX : Fortran_Integer;
Y : Vector;
INCY : Fortran_Integer
) return Float_Type'Base;
pragma Import (Fortran, SDSDOT, Name_Prepend & "sdsdot" & Name_Append);
function DSDOT (
N : Natural;
ALPHA : Float_Type'Base;
X : Vector;
INCX : Integer;
Y : Vector;
INCY : Integer
) return Float_Type'Base is
begin
if N > 0 and (
(N - 1) * abs INCX >= X'Length or (N - 1) * abs INCY >= Y'Length
) then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
return SDSDOT (
Fortran_Integer (N),
ALPHA,
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Double =>
return ALPHA + DDOT (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end DSDOT;
function DSDOT (
ALPHA : Float_Type'Base;
X, Y : Vector
) return Float_Type'Base is
begin
if X'Length /= Y'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
return SDSDOT (X'Length, ALPHA, X, 1, Y, 1);
when Double =>
return ALPHA + DDOT (X'Length, X, 1, Y, 1);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end DSDOT;
----------
-- NRM2 --
----------
function SNRM2 (
N : Fortran_Integer;
X : Vector;
INCX : Fortran_Integer
) return Float_Type'Base;
pragma Import (Fortran, SNRM2, Name_Prepend & "snrm2" & Name_Append);
function DNRM2 (
N : Fortran_Integer;
X : Vector;
INCX : Fortran_Integer
) return Float_Type'Base;
pragma Import (Fortran, DNRM2, Name_Prepend & "dnrm2" & Name_Append);
function NRM2 (
N : Natural;
X : Vector;
INCX : Integer
) return Float_Type'Base is
begin
-- Use of INCX rather than abs INCX is deliberate
if N > 0 and (N - 1) * INCX >= X'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
return SNRM2 (Fortran_Integer (N), X, Fortran_Integer (INCX));
when Double =>
return DNRM2 (Fortran_Integer (N), X, Fortran_Integer (INCX));
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end NRM2;
function NRM2 (X : Vector) return Float_Type'Base is
begin
case BLAS_Precision is
when Single =>
return SNRM2 (X'Length, X, 1);
when Double =>
return DNRM2 (X'Length, X, 1);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end NRM2;
---------
-- ROT --
---------
procedure SROT (
N : in Fortran_Integer;
X : in out Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer;
C, S : in Float_Type'Base
);
pragma Import (Fortran, SROT, Name_Prepend & "srot" & Name_Append);
procedure DROT (
N : in Fortran_Integer;
X : in out Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer;
C, S : in Float_Type'Base
);
pragma Import (Fortran, DROT, Name_Prepend & "drot" & Name_Append);
procedure ROT (
N : in Natural;
X : in out Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer;
C, S : in Float_Type'Base
) is
begin
if N > 0 and (
(N - 1) * abs INCX >= X'Length or (N - 1) * abs INCY >= Y'Length
) then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SROT (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY),
C, S
);
when Double =>
DROT (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY),
C, S
);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end ROT;
procedure ROT (
X, Y : in out Vector;
C, S : in Float_Type'Base
) is
begin
if X'Length /= Y'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SROT (X'Length, X, 1, Y, 1, C, S);
when Double =>
DROT (X'Length, X, 1, Y, 1, C, S);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end ROT;
----------
-- ROTM --
----------
procedure SROTM (
N : in Fortran_Integer;
X : in out Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer;
PARAM : in Modified_Givens
);
pragma Import (Fortran, SROTM, Name_Prepend & "srotm" & Name_Append);
procedure DROTM (
N : in Fortran_Integer;
X : in out Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer;
PARAM : in Modified_Givens
);
pragma Import (Fortran, DROTM, Name_Prepend & "drotm" & Name_Append);
procedure ROTM (
N : in Natural;
X : in out Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer;
PARAM : in Modified_Givens
) is
begin
if N > 0 and (
(N - 1) * abs INCX >= X'Length or (N - 1) * abs INCY >= Y'Length
) then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SROTM (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY),
PARAM
);
when Double =>
DROTM (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY),
PARAM
);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end ROTM;
procedure ROTM (
X, Y : in out Vector;
PARAM : in Modified_Givens
) is
begin
if X'Length /= Y'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SROTM (X'Length, X, 1, Y, 1, PARAM);
when Double =>
DROTM (X'Length, X, 1, Y, 1, PARAM);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end ROTM;
----------
-- ROTG --
----------
procedure SROTG (
A, B : in out Float_Type'Base;
C, S : out Float_Type'Base
);
pragma Import (Fortran, SROTG, Name_Prepend & "srotg" & Name_Append);
procedure DROTG (
A, B : in out Float_Type'Base;
C, S : out Float_Type'Base
);
pragma Import (Fortran, DROTG, Name_Prepend & "drotg" & Name_Append);
procedure ROTG (
A, B : in out Float_Type'Base;
C, S : out Float_Type'Base
) is
begin
case BLAS_Precision is
when Single =>
SROTG (A, B, C, S);
when Double =>
DROTG (A, B, C, S);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end ROTG;
-----------
-- ROTMG --
-----------
procedure SROTMG (
D1, D2 : in out Float_Type'Base;
A : in out Float_Type'Base;
B : in Float_Type'Base;
PARAM : out Modified_Givens
);
pragma Import (Fortran, SROTMG, Name_Prepend & "srotmg" & Name_Append);
procedure DROTMG (
D1, D2 : in out Float_Type'Base;
A : in out Float_Type'Base;
B : in Float_Type'Base;
PARAM : out Modified_Givens
);
pragma Import (Fortran, DROTMG, Name_Prepend & "drotmg" & Name_Append);
procedure ROTMG (
D1, D2 : in out Float_Type'Base;
A : in out Float_Type'Base;
B : in Float_Type'Base;
PARAM : out Modified_Givens
) is
begin
case BLAS_Precision is
when Single =>
SROTMG (D1, D2, A, B, PARAM);
when Double =>
DROTMG (D1, D2, A, B, PARAM);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end ROTMG;
----------
-- SCAL --
----------
procedure SSCAL (
N : in Fortran_Integer;
ALPHA : in Float_Type'Base;
X : in out Vector;
INCX : in Fortran_Integer
);
pragma Import (Fortran, SSCAL, Name_Prepend & "sscal" & Name_Append);
procedure DSCAL (
N : in Fortran_Integer;
ALPHA : in Float_Type'Base;
X : in out Vector;
INCX : in Fortran_Integer
);
pragma Import (Fortran, DSCAL, Name_Prepend & "dscal" & Name_Append);
procedure SCAL (
N : in Natural;
ALPHA : in Float_Type'Base;
X : in out Vector;
INCX : in Integer
) is
begin
-- Use of INCX rather than abs INCX is deliberate
if N > 0 and (N - 1) * INCX >= X'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SSCAL (Fortran_Integer (N), ALPHA, X, Fortran_Integer (INCX));
when Double =>
DSCAL (Fortran_Integer (N), ALPHA, X, Fortran_Integer (INCX));
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end SCAL;
procedure SCAL (
ALPHA : in Float_Type'Base;
X : in out Vector
) is
begin
case BLAS_Precision is
when Single =>
SSCAL (X'Length, ALPHA, X, 1);
when Double =>
DSCAL (X'Length, ALPHA, X, 1);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end SCAL;
----------
-- SWAP --
----------
procedure SSWAP (
N : in Fortran_Integer;
X : in out Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer
);
pragma Import (Fortran, SSWAP, Name_Prepend & "sswap" & Name_Append);
procedure DSWAP (
N : in Fortran_Integer;
X : in out Vector;
INCX : in Fortran_Integer;
Y : in out Vector;
INCY : in Fortran_Integer
);
pragma Import (Fortran, DSWAP, Name_Prepend & "dswap" & Name_Append);
procedure SWAP (
N : in Natural;
X : in out Vector;
INCX : in Integer;
Y : in out Vector;
INCY : in Integer
) is
begin
if N > 0 and (
(N - 1) * abs INCX >= X'Length or (N - 1) * abs INCY >= Y'Length
) then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SSWAP (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Double =>
DSWAP (
Fortran_Integer (N),
X,
Fortran_Integer (INCX),
Y,
Fortran_Integer (INCY)
);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end SWAP;
procedure SWAP (X, Y : in out Vector) is
begin
if X'Length /= Y'Length then
raise Argument_Error;
end if;
case BLAS_Precision is
when Single =>
SSWAP (X'Length, X, 1, Y, 1);
when Double =>
DSWAP (X'Length, X, 1, Y, 1);
when Unsupported =>
raise Unsupported_Precision_Error;
end case;
end SWAP;
end Real_BLAS;
^ permalink raw reply [flat|nested] 2+ messages in thread
* Re: BLAS binding
2000-05-11 0:00 BLAS binding Duncan Sands
@ 2000-05-11 0:00 ` Gautier
0 siblings, 0 replies; 2+ messages in thread
From: Gautier @ 2000-05-11 0:00 UTC (permalink / raw)
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;
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2000-05-11 0:00 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2000-05-11 0:00 BLAS binding Duncan Sands
2000-05-11 0:00 ` Gautier
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox