===================================================================
@@ -6,7 +6,7 @@
-- --
-- B o d y --
-- --
+-- Copyright (C) 2006-2011, Free Software Foundation, Inc. --
-- --
-- GNAT is free software; you can redistribute it and/or modify it under --
-- terms of the GNU General Public License as published by the Free Soft- --
@@ -30,66 +30,35 @@
------------------------------------------------------------------------------
with System.Generic_Array_Operations; use System.Generic_Array_Operations;
-with System.Generic_Complex_BLAS;
-with System.Generic_Complex_LAPACK;
+with Ada.Numerics; use Ada.Numerics;
package body Ada.Numerics.Generic_Complex_Arrays is
- -- Operations involving inner products use BLAS library implementations.
- -- This allows larger matrices and vectors to be computed efficiently,
- -- taking into account memory hierarchy issues and vector instructions
- -- that vary widely between machines.
-
-- Operations that are defined in terms of operations on the type Real,
-- such as addition, subtraction and scaling, are computed in the canonical
-- way looping over all elements.
- -- Operations for solving linear systems and computing determinant,
- -- eigenvalues, eigensystem and inverse, are implemented using the
- -- LAPACK library.
+ package Ops renames System.Generic_Array_Operations;
- type BLAS_Real_Vector is array (Integer range <>) of Real;
-
- package BLAS is new System.Generic_Complex_BLAS
- (Real => Real,
- Complex_Types => Complex_Types,
- Complex_Vector => Complex_Vector,
- Complex_Matrix => Complex_Matrix);
-
- package LAPACK is new System.Generic_Complex_LAPACK
- (Real => Real,
- Real_Vector => BLAS_Real_Vector,
- Complex_Types => Complex_Types,
- Complex_Vector => Complex_Vector,
- Complex_Matrix => Complex_Matrix);
-
subtype Real is Real_Arrays.Real;
-- Work around visibility bug ???
- use BLAS, LAPACK;
+ function Is_Non_Zero (X : Complex) return Boolean is (X /= (0.0, 0.0));
+ -- Needed by Back_Substitute
- -- Procedure versions of functions returning unconstrained values.
- -- This allows for inlining the function wrapper.
+ procedure Back_Substitute is new Ops.Back_Substitute
+ (Scalar => Complex,
+ Matrix => Complex_Matrix,
+ Is_Non_Zero => Is_Non_Zero);
- procedure Eigenvalues
- (A : Complex_Matrix;
- Values : out Real_Vector);
+ procedure Forward_Eliminate is new Ops.Forward_Eliminate
+ (Scalar => Complex,
+ Real => Real'Base,
+ Matrix => Complex_Matrix,
+ Zero => (0.0, 0.0),
+ One => (1.0, 0.0));
- procedure Inverse
- (A : Complex_Matrix;
- R : out Complex_Matrix);
-
- procedure Solve
- (A : Complex_Matrix;
- X : Complex_Vector;
- B : out Complex_Vector);
-
- procedure Solve
- (A : Complex_Matrix;
- X : Complex_Matrix;
- B : out Complex_Matrix);
-
- procedure Transpose is new System.Generic_Array_Operations.Transpose
+ procedure Transpose is new Ops.Transpose
(Scalar => Complex,
Matrix => Complex_Matrix);
@@ -98,6 +67,12 @@
function Length is new Square_Matrix_Length (Complex, Complex_Matrix);
+ -- Instant a generic square root implementation here, in order to avoid
+ -- instantiating a complete copy of Generic_Elementary_Functions.
+ -- Speed of the square root is not a big concern here.
+
+ function Sqrt is new Ops.Sqrt (Real'Base);
+
-- Instantiating the following subprograms directly would lead to
-- name clashes, so use a local package.
@@ -155,6 +130,14 @@
Right_Vector => Complex_Vector,
Zero => (0.0, 0.0));
+ function "*" is new Inner_Product
+ (Left_Scalar => Complex,
+ Right_Scalar => Complex,
+ Result_Scalar => Complex,
+ Left_Vector => Complex_Vector,
+ Right_Vector => Complex_Vector,
+ Zero => (0.0, 0.0));
+
function "*" is new Outer_Product
(Left_Scalar => Complex,
Right_Scalar => Complex,
@@ -229,6 +212,15 @@
Result_Vector => Complex_Vector,
Zero => (0.0, 0.0));
+ function "*" is new Matrix_Vector_Product
+ (Left_Scalar => Complex,
+ Right_Scalar => Complex,
+ Result_Scalar => Complex,
+ Matrix => Complex_Matrix,
+ Right_Vector => Complex_Vector,
+ Result_Vector => Complex_Vector,
+ Zero => (0.0, 0.0));
+
function "*" is new Vector_Matrix_Product
(Left_Scalar => Real'Base,
Right_Scalar => Complex,
@@ -247,7 +239,25 @@
Result_Vector => Complex_Vector,
Zero => (0.0, 0.0));
+ function "*" is new Vector_Matrix_Product
+ (Left_Scalar => Complex,
+ Right_Scalar => Complex,
+ Result_Scalar => Complex,
+ Left_Vector => Complex_Vector,
+ Matrix => Complex_Matrix,
+ Result_Vector => Complex_Vector,
+ Zero => (0.0, 0.0));
+
function "*" is new Matrix_Matrix_Product
+ (Left_Scalar => Complex,
+ Right_Scalar => Complex,
+ Result_Scalar => Complex,
+ Left_Matrix => Complex_Matrix,
+ Right_Matrix => Complex_Matrix,
+ Result_Matrix => Complex_Matrix,
+ Zero => (0.0, 0.0));
+
+ function "*" is new Matrix_Matrix_Product
(Left_Scalar => Real'Base,
Right_Scalar => Complex,
Result_Scalar => Complex,
@@ -445,6 +455,15 @@
Result_Matrix => Complex_Matrix,
Operation => "/");
+ -----------
+ -- "abs" --
+ -----------
+
+ function "abs" is new L2_Norm
+ (X_Scalar => Complex,
+ Result_Real => Real'Base,
+ X_Vector => Complex_Vector);
+
--------------
-- Argument --
--------------
@@ -671,6 +690,16 @@
Y_Matrix => Real_Matrix,
Update => Set_Re);
+ -----------
+ -- Solve --
+ -----------
+
+ function Solve is
+ new Matrix_Vector_Solution (Complex, Complex_Vector, Complex_Matrix);
+
+ function Solve is
+ new Matrix_Matrix_Solution (Complex, Complex_Matrix);
+
-----------------
-- Unit_Matrix --
-----------------
@@ -686,7 +715,6 @@
Vector => Complex_Vector,
Zero => (0.0, 0.0),
One => (1.0, 0.0));
-
end Instantiations;
---------
@@ -696,16 +724,8 @@
function "*"
(Left : Complex_Vector;
Right : Complex_Vector) return Complex
- is
- begin
- if Left'Length /= Right'Length then
- raise Constraint_Error with
- "vectors are of different length in inner product";
- end if;
+ renames Instantiations."*";
- return dot (Left'Length, X => Left, Y => Right);
- end "*";
-
function "*"
(Left : Real_Vector;
Right : Complex_Vector) return Complex
@@ -738,32 +758,9 @@
function "*"
(Left : Complex_Matrix;
- Right : Complex_Matrix)
- return Complex_Matrix
- is
- R : Complex_Matrix (Left'Range (1), Right'Range (2));
+ Right : Complex_Matrix) return Complex_Matrix
+ renames Instantiations."*";
- begin
- if Left'Length (2) /= Right'Length (1) then
- raise Constraint_Error with
- "incompatible dimensions in matrix-matrix multiplication";
- end if;
-
- gemm (Trans_A => No_Trans'Access,
- Trans_B => No_Trans'Access,
- M => Right'Length (2),
- N => Left'Length (1),
- K => Right'Length (1),
- A => Right,
- Ld_A => Right'Length (2),
- B => Left,
- Ld_B => Left'Length (2),
- C => R,
- Ld_C => R'Length (2));
-
- return R;
- end "*";
-
function "*"
(Left : Complex_Vector;
Right : Complex_Vector) return Complex_Matrix
@@ -772,49 +769,13 @@
function "*"
(Left : Complex_Vector;
Right : Complex_Matrix) return Complex_Vector
- is
- R : Complex_Vector (Right'Range (2));
+ renames Instantiations."*";
- begin
- if Left'Length /= Right'Length (1) then
- raise Constraint_Error with
- "incompatible dimensions in vector-matrix multiplication";
- end if;
-
- gemv (Trans => No_Trans'Access,
- M => Right'Length (2),
- N => Right'Length (1),
- A => Right,
- Ld_A => Right'Length (2),
- X => Left,
- Y => R);
-
- return R;
- end "*";
-
function "*"
(Left : Complex_Matrix;
Right : Complex_Vector) return Complex_Vector
- is
- R : Complex_Vector (Left'Range (1));
+ renames Instantiations."*";
- begin
- if Left'Length (2) /= Right'Length then
- raise Constraint_Error with
- "incompatible dimensions in matrix-vector multiplication";
- end if;
-
- gemv (Trans => Trans'Access,
- M => Left'Length (2),
- N => Left'Length (1),
- A => Left,
- Ld_A => Left'Length (2),
- X => Right,
- Y => R);
-
- return R;
- end "*";
-
function "*"
(Left : Real_Matrix;
Right : Complex_Matrix) return Complex_Matrix
@@ -984,10 +945,8 @@
-- "abs" --
-----------
- function "abs" (Right : Complex_Vector) return Complex is
- begin
- return (nrm2 (Right'Length, Right), 0.0);
- end "abs";
+ function "abs" (Right : Complex_Vector) return Real'Base
+ renames Instantiations."abs";
--------------
-- Argument --
@@ -1070,38 +1029,12 @@
-----------------
function Determinant (A : Complex_Matrix) return Complex is
- N : constant Integer := Length (A);
- LU : Complex_Matrix (1 .. N, 1 .. N) := A;
- Piv : Integer_Vector (1 .. N);
- Info : aliased Integer := -1;
- Neg : Boolean;
- Det : Complex;
-
+ M : Complex_Matrix := A;
+ B : Complex_Matrix (A'Range (1), 1 .. 0);
+ R : Complex;
begin
- if N = 0 then
- return (0.0, 0.0);
- end if;
-
- getrf (N, N, LU, N, Piv, Info'Access);
-
- if Info /= 0 then
- raise Constraint_Error with "ill-conditioned matrix";
- end if;
-
- Det := LU (1, 1);
- Neg := Piv (1) /= 1;
-
- for J in 2 .. N loop
- Det := Det * LU (J, J);
- Neg := Neg xor (Piv (J) /= J);
- end loop;
-
- if Neg then
- return -Det;
-
- else
- return Det;
- end if;
+ Forward_Eliminate (M, B, R);
+ return R;
end Determinant;
-----------------
@@ -1113,174 +1046,96 @@
Values : out Real_Vector;
Vectors : out Complex_Matrix)
is
- Job_Z : aliased Character := 'V';
- Rng : aliased Character := 'A';
- Uplo : aliased Character := 'U';
+ N : constant Natural := Length (A);
- N : constant Natural := Length (A);
- W : BLAS_Real_Vector (Values'Range);
- M : Integer;
- B : Complex_Matrix (1 .. N, 1 .. N);
- L_Work : Complex_Vector (1 .. 1);
- LR_Work : BLAS_Real_Vector (1 .. 1);
- LI_Work : Integer_Vector (1 .. 1);
- I_Supp_Z : Integer_Vector (1 .. 2 * N);
- Info : aliased Integer;
+ -- For a Hermitian matrix C, we convert the eigenvalue problem to a
+ -- real symmetric one: if C = A + i * B, then the (N, N) complex
+ -- eigenvalue problem:
+ -- (A + i * B) * (u + i * v) = Lambda * (u + i * v)
+ --
+ -- is equivalent to the (2 * N, 2 * N) real eigenvalue problem:
+ -- [ A, B ] [ u ] = Lambda * [ u ]
+ -- [ -B, A ] [ v ] [ v ]
+ --
+ -- Note that the (2 * N, 2 * N) matrix above is symmetric, as
+ -- Transpose (A) = A and Transpose (B) = -B if C is Hermitian.
- begin
- if Values'Length /= N then
- raise Constraint_Error with "wrong length for output vector";
- end if;
+ -- We solve this eigensystem using the real-valued algorithms. The final
+ -- result will have every eigenvalue twice, so in the sorted output we
+ -- just pick every second value, with associated eigenvector u + i * v.
- if Vectors'First (1) /= A'First (1)
- or else Vectors'Last (1) /= A'Last (1)
- or else Vectors'First (2) /= A'First (2)
- or else Vectors'Last (2) /= A'Last (2)
- then
- raise Constraint_Error with "wrong dimensions for output matrix";
- end if;
+ M : Real_Matrix (1 .. 2 * N, 1 .. 2 * N);
+ Vals : Real_Vector (1 .. 2 * N);
+ Vecs : Real_Matrix (1 .. 2 * N, 1 .. 2 * N);
- if N = 0 then
- return;
- end if;
+ begin
+ for J in 1 .. N loop
+ for K in 1 .. N loop
+ declare
+ C : constant Complex :=
+ (A (A'First (1) + (J - 1), A'First (2) + (K - 1)));
+ begin
+ M (J, K) := Re (C);
+ M (J + N, K + N) := Re (C);
+ M (J + N, K) := Im (C);
+ M (J, K + N) := -Im (C);
+ end;
+ end loop;
+ end loop;
- -- Check for hermitian matrix ???
- -- Copy only required triangle ???
+ Eigensystem (M, Vals, Vecs);
- B := A;
+ for J in 1 .. N loop
+ declare
+ Col : constant Integer := Values'First + (J - 1);
+ begin
+ Values (Col) := Vals (2 * J);
- -- Find size of work area
-
- heevr
- (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
- M => M,
- W => W,
- Z => Vectors,
- Ld_Z => N,
- I_Supp_Z => I_Supp_Z,
- Work => L_Work,
- L_Work => -1,
- R_Work => LR_Work,
- LR_Work => -1,
- I_Work => LI_Work,
- LI_Work => -1,
- Info => Info'Access);
-
- if Info /= 0 then
- raise Constraint_Error;
- end if;
-
- declare
- Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
- R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
- I_Work : Integer_Vector (1 .. LI_Work (1));
-
- begin
- heevr
- (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
- M => M,
- W => W,
- Z => Vectors,
- Ld_Z => N,
- I_Supp_Z => I_Supp_Z,
- Work => Work,
- L_Work => Work'Length,
- R_Work => R_Work,
- LR_Work => LR_Work'Length,
- I_Work => I_Work,
- LI_Work => LI_Work'Length,
- Info => Info'Access);
-
- if Info /= 0 then
- raise Constraint_Error with "inverting non-Hermitian matrix";
- end if;
-
- for J in Values'Range loop
- Values (J) := W (J);
- end loop;
- end;
+ for K in 1 .. N loop
+ declare
+ Row : constant Integer := Vectors'First (2) + (K - 1);
+ begin
+ Vectors (Row, Col)
+ := (Vecs (J * 2, Col), Vecs (J * 2, Col + N));
+ end;
+ end loop;
+ end;
+ end loop;
end Eigensystem;
-----------------
-- Eigenvalues --
-----------------
- procedure Eigenvalues
- (A : Complex_Matrix;
- Values : out Real_Vector)
- is
- Job_Z : aliased Character := 'N';
- Rng : aliased Character := 'A';
- Uplo : aliased Character := 'U';
- N : constant Natural := Length (A);
- B : Complex_Matrix (1 .. N, 1 .. N) := A;
- Z : Complex_Matrix (1 .. 1, 1 .. 1);
- W : BLAS_Real_Vector (Values'Range);
- L_Work : Complex_Vector (1 .. 1);
- LR_Work : BLAS_Real_Vector (1 .. 1);
- LI_Work : Integer_Vector (1 .. 1);
- I_Supp_Z : Integer_Vector (1 .. 2 * N);
- M : Integer;
- Info : aliased Integer;
+ function Eigenvalues (A : Complex_Matrix) return Real_Vector is
+ -- See Eigensystem for a description of the algorithm
+ N : constant Natural := Length (A);
+ R : Real_Vector (A'Range (1));
+
+ M : Real_Matrix (1 .. 2 * N, 1 .. 2 * N);
+ Vals : Real_Vector (1 .. 2 * N);
begin
- if Values'Length /= N then
- raise Constraint_Error with "wrong length for output vector";
- end if;
+ for J in 1 .. N loop
+ for K in 1 .. N loop
+ declare
+ C : constant Complex :=
+ (A (A'First (1) + (J - 1), A'First (2) + (K - 1)));
+ begin
+ M (J, K) := Re (C);
+ M (J + N, K + N) := Re (C);
+ M (J + N, K) := Im (C);
+ M (J, K + N) := -Im (C);
+ end;
+ end loop;
+ end loop;
- if N = 0 then
- return;
- end if;
+ Vals := Eigenvalues (M);
- -- Check for hermitian matrix ???
+ for J in 1 .. N loop
+ R (A'First (1) + (J - 1)) := Vals (2 * J);
+ end loop;
- -- Find size of work area
-
- heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
- M => M,
- W => W,
- Z => Z,
- Ld_Z => 1,
- I_Supp_Z => I_Supp_Z,
- Work => L_Work, L_Work => -1,
- R_Work => LR_Work, LR_Work => -1,
- I_Work => LI_Work, LI_Work => -1,
- Info => Info'Access);
-
- if Info /= 0 then
- raise Constraint_Error;
- end if;
-
- declare
- Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
- R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
- I_Work : Integer_Vector (1 .. LI_Work (1));
- begin
- heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
- M => M,
- W => W,
- Z => Z,
- Ld_Z => 1,
- I_Supp_Z => I_Supp_Z,
- Work => Work, L_Work => Work'Length,
- R_Work => R_Work, LR_Work => R_Work'Length,
- I_Work => I_Work, LI_Work => I_Work'Length,
- Info => Info'Access);
-
- if Info /= 0 then
- raise Constraint_Error with "inverting singular matrix";
- end if;
-
- for J in Values'Range loop
- Values (J) := W (J);
- end loop;
- end;
- end Eigenvalues;
-
- function Eigenvalues (A : Complex_Matrix) return Real_Vector is
- R : Real_Vector (A'Range (1));
- begin
- Eigenvalues (A, R);
return R;
end Eigenvalues;
@@ -1298,73 +1153,8 @@
-- Inverse --
-------------
- procedure Inverse (A : Complex_Matrix; R : out Complex_Matrix) is
- N : constant Integer := Length (A);
- Piv : Integer_Vector (1 .. N);
- L_Work : Complex_Vector (1 .. 1);
- Info : aliased Integer := -1;
-
- begin
- -- All computations are done using column-major order, but this works
- -- out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
-
- R := A;
-
- -- Compute LU decomposition
-
- getrf (M => N,
- N => N,
- A => R,
- Ld_A => N,
- I_Piv => Piv,
- Info => Info'Access);
-
- if Info /= 0 then
- raise Constraint_Error with "inverting singular matrix";
- end if;
-
- -- Determine size of work area
-
- getri (N => N,
- A => R,
- Ld_A => N,
- I_Piv => Piv,
- Work => L_Work,
- L_Work => -1,
- Info => Info'Access);
-
- if Info /= 0 then
- raise Constraint_Error;
- end if;
-
- declare
- Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
-
- begin
- -- Compute inverse from LU decomposition
-
- getri (N => N,
- A => R,
- Ld_A => N,
- I_Piv => Piv,
- Work => Work,
- L_Work => Work'Length,
- Info => Info'Access);
-
- if Info /= 0 then
- raise Constraint_Error with "inverting singular matrix";
- end if;
-
- -- ??? Should iterate with gerfs, based on implementation advice
- end;
- end Inverse;
-
function Inverse (A : Complex_Matrix) return Complex_Matrix is
- R : Complex_Matrix (A'Range (2), A'Range (1));
- begin
- Inverse (A, R);
- return R;
- end Inverse;
+ (Solve (A, Unit_Matrix (Length (A))));
-------------
-- Modulus --
@@ -1418,54 +1208,16 @@
-- Solve --
-----------
- procedure Solve
+ function Solve
(A : Complex_Matrix;
- X : Complex_Vector;
- B : out Complex_Vector)
- is
- begin
- if Length (A) /= X'Length then
- raise Constraint_Error with
- "incompatible matrix and vector dimensions";
- end if;
+ X : Complex_Vector) return Complex_Vector
+ renames Instantiations.Solve;
- -- ??? Should solve directly, is faster and more accurate
-
- B := Inverse (A) * X;
- end Solve;
-
- procedure Solve
- (A : Complex_Matrix;
- X : Complex_Matrix;
- B : out Complex_Matrix)
- is
- begin
- if Length (A) /= X'Length (1) then
- raise Constraint_Error with "incompatible matrix dimensions";
- end if;
-
- -- ??? Should solve directly, is faster and more accurate
-
- B := Inverse (A) * X;
- end Solve;
-
function Solve
(A : Complex_Matrix;
- X : Complex_Vector) return Complex_Vector
- is
- B : Complex_Vector (A'Range (2));
- begin
- Solve (A, X, B);
- return B;
- end Solve;
+ X : Complex_Matrix) return Complex_Matrix
+ renames Instantiations.Solve;
- function Solve (A, X : Complex_Matrix) return Complex_Matrix is
- B : Complex_Matrix (A'Range (2), X'Range (2));
- begin
- Solve (A, X, B);
- return B;
- end Solve;
-
---------------
-- Transpose --
---------------
===================================================================
@@ -66,7 +66,7 @@
function "+" (Left, Right : Complex_Vector) return Complex_Vector;
function "-" (Left, Right : Complex_Vector) return Complex_Vector;
function "*" (Left, Right : Complex_Vector) return Complex;
- function "abs" (Right : Complex_Vector) return Complex;
+ function "abs" (Right : Complex_Vector) return Real'Base;
-- Mixed Real_Vector and Complex_Vector arithmetic operations
===================================================================
@@ -33,7 +33,7 @@
-- reason for this is new Ada 2012 requirements that prohibit algorithms such
-- as Strassen's algorithm, which may be used by some BLAS implementations. In
-- addition, some platforms lacked suitable compilers to compile the reference
+-- BLAS/LAPACK implementation. Finally, on some platforms there are more
-- floating point types than supported by BLAS/LAPACK.
with Ada.Containers.Generic_Anonymous_Array_Sort; use Ada.Containers;
@@ -337,6 +337,11 @@
Result_Matrix => Real_Matrix,
Operation => "abs");
+ function Solve is
+ new Matrix_Vector_Solution (Real'Base, Real_Vector, Real_Matrix);
+
+ function Solve is new Matrix_Matrix_Solution (Real'Base, Real_Matrix);
+
function Unit_Matrix is new
Generic_Array_Operations.Unit_Matrix
(Scalar => Real'Base,
@@ -696,59 +701,12 @@
-- Solve --
-----------
- function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector is
- N : constant Natural := Length (A);
- MA : Real_Matrix := A;
- MX : Real_Matrix (A'Range (1), 1 .. 1);
- R : Real_Vector (A'Range (2));
- Det : Real'Base;
+ function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector
+ renames Instantiations.Solve;
- begin
- if X'Length /= N then
- raise Constraint_Error with "incompatible vector length";
- end if;
+ function Solve (A, X : Real_Matrix) return Real_Matrix
+ renames Instantiations.Solve;
- for J in 0 .. MX'Length (1) - 1 loop
- MX (MX'First (1) + J, 1) := X (X'First + J);
- end loop;
-
- Forward_Eliminate (MA, MX, Det);
- Back_Substitute (MA, MX);
-
- for J in 0 .. R'Length - 1 loop
- R (R'First + J) := MX (MX'First (1) + J, 1);
- end loop;
-
- return R;
- end Solve;
-
- function Solve (A, X : Real_Matrix) return Real_Matrix is
- N : constant Natural := Length (A);
- MA : Real_Matrix (A'Range (2), A'Range (2));
- MB : Real_Matrix (A'Range (2), X'Range (2));
- Det : Real'Base;
-
- begin
- if X'Length (1) /= N then
- raise Constraint_Error with "matrices have unequal number of rows";
- end if;
-
- for J in 0 .. A'Length (1) - 1 loop
- for K in MA'Range (2) loop
- MA (MA'First (1) + J, K) := A (A'First (1) + J, K);
- end loop;
-
- for K in MB'Range (2) loop
- MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
- end loop;
- end loop;
-
- Forward_Eliminate (MA, MB, Det);
- Back_Substitute (MA, MB);
-
- return MB;
- end Solve;
-
----------------------
-- Sort_Eigensystem --
----------------------
===================================================================
@@ -651,6 +651,75 @@
return R;
end Matrix_Matrix_Product;
+ ----------------------------
+ -- Matrix_Vector_Solution --
+ ----------------------------
+
+ function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is
+ N : constant Natural := A'Length (1);
+ MA : Matrix := A;
+ MX : Matrix (A'Range (1), 1 .. 1);
+ R : Vector (A'Range (2));
+ Det : Scalar;
+
+ begin
+ if A'Length (2) /= N then
+ raise Constraint_Error with "matrix is not square";
+ end if;
+
+ if X'Length /= N then
+ raise Constraint_Error with "incompatible vector length";
+ end if;
+
+ for J in 0 .. MX'Length (1) - 1 loop
+ MX (MX'First (1) + J, 1) := X (X'First + J);
+ end loop;
+
+ Forward_Eliminate (MA, MX, Det);
+ Back_Substitute (MA, MX);
+
+ for J in 0 .. R'Length - 1 loop
+ R (R'First + J) := MX (MX'First (1) + J, 1);
+ end loop;
+
+ return R;
+ end Matrix_Vector_Solution;
+
+ ----------------------------
+ -- Matrix_Matrix_Solution --
+ ----------------------------
+
+ function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is
+ N : constant Natural := A'Length (1);
+ MA : Matrix (A'Range (2), A'Range (2));
+ MB : Matrix (A'Range (2), X'Range (2));
+ Det : Scalar;
+
+ begin
+ if A'Length (2) /= N then
+ raise Constraint_Error with "matrix is not square";
+ end if;
+
+ if X'Length (1) /= N then
+ raise Constraint_Error with "matrices have unequal number of rows";
+ end if;
+
+ for J in 0 .. A'Length (1) - 1 loop
+ for K in MA'Range (2) loop
+ MA (MA'First (1) + J, K) := A (A'First (1) + J, K);
+ end loop;
+
+ for K in MB'Range (2) loop
+ MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
+ end loop;
+ end loop;
+
+ Forward_Eliminate (MA, MB, Det);
+ Back_Substitute (MA, MB);
+
+ return MB;
+ end Matrix_Matrix_Solution;
+
---------------------------
-- Matrix_Vector_Product --
---------------------------
===================================================================
@@ -390,6 +390,35 @@
(Left : Left_Matrix;
Right : Right_Matrix) return Result_Matrix;
+ ----------------------------
+ -- Matrix_Vector_Solution --
+ ----------------------------
+
+ generic
+ type Scalar is private;
+ type Vector is array (Integer range <>) of Scalar;
+ type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+ with procedure Back_Substitute (M, N : in out Matrix) is <>;
+ with procedure Forward_Eliminate
+ (M : in out Matrix;
+ N : in out Matrix;
+ Det : out Scalar) is <>;
+ function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector;
+
+ ----------------------------
+ -- Matrix_Matrix_Solution --
+ ----------------------------
+
+ generic
+ type Scalar is private;
+ type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+ with procedure Back_Substitute (M, N : in out Matrix) is <>;
+ with procedure Forward_Eliminate
+ (M : in out Matrix;
+ N : in out Matrix;
+ Det : out Scalar) is <>;
+ function Matrix_Matrix_Solution (A : Matrix; X : Matrix) return Matrix;
+
----------
-- Sqrt --
----------
This completes the removal of dependencies on BLAS and LAPACK for these packages. The main reason for this is limited availability of these libraries on some platforms and for some types, in particular types wider than 64 bits. Furthermore, some BLAS implementations may use sub-cubic implementations for matrix multiplication, which would violate Ada 2012 accuracy requirements. Mostly the complex versions of the various routines could use generalized versions of the existing implementation of Generic_Real_Arrays. For solving eigensystems on Hermitian matrices, the implementation for symmetric real matrices is reused, through the use of an augmented matrix with doubled dimensions. This change also addresses AI05-0047, which fixes return type of the "abs" function to be real instead of complex. 2011-10-13 Geert Bosch <bosch@adacore.com> * a-ngrear.adb (Solve): Make generic and move to System.Generic_Array_Operations. * s-gearop.ads (Matrix_Vector_Solution, Matrix_Matrix_Solution): New generic solvers to compute a vector resp. matrix Y such that A * Y = X, approximately. * s-gearop.adb (Matrix_Vector_Solution, Matrix_Matrix_Solution): Implement using Forward_Eliminate and Back_Substitute * a-ngcoar.adb: Reimplement in pure Ada to remove dependencies on BLAS and LAPACK. * a-ngcoar.ads ("abs"): Fix return type to be real.