785 lines
23 KiB
Ada
785 lines
23 KiB
Ada
------------------------------------------------------------------------------
|
|
-- --
|
|
-- GNAT RUN-TIME COMPONENTS --
|
|
-- --
|
|
-- ADA.NUMERICS.GENERIC_REAL_ARRAYS --
|
|
-- --
|
|
-- B o d y --
|
|
-- --
|
|
-- Copyright (C) 2006-2009, 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- --
|
|
-- ware Foundation; either version 3, or (at your option) any later ver- --
|
|
-- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
|
|
-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
|
|
-- or FITNESS FOR A PARTICULAR PURPOSE. --
|
|
-- --
|
|
-- As a special exception under Section 7 of GPL version 3, you are granted --
|
|
-- additional permissions described in the GCC Runtime Library Exception, --
|
|
-- version 3.1, as published by the Free Software Foundation. --
|
|
-- --
|
|
-- You should have received a copy of the GNU General Public License and --
|
|
-- a copy of the GCC Runtime Library Exception along with this program; --
|
|
-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
|
|
-- <http://www.gnu.org/licenses/>. --
|
|
-- --
|
|
-- GNAT was originally developed by the GNAT team at New York University. --
|
|
-- Extensive contributions were provided by Ada Core Technologies Inc. --
|
|
-- --
|
|
------------------------------------------------------------------------------
|
|
|
|
with System; use System;
|
|
with System.Generic_Real_BLAS;
|
|
with System.Generic_Real_LAPACK;
|
|
with System.Generic_Array_Operations; use System.Generic_Array_Operations;
|
|
|
|
package body Ada.Numerics.Generic_Real_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 BLAS is
|
|
new Generic_Real_BLAS (Real'Base, Real_Vector, Real_Matrix);
|
|
|
|
package LAPACK is
|
|
new Generic_Real_LAPACK (Real'Base, Real_Vector, Real_Matrix);
|
|
|
|
use BLAS, LAPACK;
|
|
|
|
-- Procedure versions of functions returning unconstrained values.
|
|
-- This allows for inlining the function wrapper.
|
|
|
|
procedure Eigenvalues (A : Real_Matrix; Values : out Real_Vector);
|
|
procedure Inverse (A : Real_Matrix; R : out Real_Matrix);
|
|
procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector);
|
|
procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix);
|
|
|
|
procedure Transpose is new
|
|
Generic_Array_Operations.Transpose
|
|
(Scalar => Real'Base,
|
|
Matrix => Real_Matrix);
|
|
|
|
-- Helper function that raises a Constraint_Error is the argument is
|
|
-- not a square matrix, and otherwise returns its length.
|
|
|
|
function Length is new Square_Matrix_Length (Real'Base, Real_Matrix);
|
|
|
|
-- Instantiating the following subprograms directly would lead to
|
|
-- name clashes, so use a local package.
|
|
|
|
package Instantiations is
|
|
|
|
function "+" is new
|
|
Vector_Elementwise_Operation
|
|
(X_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
X_Vector => Real_Vector,
|
|
Result_Vector => Real_Vector,
|
|
Operation => "+");
|
|
|
|
function "+" is new
|
|
Matrix_Elementwise_Operation
|
|
(X_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
X_Matrix => Real_Matrix,
|
|
Result_Matrix => Real_Matrix,
|
|
Operation => "+");
|
|
|
|
function "+" is new
|
|
Vector_Vector_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Left_Vector => Real_Vector,
|
|
Right_Vector => Real_Vector,
|
|
Result_Vector => Real_Vector,
|
|
Operation => "+");
|
|
|
|
function "+" is new
|
|
Matrix_Matrix_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Left_Matrix => Real_Matrix,
|
|
Right_Matrix => Real_Matrix,
|
|
Result_Matrix => Real_Matrix,
|
|
Operation => "+");
|
|
|
|
function "-" is new
|
|
Vector_Elementwise_Operation
|
|
(X_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
X_Vector => Real_Vector,
|
|
Result_Vector => Real_Vector,
|
|
Operation => "-");
|
|
|
|
function "-" is new
|
|
Matrix_Elementwise_Operation
|
|
(X_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
X_Matrix => Real_Matrix,
|
|
Result_Matrix => Real_Matrix,
|
|
Operation => "-");
|
|
|
|
function "-" is new
|
|
Vector_Vector_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Left_Vector => Real_Vector,
|
|
Right_Vector => Real_Vector,
|
|
Result_Vector => Real_Vector,
|
|
Operation => "-");
|
|
|
|
function "-" is new
|
|
Matrix_Matrix_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Left_Matrix => Real_Matrix,
|
|
Right_Matrix => Real_Matrix,
|
|
Result_Matrix => Real_Matrix,
|
|
Operation => "-");
|
|
|
|
function "*" is new
|
|
Scalar_Vector_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Right_Vector => Real_Vector,
|
|
Result_Vector => Real_Vector,
|
|
Operation => "*");
|
|
|
|
function "*" is new
|
|
Scalar_Matrix_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Right_Matrix => Real_Matrix,
|
|
Result_Matrix => Real_Matrix,
|
|
Operation => "*");
|
|
|
|
function "*" is new
|
|
Vector_Scalar_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Left_Vector => Real_Vector,
|
|
Result_Vector => Real_Vector,
|
|
Operation => "*");
|
|
|
|
function "*" is new
|
|
Matrix_Scalar_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Left_Matrix => Real_Matrix,
|
|
Result_Matrix => Real_Matrix,
|
|
Operation => "*");
|
|
|
|
function "*" is new
|
|
Outer_Product
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Left_Vector => Real_Vector,
|
|
Right_Vector => Real_Vector,
|
|
Matrix => Real_Matrix);
|
|
|
|
function "/" is new
|
|
Vector_Scalar_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Left_Vector => Real_Vector,
|
|
Result_Vector => Real_Vector,
|
|
Operation => "/");
|
|
|
|
function "/" is new
|
|
Matrix_Scalar_Elementwise_Operation
|
|
(Left_Scalar => Real'Base,
|
|
Right_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
Left_Matrix => Real_Matrix,
|
|
Result_Matrix => Real_Matrix,
|
|
Operation => "/");
|
|
|
|
function "abs" is new
|
|
Vector_Elementwise_Operation
|
|
(X_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
X_Vector => Real_Vector,
|
|
Result_Vector => Real_Vector,
|
|
Operation => "abs");
|
|
|
|
function "abs" is new
|
|
Matrix_Elementwise_Operation
|
|
(X_Scalar => Real'Base,
|
|
Result_Scalar => Real'Base,
|
|
X_Matrix => Real_Matrix,
|
|
Result_Matrix => Real_Matrix,
|
|
Operation => "abs");
|
|
|
|
function Unit_Matrix is new
|
|
Generic_Array_Operations.Unit_Matrix
|
|
(Scalar => Real'Base,
|
|
Matrix => Real_Matrix,
|
|
Zero => 0.0,
|
|
One => 1.0);
|
|
|
|
function Unit_Vector is new
|
|
Generic_Array_Operations.Unit_Vector
|
|
(Scalar => Real'Base,
|
|
Vector => Real_Vector,
|
|
Zero => 0.0,
|
|
One => 1.0);
|
|
|
|
end Instantiations;
|
|
|
|
---------
|
|
-- "+" --
|
|
---------
|
|
|
|
function "+" (Right : Real_Vector) return Real_Vector
|
|
renames Instantiations."+";
|
|
|
|
function "+" (Right : Real_Matrix) return Real_Matrix
|
|
renames Instantiations."+";
|
|
|
|
function "+" (Left, Right : Real_Vector) return Real_Vector
|
|
renames Instantiations."+";
|
|
|
|
function "+" (Left, Right : Real_Matrix) return Real_Matrix
|
|
renames Instantiations."+";
|
|
|
|
---------
|
|
-- "-" --
|
|
---------
|
|
|
|
function "-" (Right : Real_Vector) return Real_Vector
|
|
renames Instantiations."-";
|
|
|
|
function "-" (Right : Real_Matrix) return Real_Matrix
|
|
renames Instantiations."-";
|
|
|
|
function "-" (Left, Right : Real_Vector) return Real_Vector
|
|
renames Instantiations."-";
|
|
|
|
function "-" (Left, Right : Real_Matrix) return Real_Matrix
|
|
renames Instantiations."-";
|
|
|
|
---------
|
|
-- "*" --
|
|
---------
|
|
|
|
-- Scalar multiplication
|
|
|
|
function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector
|
|
renames Instantiations."*";
|
|
|
|
function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector
|
|
renames Instantiations."*";
|
|
|
|
function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix
|
|
renames Instantiations."*";
|
|
|
|
function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
|
|
renames Instantiations."*";
|
|
|
|
-- Vector multiplication
|
|
|
|
function "*" (Left, Right : Real_Vector) return Real'Base is
|
|
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 "*" (Left, Right : Real_Vector) return Real_Matrix
|
|
renames Instantiations."*";
|
|
|
|
function "*"
|
|
(Left : Real_Vector;
|
|
Right : Real_Matrix) return Real_Vector
|
|
is
|
|
R : Real_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 "*"
|
|
(Left : Real_Matrix;
|
|
Right : Real_Vector) return Real_Vector
|
|
is
|
|
R : Real_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 "*";
|
|
|
|
-- Matrix Multiplication
|
|
|
|
function "*" (Left, Right : Real_Matrix) return Real_Matrix is
|
|
R : Real_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 "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector
|
|
renames Instantiations."/";
|
|
|
|
function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
|
|
renames Instantiations."/";
|
|
|
|
-----------
|
|
-- "abs" --
|
|
-----------
|
|
|
|
function "abs" (Right : Real_Vector) return Real'Base is
|
|
begin
|
|
return nrm2 (Right'Length, Right);
|
|
end "abs";
|
|
|
|
function "abs" (Right : Real_Vector) return Real_Vector
|
|
renames Instantiations."abs";
|
|
|
|
function "abs" (Right : Real_Matrix) return Real_Matrix
|
|
renames Instantiations."abs";
|
|
|
|
-----------------
|
|
-- Determinant --
|
|
-----------------
|
|
|
|
function Determinant (A : Real_Matrix) return Real'Base is
|
|
N : constant Integer := Length (A);
|
|
LU : Real_Matrix (1 .. N, 1 .. N) := A;
|
|
Piv : Integer_Vector (1 .. N);
|
|
Info : aliased Integer := -1;
|
|
Det : Real := 1.0;
|
|
|
|
begin
|
|
getrf (M => N,
|
|
N => N,
|
|
A => LU,
|
|
Ld_A => N,
|
|
I_Piv => Piv,
|
|
Info => Info'Access);
|
|
|
|
if Info /= 0 then
|
|
raise Constraint_Error with "ill-conditioned matrix";
|
|
end if;
|
|
|
|
for J in 1 .. N loop
|
|
Det := (if Piv (J) /= J then -Det * LU (J, J) else Det * LU (J, J));
|
|
end loop;
|
|
|
|
return Det;
|
|
end Determinant;
|
|
|
|
-----------------
|
|
-- Eigensystem --
|
|
-----------------
|
|
|
|
procedure Eigensystem
|
|
(A : Real_Matrix;
|
|
Values : out Real_Vector;
|
|
Vectors : out Real_Matrix)
|
|
is
|
|
N : constant Natural := Length (A);
|
|
Tau : Real_Vector (1 .. N);
|
|
L_Work : Real_Vector (1 .. 1);
|
|
Info : aliased Integer;
|
|
|
|
E : Real_Vector (1 .. N);
|
|
pragma Warnings (Off, E);
|
|
|
|
begin
|
|
if Values'Length /= N then
|
|
raise Constraint_Error with "wrong length for output vector";
|
|
end if;
|
|
|
|
if N = 0 then
|
|
return;
|
|
end if;
|
|
|
|
-- Initialize working matrix and check for symmetric input matrix
|
|
|
|
Transpose (A, Vectors);
|
|
|
|
if A /= Vectors then
|
|
raise Argument_Error with "matrix not symmetric";
|
|
end if;
|
|
|
|
-- Compute size of additional working space
|
|
|
|
sytrd (Uplo => Lower'Access,
|
|
N => N,
|
|
A => Vectors,
|
|
Ld_A => N,
|
|
D => Values,
|
|
E => E,
|
|
Tau => Tau,
|
|
Work => L_Work,
|
|
L_Work => -1,
|
|
Info => Info'Access);
|
|
|
|
declare
|
|
Work : Real_Vector (1 .. Integer'Max (Integer (L_Work (1)), 2 * N));
|
|
pragma Warnings (Off, Work);
|
|
|
|
Comp_Z : aliased constant Character := 'V';
|
|
|
|
begin
|
|
-- Reduce matrix to tridiagonal form
|
|
|
|
sytrd (Uplo => Lower'Access,
|
|
N => N,
|
|
A => Vectors,
|
|
Ld_A => A'Length (1),
|
|
D => Values,
|
|
E => E,
|
|
Tau => Tau,
|
|
Work => Work,
|
|
L_Work => Work'Length,
|
|
Info => Info'Access);
|
|
|
|
if Info /= 0 then
|
|
raise Program_Error;
|
|
end if;
|
|
|
|
-- Generate the real orthogonal matrix determined by sytrd
|
|
|
|
orgtr (Uplo => Lower'Access,
|
|
N => N,
|
|
A => Vectors,
|
|
Ld_A => N,
|
|
Tau => Tau,
|
|
Work => Work,
|
|
L_Work => Work'Length,
|
|
Info => Info'Access);
|
|
|
|
if Info /= 0 then
|
|
raise Program_Error;
|
|
end if;
|
|
|
|
-- Compute all eigenvalues and eigenvectors using QR algorithm
|
|
|
|
steqr (Comp_Z => Comp_Z'Access,
|
|
N => N,
|
|
D => Values,
|
|
E => E,
|
|
Z => Vectors,
|
|
Ld_Z => N,
|
|
Work => Work,
|
|
Info => Info'Access);
|
|
|
|
if Info /= 0 then
|
|
raise Constraint_Error with
|
|
"eigensystem computation failed to converge";
|
|
end if;
|
|
end;
|
|
end Eigensystem;
|
|
|
|
-----------------
|
|
-- Eigenvalues --
|
|
-----------------
|
|
|
|
procedure Eigenvalues
|
|
(A : Real_Matrix;
|
|
Values : out Real_Vector)
|
|
is
|
|
N : constant Natural := Length (A);
|
|
L_Work : Real_Vector (1 .. 1);
|
|
Info : aliased Integer;
|
|
|
|
B : Real_Matrix (1 .. N, 1 .. N);
|
|
Tau : Real_Vector (1 .. N);
|
|
E : Real_Vector (1 .. N);
|
|
pragma Warnings (Off, B);
|
|
pragma Warnings (Off, Tau);
|
|
pragma Warnings (Off, E);
|
|
|
|
begin
|
|
if Values'Length /= N then
|
|
raise Constraint_Error with "wrong length for output vector";
|
|
end if;
|
|
|
|
if N = 0 then
|
|
return;
|
|
end if;
|
|
|
|
-- Initialize working matrix and check for symmetric input matrix
|
|
|
|
Transpose (A, B);
|
|
|
|
if A /= B then
|
|
raise Argument_Error with "matrix not symmetric";
|
|
end if;
|
|
|
|
-- Find size of work area
|
|
|
|
sytrd (Uplo => Lower'Access,
|
|
N => N,
|
|
A => B,
|
|
Ld_A => N,
|
|
D => Values,
|
|
E => E,
|
|
Tau => Tau,
|
|
Work => L_Work,
|
|
L_Work => -1,
|
|
Info => Info'Access);
|
|
|
|
declare
|
|
Work : Real_Vector (1 .. Integer'Min (Integer (L_Work (1)), 4 * N));
|
|
pragma Warnings (Off, Work);
|
|
|
|
begin
|
|
-- Reduce matrix to tridiagonal form
|
|
|
|
sytrd (Uplo => Lower'Access,
|
|
N => N,
|
|
A => B,
|
|
Ld_A => A'Length (1),
|
|
D => Values,
|
|
E => E,
|
|
Tau => Tau,
|
|
Work => Work,
|
|
L_Work => Work'Length,
|
|
Info => Info'Access);
|
|
|
|
if Info /= 0 then
|
|
raise Constraint_Error;
|
|
end if;
|
|
|
|
-- Compute all eigenvalues using QR algorithm
|
|
|
|
sterf (N => N,
|
|
D => Values,
|
|
E => E,
|
|
Info => Info'Access);
|
|
|
|
if Info /= 0 then
|
|
raise Constraint_Error with
|
|
"eigenvalues computation failed to converge";
|
|
end if;
|
|
end;
|
|
end Eigenvalues;
|
|
|
|
function Eigenvalues (A : Real_Matrix) return Real_Vector is
|
|
R : Real_Vector (A'Range (1));
|
|
begin
|
|
Eigenvalues (A, R);
|
|
return R;
|
|
end Eigenvalues;
|
|
|
|
-------------
|
|
-- Inverse --
|
|
-------------
|
|
|
|
procedure Inverse (A : Real_Matrix; R : out Real_Matrix) is
|
|
N : constant Integer := Length (A);
|
|
Piv : Integer_Vector (1 .. N);
|
|
L_Work : Real_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 : Real_Vector (1 .. Integer (L_Work (1)));
|
|
pragma Warnings (Off, Work);
|
|
|
|
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 : Real_Matrix) return Real_Matrix is
|
|
R : Real_Matrix (A'Range (2), A'Range (1));
|
|
begin
|
|
Inverse (A, R);
|
|
return R;
|
|
end Inverse;
|
|
|
|
-----------
|
|
-- Solve --
|
|
-----------
|
|
|
|
procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_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 : Real_Matrix; X : Real_Matrix; B : out Real_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 : Real_Matrix; X : Real_Vector) return Real_Vector is
|
|
B : Real_Vector (A'Range (2));
|
|
begin
|
|
Solve (A, X, B);
|
|
return B;
|
|
end Solve;
|
|
|
|
function Solve (A, X : Real_Matrix) return Real_Matrix is
|
|
B : Real_Matrix (A'Range (2), X'Range (2));
|
|
begin
|
|
Solve (A, X, B);
|
|
return B;
|
|
end Solve;
|
|
|
|
---------------
|
|
-- Transpose --
|
|
---------------
|
|
|
|
function Transpose (X : Real_Matrix) return Real_Matrix is
|
|
R : Real_Matrix (X'Range (2), X'Range (1));
|
|
begin
|
|
Transpose (X, R);
|
|
|
|
return R;
|
|
end Transpose;
|
|
|
|
-----------------
|
|
-- Unit_Matrix --
|
|
-----------------
|
|
|
|
function Unit_Matrix
|
|
(Order : Positive;
|
|
First_1 : Integer := 1;
|
|
First_2 : Integer := 1) return Real_Matrix
|
|
renames Instantiations.Unit_Matrix;
|
|
|
|
-----------------
|
|
-- Unit_Vector --
|
|
-----------------
|
|
|
|
function Unit_Vector
|
|
(Index : Integer;
|
|
Order : Positive;
|
|
First : Integer := 1) return Real_Vector
|
|
renames Instantiations.Unit_Vector;
|
|
|
|
end Ada.Numerics.Generic_Real_Arrays;
|