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.9 required=5.0 tests=BAYES_00 autolearn=ham
	autolearn_force=no version=3.4.4
X-Google-Language: ENGLISH,ASCII-7-bit
X-Google-Thread: 103376,95f6e17d8964ca3c,start
X-Google-Attributes: gid103376,public
From: sands@topo.math.u-psud.fr (Duncan Sands)
Subject: BLAS binding
Date: 2000/05/11
Message-ID: <8fdrdf$pag$1@upsn21.u-psud.fr>
X-Deja-AN: 621971457
Organization: Equipe de Topologie et Dynamique (CNRS Orsay)
Newsgroups: comp.lang.ada
Date: 2000-05-11T00:00:00+00:00
List-Id: <comp.lang.ada>

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;