| ------------------------------------------------------------------------------ |
| -- -- |
| -- 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; |