Skip to content

R
riscv-gcc-1
Commit a5917ffb by Geert Bosch Committed by Arnaud Charlet
Options

a-ngrear.adb (Solve): Make generic and move to System.Generic_Array_Operations. 

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.

From-SVN: r179912
parent 574ec945
Showing with 293 additions and 472 deletions
gcc/ada/ChangeLog
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.
2011-10-13 Eric Botcazou <ebotcazou@adacore.com> 2011-10-13 Eric Botcazou <ebotcazou@adacore.com>
PR ada/50589 PR ada/50589
......
gcc/ada/a-ngcoar.adb
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
-- -- -- --
-- B o d y -- -- B o d y --
-- -- -- --
-- Copyright (C) 2006-2009, Free Software Foundation, Inc. -- -- Copyright (C) 2006-2011, Free Software Foundation, Inc. --
-- -- -- --
-- GNAT is free software; you can redistribute it and/or modify it under -- -- 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- -- -- terms of the GNU General Public License as published by the Free Soft- --
...@@ -30,66 +30,35 @@ ...@@ -30,66 +30,35 @@
------------------------------------------------------------------------------ ------------------------------------------------------------------------------
with System.Generic_Array_Operations; use System.Generic_Array_Operations; with System.Generic_Array_Operations; use System.Generic_Array_Operations;
with System.Generic_Complex_BLAS; with Ada.Numerics; use Ada.Numerics;
with System.Generic_Complex_LAPACK;
package body Ada.Numerics.Generic_Complex_Arrays is 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, -- Operations that are defined in terms of operations on the type Real,
-- such as addition, subtraction and scaling, are computed in the canonical -- such as addition, subtraction and scaling, are computed in the canonical
-- way looping over all elements. -- way looping over all elements.
-- Operations for solving linear systems and computing determinant, package Ops renames System.Generic_Array_Operations;
-- eigenvalues, eigensystem and inverse, are implemented using the
-- LAPACK library.
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; subtype Real is Real_Arrays.Real;
-- Work around visibility bug ??? -- 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 Eigenvalues procedure Back_Substitute is new Ops.Back_Substitute
(A : Complex_Matrix; (Scalar => Complex,
Values : out Real_Vector); Matrix => Complex_Matrix,
Is_Non_Zero => Is_Non_Zero);
procedure Inverse procedure Forward_Eliminate is new Ops.Forward_Eliminate
(A : Complex_Matrix; (Scalar => Complex,
R : out Complex_Matrix); Real => Real'Base,
Matrix => Complex_Matrix,
Zero => (0.0, 0.0),
One => (1.0, 0.0));
procedure Solve procedure Transpose is new Ops.Transpose
(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
(Scalar => Complex, (Scalar => Complex,
Matrix => Complex_Matrix); Matrix => Complex_Matrix);
...@@ -98,6 +67,12 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -98,6 +67,12 @@ package body Ada.Numerics.Generic_Complex_Arrays is
function Length is new Square_Matrix_Length (Complex, Complex_Matrix); 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 -- Instantiating the following subprograms directly would lead to
-- name clashes, so use a local package. -- name clashes, so use a local package.
...@@ -155,6 +130,14 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -155,6 +130,14 @@ package body Ada.Numerics.Generic_Complex_Arrays is
Right_Vector => Complex_Vector, Right_Vector => Complex_Vector,
Zero => (0.0, 0.0)); 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 function "*" is new Outer_Product
(Left_Scalar => Complex, (Left_Scalar => Complex,
Right_Scalar => Complex, Right_Scalar => Complex,
...@@ -229,6 +212,15 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -229,6 +212,15 @@ package body Ada.Numerics.Generic_Complex_Arrays is
Result_Vector => Complex_Vector, Result_Vector => Complex_Vector,
Zero => (0.0, 0.0)); 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 function "*" is new Vector_Matrix_Product
(Left_Scalar => Real'Base, (Left_Scalar => Real'Base,
Right_Scalar => Complex, Right_Scalar => Complex,
...@@ -247,6 +239,24 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -247,6 +239,24 @@ package body Ada.Numerics.Generic_Complex_Arrays is
Result_Vector => Complex_Vector, Result_Vector => Complex_Vector,
Zero => (0.0, 0.0)); 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 function "*" is new Matrix_Matrix_Product
(Left_Scalar => Real'Base, (Left_Scalar => Real'Base,
Right_Scalar => Complex, Right_Scalar => Complex,
...@@ -445,6 +455,15 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -445,6 +455,15 @@ package body Ada.Numerics.Generic_Complex_Arrays is
Result_Matrix => Complex_Matrix, Result_Matrix => Complex_Matrix,
Operation => "/"); Operation => "/");
-----------
-- "abs" --
-----------
function "abs" is new L2_Norm
(X_Scalar => Complex,
Result_Real => Real'Base,
X_Vector => Complex_Vector);
-------------- --------------
-- Argument -- -- Argument --
-------------- --------------
...@@ -671,6 +690,16 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -671,6 +690,16 @@ package body Ada.Numerics.Generic_Complex_Arrays is
Y_Matrix => Real_Matrix, Y_Matrix => Real_Matrix,
Update => Set_Re); 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 -- -- Unit_Matrix --
----------------- -----------------
...@@ -686,7 +715,6 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -686,7 +715,6 @@ package body Ada.Numerics.Generic_Complex_Arrays is
Vector => Complex_Vector, Vector => Complex_Vector,
Zero => (0.0, 0.0), Zero => (0.0, 0.0),
One => (1.0, 0.0)); One => (1.0, 0.0));
end Instantiations; end Instantiations;
--------- ---------
...@@ -696,15 +724,7 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -696,15 +724,7 @@ package body Ada.Numerics.Generic_Complex_Arrays is
function "*" function "*"
(Left : Complex_Vector; (Left : Complex_Vector;
Right : Complex_Vector) return Complex Right : Complex_Vector) return Complex
is renames Instantiations."*";
begin
if Left'Length /= Right'Length then
raise Constraint_Error with
"vectors are of different length in inner product";
end if;
return dot (Left'Length, X => Left, Y => Right);
end "*";
function "*" function "*"
(Left : Real_Vector; (Left : Real_Vector;
...@@ -738,31 +758,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -738,31 +758,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is
function "*" function "*"
(Left : Complex_Matrix; (Left : Complex_Matrix;
Right : Complex_Matrix) Right : Complex_Matrix) return Complex_Matrix
return Complex_Matrix renames Instantiations."*";
is
R : Complex_Matrix (Left'Range (1), Right'Range (2));
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 "*" function "*"
(Left : Complex_Vector; (Left : Complex_Vector;
...@@ -772,48 +769,12 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -772,48 +769,12 @@ package body Ada.Numerics.Generic_Complex_Arrays is
function "*" function "*"
(Left : Complex_Vector; (Left : Complex_Vector;
Right : Complex_Matrix) return Complex_Vector Right : Complex_Matrix) return Complex_Vector
is renames Instantiations."*";
R : Complex_Vector (Right'Range (2));
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 "*" function "*"
(Left : Complex_Matrix; (Left : Complex_Matrix;
Right : Complex_Vector) return Complex_Vector Right : Complex_Vector) return Complex_Vector
is renames Instantiations."*";
R : Complex_Vector (Left'Range (1));
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 "*" function "*"
(Left : Real_Matrix; (Left : Real_Matrix;
...@@ -984,10 +945,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -984,10 +945,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is
-- "abs" -- -- "abs" --
----------- -----------
function "abs" (Right : Complex_Vector) return Complex is function "abs" (Right : Complex_Vector) return Real'Base
begin renames Instantiations."abs";
return (nrm2 (Right'Length, Right), 0.0);
end "abs";
-------------- --------------
-- Argument -- -- Argument --
...@@ -1070,38 +1029,12 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -1070,38 +1029,12 @@ package body Ada.Numerics.Generic_Complex_Arrays is
----------------- -----------------
function Determinant (A : Complex_Matrix) return Complex is function Determinant (A : Complex_Matrix) return Complex is
N : constant Integer := Length (A); M : Complex_Matrix := A;
LU : Complex_Matrix (1 .. N, 1 .. N) := A; B : Complex_Matrix (A'Range (1), 1 .. 0);
Piv : Integer_Vector (1 .. N); R : Complex;
Info : aliased Integer := -1;
Neg : Boolean;
Det : Complex;
begin begin
if N = 0 then Forward_Eliminate (M, B, R);
return (0.0, 0.0); return R;
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;
end Determinant; end Determinant;
----------------- -----------------
...@@ -1113,174 +1046,96 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -1113,174 +1046,96 @@ package body Ada.Numerics.Generic_Complex_Arrays is
Values : out Real_Vector; Values : out Real_Vector;
Vectors : out Complex_Matrix) Vectors : out Complex_Matrix)
is is
Job_Z : aliased Character := 'V'; N : constant Natural := Length (A);
Rng : aliased Character := 'A';
Uplo : aliased Character := 'U'; -- 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
N : constant Natural := Length (A); -- eigenvalue problem:
W : BLAS_Real_Vector (Values'Range); -- (A + i * B) * (u + i * v) = Lambda * (u + i * v)
M : Integer; --
B : Complex_Matrix (1 .. N, 1 .. N); -- is equivalent to the (2 * N, 2 * N) real eigenvalue problem:
L_Work : Complex_Vector (1 .. 1); -- [ A, B ] [ u ] = Lambda * [ u ]
LR_Work : BLAS_Real_Vector (1 .. 1); -- [ -B, A ] [ v ] [ v ]
LI_Work : Integer_Vector (1 .. 1); --
I_Supp_Z : Integer_Vector (1 .. 2 * N); -- Note that the (2 * N, 2 * N) matrix above is symmetric, as
Info : aliased Integer; -- Transpose (A) = A and Transpose (B) = -B if C is Hermitian.
-- 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.
M : Real_Matrix (1 .. 2 * N, 1 .. 2 * N);
Vals : Real_Vector (1 .. 2 * N);
Vecs : Real_Matrix (1 .. 2 * N, 1 .. 2 * N);
begin begin
if Values'Length /= N then for J in 1 .. N loop
raise Constraint_Error with "wrong length for output vector"; for K in 1 .. N loop
end if; declare
C : constant Complex :=
if Vectors'First (1) /= A'First (1) (A (A'First (1) + (J - 1), A'First (2) + (K - 1)));
or else Vectors'Last (1) /= A'Last (1) begin
or else Vectors'First (2) /= A'First (2) M (J, K) := Re (C);
or else Vectors'Last (2) /= A'Last (2) M (J + N, K + N) := Re (C);
then M (J + N, K) := Im (C);
raise Constraint_Error with "wrong dimensions for output matrix"; M (J, K + N) := -Im (C);
end if; end;
if N = 0 then
return;
end if;
-- Check for hermitian matrix ???
-- Copy only required triangle ???
B := A;
-- 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 loop;
end; end loop;
Eigensystem (M, Vals, Vecs);
for J in 1 .. N loop
declare
Col : constant Integer := Values'First + (J - 1);
begin
Values (Col) := Vals (2 * J);
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; end Eigensystem;
----------------- -----------------
-- Eigenvalues -- -- Eigenvalues --
----------------- -----------------
procedure Eigenvalues function Eigenvalues (A : Complex_Matrix) return Real_Vector is
(A : Complex_Matrix; -- See Eigensystem for a description of the algorithm
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;
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 begin
if Values'Length /= N then for J in 1 .. N loop
raise Constraint_Error with "wrong length for output vector"; for K in 1 .. N loop
end if; declare
C : constant Complex :=
if N = 0 then (A (A'First (1) + (J - 1), A'First (2) + (K - 1)));
return; begin
end if; M (J, K) := Re (C);
M (J + N, K + N) := Re (C);
-- Check for hermitian matrix ??? M (J + N, K) := Im (C);
M (J, K + N) := -Im (C);
-- Find size of work area end;
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 loop;
end; end loop;
end Eigenvalues;
Vals := Eigenvalues (M);
for J in 1 .. N loop
R (A'First (1) + (J - 1)) := Vals (2 * J);
end loop;
function Eigenvalues (A : Complex_Matrix) return Real_Vector is
R : Real_Vector (A'Range (1));
begin
Eigenvalues (A, R);
return R; return R;
end Eigenvalues; end Eigenvalues;
...@@ -1298,73 +1153,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -1298,73 +1153,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is
-- Inverse -- -- 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 function Inverse (A : Complex_Matrix) return Complex_Matrix is
R : Complex_Matrix (A'Range (2), A'Range (1)); (Solve (A, Unit_Matrix (Length (A))));
begin
Inverse (A, R);
return R;
end Inverse;
------------- -------------
-- Modulus -- -- Modulus --
...@@ -1418,53 +1208,15 @@ package body Ada.Numerics.Generic_Complex_Arrays is ...@@ -1418,53 +1208,15 @@ package body Ada.Numerics.Generic_Complex_Arrays is
-- Solve -- -- Solve --
----------- -----------
procedure 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;
-- ??? 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 function Solve
(A : Complex_Matrix; (A : Complex_Matrix;
X : Complex_Vector) return Complex_Vector X : Complex_Vector) return Complex_Vector
is renames Instantiations.Solve;
B : Complex_Vector (A'Range (2));
begin
Solve (A, X, B);
return B;
end Solve;
function Solve (A, X : Complex_Matrix) return Complex_Matrix is function Solve
B : Complex_Matrix (A'Range (2), X'Range (2)); (A : Complex_Matrix;
begin X : Complex_Matrix) return Complex_Matrix
Solve (A, X, B); renames Instantiations.Solve;
return B;
end Solve;
--------------- ---------------
-- Transpose -- -- Transpose --
......
gcc/ada/a-ngcoar.ads
...@@ -66,7 +66,7 @@ package Ada.Numerics.Generic_Complex_Arrays is ...@@ -66,7 +66,7 @@ package Ada.Numerics.Generic_Complex_Arrays is
function "+" (Left, Right : Complex_Vector) return Complex_Vector; function "+" (Left, Right : Complex_Vector) return Complex_Vector;
function "-" (Left, Right : Complex_Vector) return Complex_Vector; function "-" (Left, Right : Complex_Vector) return Complex_Vector;
function "*" (Left, Right : Complex_Vector) return Complex; 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 -- Mixed Real_Vector and Complex_Vector arithmetic operations
......
gcc/ada/a-ngrear.adb
...@@ -33,7 +33,7 @@ ...@@ -33,7 +33,7 @@
-- reason for this is new Ada 2012 requirements that prohibit algorithms such -- 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 -- as Strassen's algorithm, which may be used by some BLAS implementations. In
-- addition, some platforms lacked suitable compilers to compile the reference -- addition, some platforms lacked suitable compilers to compile the reference
-- BLAS/LAPACK implementation. Finally, on some platforms there are be more -- BLAS/LAPACK implementation. Finally, on some platforms there are more
-- floating point types than supported by BLAS/LAPACK. -- floating point types than supported by BLAS/LAPACK.
with Ada.Containers.Generic_Anonymous_Array_Sort; use Ada.Containers; with Ada.Containers.Generic_Anonymous_Array_Sort; use Ada.Containers;
...@@ -337,6 +337,11 @@ package body Ada.Numerics.Generic_Real_Arrays is ...@@ -337,6 +337,11 @@ package body Ada.Numerics.Generic_Real_Arrays is
Result_Matrix => Real_Matrix, Result_Matrix => Real_Matrix,
Operation => "abs"); 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 function Unit_Matrix is new
Generic_Array_Operations.Unit_Matrix Generic_Array_Operations.Unit_Matrix
(Scalar => Real'Base, (Scalar => Real'Base,
...@@ -696,58 +701,11 @@ package body Ada.Numerics.Generic_Real_Arrays is ...@@ -696,58 +701,11 @@ package body Ada.Numerics.Generic_Real_Arrays is
-- Solve -- -- Solve --
----------- -----------
function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector is function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector
N : constant Natural := Length (A); renames Instantiations.Solve;
MA : Real_Matrix := A;
MX : Real_Matrix (A'Range (1), 1 .. 1);
R : Real_Vector (A'Range (2));
Det : Real'Base;
begin
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 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; function Solve (A, X : Real_Matrix) return Real_Matrix
end Solve; renames Instantiations.Solve;
---------------------- ----------------------
-- Sort_Eigensystem -- -- Sort_Eigensystem --
......
gcc/ada/s-gearop.adb
...@@ -651,6 +651,75 @@ package body System.Generic_Array_Operations is ...@@ -651,6 +651,75 @@ package body System.Generic_Array_Operations is
return R; return R;
end Matrix_Matrix_Product; 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 -- -- Matrix_Vector_Product --
--------------------------- ---------------------------
......
gcc/ada/s-gearop.ads
...@@ -390,6 +390,35 @@ pragma Pure (Generic_Array_Operations); ...@@ -390,6 +390,35 @@ pragma Pure (Generic_Array_Operations);
(Left : Left_Matrix; (Left : Left_Matrix;
Right : Right_Matrix) return Result_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 -- -- Sqrt --
---------- ----------
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment