File : ada_blas-real.ads
-------------------------------------------------------------------------------
-- Copyright (C) 2000-2001 Centre National de la Recherche Scientifique --
-- --
-- This program is free software; you can redistribute it and/or modify --
-- it under the terms of the GNU General Public License as published by --
-- the Free Software Foundation; either version 2 of the License, or --
-- (at your option) any later version. --
-- As a special exception, if other files instantiate generics from this --
-- unit, or you link this unit with other files to produce an executable, --
-- this unit does not by itself cause the resulting executable to be --
-- covered by the GNU General Public License. This exception does not --
-- however invalidate any other reasons why the executable file might be --
-- covered by the GNU Public License. --
-- --
-- This program is distributed in the hope that it will be useful, --
-- but WITHOUT ANY WARRANTY; without even the implied warranty of --
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the --
-- GNU General Public License for more details. --
-- --
-- You should have received a copy of the GNU General Public License --
-- along with this program; if not, write to the Free Software --
-- Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA --
-- --
-- Author: Duncan Sands (Duncan.Sands@math.u-psud.fr) --
-- Departement de Mathematiques, Batiment 425, --
-- Universite de Paris-XI, Orsay, France. --
-- http://topo.math.u-psud.fr/~sands --
-------------------------------------------------------------------------------
-- Binding to the real BLAS
-- This package provides a thin binding to both the single and double
-- precision real BLAS.
-- The precision (version) of the Fortran BLAS to be used is determined
-- automatically from the generic parameter Float_Type.
generic
type Float_Type is digits <>;
-- The package determines whether to use the single precision (S) or the
-- double precision (D) Fortran BLAS for subroutine calls (or raise
-- Unsupported_Precision_Error) based on the characteristics of
-- Float_Type'Base
type Index_Type is (<>);
type Vector_Type is array (Index_Type range <>) of Float_Type;
-- Note: the previous line should read
-- type Vector_Type is array (Index_Type range <>) of Float_Type'Base;
-- but this does not compile with version 3.13p of the GNAT compiler.
-- Vector_Type may need to be of convention Fortran:
-- pragma Convention (Fortran, Vector_Type);
type Matrix_Type is array (Index_Type range <>, Index_Type range <>)
of Float_Type;
-- Note: the previous two lines should read
-- type Matrix_Type is array (Index_Type range <>, Index_Type range <>)
-- of Float_Type'Base;
-- but this does not compile with version 3.13p of the GNAT compiler.
-- Matrix_Type MUST be of convention Fortran (column-major order):
-- pragma Convention (Fortran, Matrix_Type);
package Ada_BLAS.Real is
pragma Pure (Ada_BLAS.Real);
Precision : constant Precision_Type;
-- The precision (version) of the Fortran BLAS that will be used
type Modified_Givens is array (1 .. 5) of Float_Type'Base;
pragma Convention (Fortran, Modified_Givens);
--------------
-- Level 1 --
--------------
-- ROTG: Generate plane rotation
procedure ROTG (
A, B : in out Float_Type'Base;
C, S : out Float_Type'Base
);
pragma Inline (ROTG);
-- 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_Type;
INCX : in Integer;
Y : in out Vector_Type;
INCY : in Integer;
C, S : in Float_Type'Base
);
procedure ROT (
X, Y : in out Vector_Type;
C, S : in Float_Type'Base
);
pragma Inline (ROT);
-- ROTM: Apply modified plane rotation
procedure ROTM (
N : in Natural;
X : in out Vector_Type;
INCX : in Integer;
Y : in out Vector_Type;
INCY : in Integer;
PARAM : in Modified_Givens
);
procedure ROTM (
X, Y : in out Vector_Type;
PARAM : in Modified_Givens
);
pragma Inline (ROTM);
-- SWAP: x <-> y
procedure SWAP (
N : in Natural;
X : in out Vector_Type;
INCX : in Integer;
Y : in out Vector_Type;
INCY : in Integer
);
procedure SWAP (X, Y : in out Vector_Type);
pragma Inline (SWAP);
-- SCAL: x <- alpha x
procedure SCAL (
N : in Natural;
ALPHA : in Float_Type'Base;
X : in out Vector_Type;
INCX : in Natural
);
procedure SCAL (
ALPHA : in Float_Type'Base;
X : in out Vector_Type
);
pragma Inline (SCAL);
-- COPY: y <- x
procedure COPY (
N : in Natural;
X : in Vector_Type;
INCX : in Integer;
Y : in out Vector_Type;
INCY : in Integer
);
procedure COPY (
X : in Vector_Type;
Y : in out Vector_Type
);
pragma Inline (COPY);
-- AXPY: y <- alpha x + y
procedure AXPY (
N : in Natural;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
INCX : in Integer;
Y : in out Vector_Type;
INCY : in Integer
);
procedure AXPY (
ALPHA : in Float_Type'Base;
X : in Vector_Type;
Y : in out Vector_Type
);
pragma Inline (AXPY);
-- DOT: dot <- x^T y
function DOT (
N : Natural;
X : Vector_Type;
INCX : Integer;
Y : Vector_Type;
INCY : Integer
) return Float_Type'Base;
function DOT (X, Y : Vector_Type) return Float_Type'Base;
pragma Inline (DOT);
-- DSDOT: dot <- alpha + x^T y (double precision accumulation)
-- Either SDSDOT (single precision) or ALPHA + DDOT (double precision)
function DSDOT (
N : Natural;
ALPHA : Float_Type'Base;
X : Vector_Type;
INCX : Integer;
Y : Vector_Type;
INCY : Integer
) return Float_Type'Base;
function DSDOT (
ALPHA : Float_Type'Base;
X, Y : Vector_Type
) return Float_Type'Base;
pragma Inline (DSDOT);
-- NRM2: nrm2 <- ||x||_2
function NRM2 (
N : Natural;
X : Vector_Type;
INCX : Natural
) return Float_Type'Base;
function NRM2 (X : Vector_Type) return Float_Type'Base;
pragma Inline (NRM2);
-- ASUM: asum <- ||x||_1
function ASUM (
N : Natural;
X : Vector_Type;
INCX : Natural
) return Float_Type'Base;
function ASUM (X : Vector_Type) return Float_Type'Base;
pragma Inline (ASUM);
-- AMAX: amax <- 1st k: |x_k| = max(|x_i|)
-- Either ISAMAX (single precision) or IDAMAX (double precision)
-- with the result converted to the correct value of Index_Type
function AMAX (
N : Natural;
X : Vector_Type;
INCX : Natural
) return Index_Type;
function AMAX (X : Vector_Type) return Index_Type;
pragma Inline (AMAX);
-------------
-- Level 2 --
-------------
-- GEMV: y <- alpha A x + beta y, y <- alpha A^T x + beta y; A - M x N
procedure GEMV (
TRANS : in Transpose_Type;
M : in Natural;
N : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Float_Type'Base;
Y : in out Vector_Type;
INCY : in Integer
);
procedure GEMV (
TRANS : in Transpose_Type;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
X : in Vector_Type;
BETA : in Float_Type'Base;
Y : in out Vector_Type
);
pragma Inline (GEMV);
-- GBMV: y <- alpha A x + beta y, y <- alpha A^T x + beta y; A - M x N
procedure GBMV (
TRANS : in Transpose_Type;
M : in Natural;
N : in Natural;
KL : in Natural;
KU : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Float_Type'Base;
Y : in out Vector_Type;
INCY : in Integer
);
procedure GBMV (
TRANS : in Transpose_Type;
KL : in Natural;
KU : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
X : in Vector_Type;
BETA : in Float_Type'Base;
Y : in out Vector_Type
);
pragma Inline (GBMV);
-- SYMV: y <- alpha A x + beta y
procedure SYMV (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Float_Type'Base;
Y : in out Vector_Type;
INCY : in Integer
);
procedure SYMV (
UPLO : in Triangle_Type;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
X : in Vector_Type;
BETA : in Float_Type'Base;
Y : in out Vector_Type
);
pragma Inline (SYMV);
-- SBMV: y <- alpha A x + beta y
procedure SBMV (
UPLO : in Triangle_Type;
N : in Natural;
K : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Float_Type'Base;
Y : in out Vector_Type;
INCY : in Integer
);
procedure SBMV (
UPLO : in Triangle_Type;
K : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
X : in Vector_Type;
BETA : in Float_Type'Base;
Y : in out Vector_Type
);
pragma Inline (SBMV);
-- SPMV: y <- alpha A x + beta y
procedure SPMV (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Float_Type'Base;
AP : in Vector_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Float_Type'Base;
Y : in out Vector_Type;
INCY : in Integer
);
procedure SPMV (
UPLO : in Triangle_Type;
ALPHA : in Float_Type'Base;
AP : in Vector_Type;
X : in Vector_Type;
BETA : in Float_Type'Base;
Y : in out Vector_Type
);
pragma Inline (SPMV);
-- TRMV: x <- A x, x <- A^T x
procedure TRMV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
N : in Natural;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer
);
procedure TRMV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
A : in Matrix_Type;
X : in Vector_Type
);
pragma Inline (TRMV);
-- TBMV: x <- A x, x <- A^T x
procedure TBMV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
N : in Natural;
K : in Natural;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer
);
procedure TBMV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
K : in Natural;
A : in Matrix_Type;
X : in Vector_Type
);
pragma Inline (TBMV);
-- TPMV: x <- A x, x <- A^T x
procedure TPMV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
N : in Natural;
AP : in Vector_Type;
X : in Vector_Type;
INCX : in Integer
);
procedure TPMV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
AP : in Vector_Type;
X : in Vector_Type
);
pragma Inline (TPMV);
-- TRSV: x <- A^-1 x, x <- A^-T x
procedure TRSV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
N : in Natural;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer
);
procedure TRSV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
A : in Matrix_Type;
X : in Vector_Type
);
pragma Inline (TRSV);
-- TBSV: x <- A^-1 x, x <- A^-T x
procedure TBSV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
N : in Natural;
K : in Natural;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer
);
procedure TBSV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
K : in Natural;
A : in Matrix_Type;
X : in Vector_Type
);
pragma Inline (TBSV);
-- TPSV: x <- A^-1 x, x <- A^-T x
procedure TPSV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
N : in Natural;
AP : in Vector_Type;
X : in Vector_Type;
INCX : in Integer
);
procedure TPSV (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
DIAG : in Diagonal_Type;
AP : in Vector_Type;
X : in Vector_Type
);
pragma Inline (TPSV);
-- GER: A <- alpha x y^T + A; A - M x N
procedure GER (
M : in Natural;
N : in Natural;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
INCX : in Integer;
Y : in Vector_Type;
INCY : in Integer;
A : in out Matrix_Type
);
procedure GER (
ALPHA : in Float_Type'Base;
X : in Vector_Type;
Y : in Vector_Type;
A : in out Matrix_Type
);
pragma Inline (GER);
-- SYR: A <- alpha x x^T + A
procedure SYR (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
INCX : in Integer;
A : in out Matrix_Type
);
procedure SYR (
UPLO : in Triangle_Type;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
A : in out Matrix_Type
);
pragma Inline (SYR);
-- SPR: A <- alpha x x^T + A
procedure SPR (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
INCX : in Integer;
AP : in out Vector_Type
);
procedure SPR (
UPLO : in Triangle_Type;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
AP : in out Vector_Type
);
pragma Inline (SPR);
-- SYR2: A <- alpha x y^T + alpha y x^T + A
procedure SYR2 (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
INCX : in Integer;
Y : in Vector_Type;
INCY : in Integer;
A : in out Matrix_Type
);
procedure SYR2 (
UPLO : in Triangle_Type;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
Y : in Vector_Type;
A : in out Matrix_Type
);
pragma Inline (SYR2);
-- SPR2: A <- alpha x y^T + alpha y x^T + A
procedure SPR2 (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
INCX : in Integer;
Y : in Vector_Type;
INCY : in Integer;
AP : in out Vector_Type
);
procedure SPR2 (
UPLO : in Triangle_Type;
ALPHA : in Float_Type'Base;
X : in Vector_Type;
Y : in Vector_Type;
AP : in out Vector_Type
);
pragma Inline (SPR2);
-------------
-- Level 3 --
-------------
-- GEMM: C <- alpha op(A) op(B) + beta C;
-- op(X) = X, X^T, X^H; C - M x N
procedure GEMM (
TRANSA : in Transpose_Type;
TRANSB : in Transpose_Type;
M : in Natural;
N : in Natural;
K : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Float_Type'Base;
C : in out Matrix_Type
);
procedure GEMM (
TRANSA : in Transpose_Type;
TRANSB : in Transpose_Type;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Float_Type'Base;
C : in out Matrix_Type
);
pragma Inline (GEMM);
-- SYMM: C <- alpha A B + beta C, C <- alpha B A + beta C;
-- C - M x N; A = A^T
procedure SYMM (
SIDE : in Side_Type;
UPLO : in Triangle_Type;
M : in Natural;
N : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Float_Type'Base;
C : in out Matrix_Type
);
procedure SYMM (
SIDE : in Side_Type;
UPLO : in Triangle_Type;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Float_Type'Base;
C : in out Matrix_Type
);
pragma Inline (SYMM);
-- SYRK: C <- alpha A A^T + beta C, C <- alpha A^T A + beta C;
-- C - N x N; C = C^T
procedure SYRK (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
N : in Natural;
K : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
BETA : in Float_Type'Base;
C : in out Matrix_Type
);
procedure SYRK (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
BETA : in Float_Type'Base;
C : in out Matrix_Type
);
pragma Inline (SYRK);
-- SYR2K: C <- alpha A B^T + alpha B A^T + beta C,
-- C <- alpha A^T B + alpha B^T A + beta C;
-- C - N x N; C = C^T
procedure SYR2K (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
N : in Natural;
K : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Float_Type'Base;
C : in out Matrix_Type
);
procedure SYR2K (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Float_Type'Base;
C : in out Matrix_Type
);
pragma Inline (SYR2K);
-- TRMM: B <- alpha op(A) B, B <- alpha B op(A);
-- op(A) = A, A^T, A^H; A triangular; B - M x N
procedure TRMM (
SIDE : in Side_Type;
UPLO : in Triangle_Type;
TRANSA : in Transpose_Type;
DIAG : in Diagonal_Type;
M : in Natural;
N : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type
);
procedure TRMM (
SIDE : in Side_Type;
UPLO : in Triangle_Type;
TRANSA : in Transpose_Type;
DIAG : in Diagonal_Type;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type
);
pragma Inline (TRMM);
-- TRSM: B <- alpha op(A^-1) B, B <- alpha B op(A^-1);
-- op(A) = A, A^T, A^H; A triangular; B - M x N
procedure TRSM (
SIDE : in Side_Type;
UPLO : in Triangle_Type;
TRANSA : in Transpose_Type;
DIAG : in Diagonal_Type;
M : in Natural;
N : in Natural;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type
);
procedure TRSM (
SIDE : in Side_Type;
UPLO : in Triangle_Type;
TRANSA : in Transpose_Type;
DIAG : in Diagonal_Type;
ALPHA : in Float_Type'Base;
A : in Matrix_Type;
B : in Matrix_Type
);
pragma Inline (TRSM);
private
-- The following method of determining Precision was chosen in order to
-- encourage compile time evaluation and thus, hopefully, dead code
-- elimination. Constraint_Error will be raised if double and single
-- precision are the same (this should never happen since the Fortran
-- standard requires double precision to be more precise than single
-- precision). I appreciate that the value of Precision is not
-- guaranteed to be correct for all machines, all compilers and all values
-- of Float_Type. However, I don't know any examples of it being wrong.
-- IF PRECISION IS SET WRONGLY ON YOUR MACHINE, PLEASE LET ME KNOW!
Precision : constant Precision_Type := Precision_Type'Val (
Precision_Type'Pos (Unsupported)
-
Boolean'Pos ( -- Test for double precision
Float_Type'Base'Digits = Interfaces.Fortran.Double_Precision'Digits
and Float_Type'Base'Size = Interfaces.Fortran.Double_Precision'Size
)
- 2 *
Boolean'Pos ( -- Test for single precision
Float_Type'Base'Digits = Interfaces.Fortran.Real'Digits and
Float_Type'Base'Size = Interfaces.Fortran.Real'Size
)
);
end Ada_BLAS.Real;