File : ada_blas-complex.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 complex BLAS
-- This package provides a thin binding to both the single and double
-- precision complex BLAS.
-- The precision (version) of the Fortran BLAS to be used is determined
-- automatically from the generic parameter Float_Type.
with Ada.Numerics.Generic_Complex_Types;
generic
type Float_Type is digits <>;
-- The package determines whether to use the single precision (C) or the
-- double precision (Z) Fortran BLAS for subroutine calls (or raise
-- Unsupported_Precision_Error) based on the characteristics of
-- Float_Type'Base
with package Complex_Types is
new Ada.Numerics.Generic_Complex_Types (Float_Type);
type Complex_Type is new Complex_Types.Complex;
type Index_Type is (<>);
type Vector_Type is array (Index_Type range <>) of Complex_Type;
-- 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 Complex_Type;
-- Matrix_Type MUST be of convention Fortran (column-major order):
-- pragma Convention (Fortran, Matrix_Type);
package Ada_BLAS.Complex is
pragma Pure (Ada_BLAS.Complex);
Precision : constant Precision_Type;
-- The precision (version) of the Fortran BLAS that will be used
--------------
-- Level 1 --
--------------
-- 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 Complex_Type;
X : in out Vector_Type;
INCX : in Natural
);
procedure SCAL (
N : in Natural;
ALPHA : in Float_Type'Base;
X : in out Vector_Type;
INCX : in Natural
);
procedure SCAL (
ALPHA : in Complex_Type;
X : in out Vector_Type
);
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 Complex_Type;
X : in Vector_Type;
INCX : in Integer;
Y : in out Vector_Type;
INCY : in Integer
);
procedure AXPY (
ALPHA : in Complex_Type;
X : in Vector_Type;
Y : in out Vector_Type
);
pragma Inline (AXPY);
-- DOTU: dot <- x^T y
function DOTU (
N : Natural;
X : Vector_Type;
INCX : Integer;
Y : Vector_Type;
INCY : Integer
) return Complex_Type;
function DOTU (X, Y : Vector_Type) return Complex_Type;
pragma Inline (DOTU);
-- DOTC: dot <- x^H y
function DOTC (
N : Natural;
X : Vector_Type;
INCX : Integer;
Y : Vector_Type;
INCY : Integer
) return Complex_Type;
function DOTC (X, Y : Vector_Type) return Complex_Type;
pragma Inline (DOTC);
-- 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 <- ||re(x)||_1 + ||im(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: |re(x_k)| + |im(x_k)| = max(|re(x_i)| + |im(x_i)|)
-- Either ICAMAX (single precision) or IZAMAX (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,
-- y <- alpha A^H x + beta y; A - M x N
procedure GEMV (
TRANS : in Transpose_Type;
M : in Natural;
N : in Natural;
ALPHA : in Complex_Type;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Complex_Type;
Y : in out Vector_Type;
INCY : in Integer
);
procedure GEMV (
TRANS : in Transpose_Type;
ALPHA : in Complex_Type;
A : in Matrix_Type;
X : in Vector_Type;
BETA : in Complex_Type;
Y : in out Vector_Type
);
pragma Inline (GEMV);
-- GBMV: y <- alpha A x + beta y, y <- alpha A^T x + beta y,
-- y <- alpha A^H 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 Complex_Type;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Complex_Type;
Y : in out Vector_Type;
INCY : in Integer
);
procedure GBMV (
TRANS : in Transpose_Type;
KL : in Natural;
KU : in Natural;
ALPHA : in Complex_Type;
A : in Matrix_Type;
X : in Vector_Type;
BETA : in Complex_Type;
Y : in out Vector_Type
);
pragma Inline (GBMV);
-- HEMV: y <- alpha A x + beta y
procedure HEMV (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Complex_Type;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Complex_Type;
Y : in out Vector_Type;
INCY : in Integer
);
procedure HEMV (
UPLO : in Triangle_Type;
ALPHA : in Complex_Type;
A : in Matrix_Type;
X : in Vector_Type;
BETA : in Complex_Type;
Y : in out Vector_Type
);
pragma Inline (HEMV);
-- HBMV: y <- alpha A x + beta y
procedure HBMV (
UPLO : in Triangle_Type;
N : in Natural;
K : in Natural;
ALPHA : in Complex_Type;
A : in Matrix_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Complex_Type;
Y : in out Vector_Type;
INCY : in Integer
);
procedure HBMV (
UPLO : in Triangle_Type;
K : in Natural;
ALPHA : in Complex_Type;
A : in Matrix_Type;
X : in Vector_Type;
BETA : in Complex_Type;
Y : in out Vector_Type
);
pragma Inline (HBMV);
-- HPMV: y <- alpha A x + beta y
procedure HPMV (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Complex_Type;
AP : in Vector_Type;
X : in Vector_Type;
INCX : in Integer;
BETA : in Complex_Type;
Y : in out Vector_Type;
INCY : in Integer
);
procedure HPMV (
UPLO : in Triangle_Type;
ALPHA : in Complex_Type;
AP : in Vector_Type;
X : in Vector_Type;
BETA : in Complex_Type;
Y : in out Vector_Type
);
pragma Inline (HPMV);
-- TRMV: x <- A x, x <- A^T x, x <- A^H 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, x <- A^H 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, x <- A^H 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, x <- A^-H 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, x <- A^-H 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, x <- A^-H 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);
-- GERU: A <- alpha x y^T + A; A - M x N
procedure GERU (
M : in Natural;
N : in Natural;
ALPHA : in Complex_Type;
X : in Vector_Type;
INCX : in Integer;
Y : in Vector_Type;
INCY : in Integer;
A : in out Matrix_Type
);
procedure GERU (
ALPHA : in Complex_Type;
X : in Vector_Type;
Y : in Vector_Type;
A : in out Matrix_Type
);
pragma Inline (GERU);
-- GERC: A <- alpha x y^H + A; A - M x N
procedure GERC (
M : in Natural;
N : in Natural;
ALPHA : in Complex_Type;
X : in Vector_Type;
INCX : in Integer;
Y : in Vector_Type;
INCY : in Integer;
A : in out Matrix_Type
);
procedure GERC (
ALPHA : in Complex_Type;
X : in Vector_Type;
Y : in Vector_Type;
A : in out Matrix_Type
);
pragma Inline (GERC);
-- HER: A <- alpha x x^H + A
procedure HER (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Complex_Type;
X : in Vector_Type;
INCX : in Integer;
A : in out Matrix_Type
);
procedure HER (
UPLO : in Triangle_Type;
ALPHA : in Complex_Type;
X : in Vector_Type;
A : in out Matrix_Type
);
pragma Inline (HER);
-- HPR: A <- alpha x x^H + A
procedure HPR (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Complex_Type;
X : in Vector_Type;
INCX : in Integer;
AP : in out Vector_Type
);
procedure HPR (
UPLO : in Triangle_Type;
ALPHA : in Complex_Type;
X : in Vector_Type;
AP : in out Vector_Type
);
pragma Inline (HPR);
-- HER2: A <- alpha x y^H + y (alpha x)^H + A
procedure HER2 (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Complex_Type;
X : in Vector_Type;
INCX : in Integer;
Y : in Vector_Type;
INCY : in Integer;
A : in out Matrix_Type
);
procedure HER2 (
UPLO : in Triangle_Type;
ALPHA : in Complex_Type;
X : in Vector_Type;
Y : in Vector_Type;
A : in out Matrix_Type
);
pragma Inline (HER2);
-- HPR2: A <- alpha x y^H + y (alpha x)^H + A
procedure HPR2 (
UPLO : in Triangle_Type;
N : in Natural;
ALPHA : in Complex_Type;
X : in Vector_Type;
INCX : in Integer;
Y : in Vector_Type;
INCY : in Integer;
AP : in out Vector_Type
);
procedure HPR2 (
UPLO : in Triangle_Type;
ALPHA : in Complex_Type;
X : in Vector_Type;
Y : in Vector_Type;
AP : in out Vector_Type
);
pragma Inline (HPR2);
-------------
-- 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 Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
procedure GEMM (
TRANSA : in Transpose_Type;
TRANSB : in Transpose_Type;
ALPHA : in Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
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 Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
procedure SYMM (
SIDE : in Side_Type;
UPLO : in Triangle_Type;
ALPHA : in Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
pragma Inline (SYMM);
-- HEMM: C <- alpha A B + beta C, C <- alpha B A + beta C;
-- C - M x N; A = A^H
procedure HEMM (
SIDE : in Side_Type;
UPLO : in Triangle_Type;
M : in Natural;
N : in Natural;
ALPHA : in Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
procedure HEMM (
SIDE : in Side_Type;
UPLO : in Triangle_Type;
ALPHA : in Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
pragma Inline (HEMM);
-- 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 Complex_Type;
A : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
procedure SYRK (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
ALPHA : in Complex_Type;
A : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
pragma Inline (SYRK);
-- HERK: C <- alpha A A^H + beta C, C <- alpha A^H A + beta C;
-- C - N x N; C = C^H
procedure HERK (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
N : in Natural;
K : in Natural;
ALPHA : in Complex_Type;
A : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
procedure HERK (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
ALPHA : in Complex_Type;
A : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
pragma Inline (HERK);
-- 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 Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
procedure SYR2K (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
ALPHA : in Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
pragma Inline (SYR2K);
-- HER2K: 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 HER2K (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
N : in Natural;
K : in Natural;
ALPHA : in Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
procedure HER2K (
UPLO : in Triangle_Type;
TRANS : in Transpose_Type;
ALPHA : in Complex_Type;
A : in Matrix_Type;
B : in Matrix_Type;
BETA : in Complex_Type;
C : in out Matrix_Type
);
pragma Inline (HER2K);
-- 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 Complex_Type;
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 Complex_Type;
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 Complex_Type;
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 Complex_Type;
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.Complex;