+2011-08-29 Hristian Kirtchev <kirtchev@adacore.com>
+
+ * exp_ch4.adb (Expand_Allocator_Expression): Add code to set attribute
+ Finalize_Address of the access type's finalization master.
+ (Expand_N_Allocator): Add code to set attribute Finalize_Address of the
+ access type's finalization master. Add a guard to prevent
+ Associated_Storage_Pool from being set on .NET/JVM.
+ * exp_ch6.adb (Make_Build_In_Place_Call_In_Allocator): Add code to set
+ attribute Finalize_Address of the access type's finalization master.
+ * exp_ch7.adb (Make_Finalize_Address_Call): New routine.
+ * exp_ch7.ads (Make_Finalize_Address_Call): New routine.
+ * rtsfind.ads: Add RE_Set_Finalize_Address to tables RE_Id and
+ RE_Unit_Table.
+ * s-finmas.adb: Add with clause for System.Address_Image. Add with and
+ use clause for System.IO
+ (Detach): Relax the assertion, to be reinstated later.
+ (Finalize): Rewrite the iteration loop to avoid pointer comparison.
+ Relax the assertion on Finalize_Address, to be reinstated later.
+ (Is_Empty_List): New routine.
+ (pm): New debug routine.
+ (Set_Finalize_Address): New routine.
+ * s-finmas.ads (pm): New debug routine.
+ (Set_Finalize_Address): New routine.
+ * s-stposu.adb (Allocate_Any_Controlled): Code reformatting.
+
+2011-08-29 Tristan Gingold <gingold@adacore.com>
+
+ * a-exexpr-gcc.adb (GCC_Exception_Access, GNAT_GCC_Exception_Access):
+ Remove convention C.
+
+2011-08-29 Tristan Gingold <gingold@adacore.com>
+
+ * s-taprop-vms.adb (Get_Exc_Stack_Addr): Remove.
+ (Initialize_TCB): Remove Exc_Stack_Ptr initialization.
+ (Finalize_TCB): Remove its finalization.
+ (Initialize): Remove assignment of GET_Exc_Stack_Addr
+ * s-soflin.adb (NT_Exc_Stack): Remove
+ (Get_Exc_Stack_Addr_NT): Likewise.
+ (Get_Exc_Stack_Addr_Soft): Likewise.
+ * s-soflin.ads (Get_Exc_Stack_Addr_NT): Remove.
+ (Get_Exc_Stack_Addr): Likewise.
+ (Get_Exc_Stack_Addr_Soft): Likewise
+ * s-taspri-vms.ads (Exc_Stack_T): Remove.
+ (Exc_Stack_Ptr_T): Likewise.
+ (Private_Data): Remove Exc_Stack_Ptr component.
+
+2011-08-29 Tristan Gingold <gingold@adacore.com>
+
+ * raise-gcc.c (get_ip_from_context): New function. Factorize code.
+
+2011-08-29 Tristan Gingold <gingold@adacore.com>
+
+ * gnat_ugn.texi: Fix aix and x86-solaris info for run-time.
+
+2011-08-29 Geert Bosch <bosch@adacore.com>
+
+ * s-gearop.ads (Back_Substitute, Diagonal, Forward_Eliminate,
+ L2_Norm, Swap_Column): New generic subprograms
+ * s-gearop.adb (Back_Substitute, Diagonal, Forward_Eliminate,
+ L2_Norm, Swap_Column): Implement new subprograms in order to
+ eliminate dependency on BLAS and LAPACK libraries in
+ Ada.Numerics.Generic_Real_Arrays and eventually also the complex
+ version. Forward_Eliminate/Back_Substitute can be used to put a
+ matrix in row echelon or reduced row echelon form using partial
+ pivoting.
+ * a-ngrear.adb: (Back_Substitute, Diagonal, Forward_Eleminate,
+ Swap_Column): Instantiate from System.Generic_Array_Operations.
+ ("*", "abs"): Implement by instantiation from Generic_Array_Operations.
+ (Sqrt): Local function for simple computation of square root without
+ adding dependencies on Generic_Elementary_Functions.
+ (Swap): New subprogram to exchange floating point numbers.
+ (Inverse): Reimplement using Jordan-Gauss elimination.
+ (Jacobi): New procedure implementing Jacobi's method for computation
+ of eigensystems, based on Rutishauser's implementation.
+ (L2_Norm): Implement directly using the inner product.
+ (Sort_Eigensystem): Sort eigenvalue/eigenvector pairs in order of
+ decreasing eigenvalue as required by the Ada RM.
+ (Swap_Column): New helper procedure for Sort_Eigensystem.
+ Remove with of System.Generic_Real_BLAS and System.Generic_Real_LAPACK.
+ Add with of Ada.Containers.Generic_Anonymous_Array_Sort, for
+ Sort_Eigensystems.
+
+2011-08-29 Thomas Quinot <quinot@adacore.com>
+
+ * put_scos.adb (Put_SCOs): Do not emit a newline for an empty
+ statements line.
+
2011-08-29 Hristian Kirtchev <kirtchev@adacore.com>
* s-finmas.adb (Finalize): Check Finalize_Address of the master rather
-- alignment below.
type GCC_Exception_Access is access all Unwind_Exception;
- pragma Convention (C, GCC_Exception_Access);
- -- Pointer to a GCC exception
+ -- Pointer to a GCC exception. Do not use convention C as on VMS this
+ -- would imply the use of 32-bits pointers.
procedure Unwind_DeleteException (Excp : not null GCC_Exception_Access);
pragma Import (C, Unwind_DeleteException, "_Unwind_DeleteException");
-- to maintain anyway.
type GNAT_GCC_Exception_Access is access all GNAT_GCC_Exception;
- pragma Convention (C, GNAT_GCC_Exception_Access);
function To_GCC_Exception is new
Unchecked_Conversion (GNAT_GCC_Exception_Access, GCC_Exception_Access);
-- --
-- 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 --
-- terms of the GNU General Public License as published by the Free Soft- --
-- --
------------------------------------------------------------------------------
+-- This version of Generic_Real_Arrays avoids the use of BLAS and LAPACK. One
+-- 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 many platforms there may be more
+-- floating point types than supported by BLAS/LAPACK.
+
+with Ada.Containers.Generic_Anonymous_Array_Sort; use Ada.Containers;
+
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 Ops renames System.Generic_Array_Operations;
- package BLAS is
- new Generic_Real_BLAS (Real'Base, Real_Vector, Real_Matrix);
+ function Is_Non_Zero (X : Real'Base) return Boolean is (X /= 0.0);
- package LAPACK is
- new Generic_Real_LAPACK (Real'Base, Real_Vector, Real_Matrix);
+ procedure Back_Substitute is new Ops.Back_Substitute
+ (Scalar => Real'Base,
+ Matrix => Real_Matrix,
+ Is_Non_Zero => Is_Non_Zero);
- use BLAS, LAPACK;
+ function Diagonal is new Ops.Diagonal
+ (Scalar => Real'Base,
+ Vector => Real_Vector,
+ Matrix => Real_Matrix);
- -- Procedure versions of functions returning unconstrained values.
- -- This allows for inlining the function wrapper.
+ procedure Forward_Eliminate is new Ops.Forward_Eliminate
+ (Scalar => Real'Base,
+ Matrix => Real_Matrix,
+ Zero => 0.0,
+ One => 1.0);
- 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 Swap_Column is new Ops.Swap_Column
+ (Scalar => Real'Base,
+ Matrix => Real_Matrix);
- procedure Transpose is new
- Generic_Array_Operations.Transpose
+ procedure Transpose is new Ops.Transpose
(Scalar => Real'Base,
Matrix => Real_Matrix);
+ function Is_Symmetric (A : Real_Matrix) return Boolean is
+ (Transpose (A) = A);
+ -- Return True iff A is symmetric, see RM G.3.1 (90).
+
+ function Is_Tiny (Value, Compared_To : Real) return Boolean is
+ (abs Compared_To + 100.0 * abs (Value) = abs Compared_To);
+ -- Return True iff the Value is much smaller in magnitude than the least
+ -- significant digit of Compared_To.
+
+ procedure Jacobi
+ (A : Real_Matrix;
+ Values : out Real_Vector;
+ Vectors : out Real_Matrix;
+ Compute_Vectors : Boolean := True);
+ -- Perform Jacobi's eigensystem algorithm on real symmetric matrix A
+
+ function Length is new Square_Matrix_Length (Real'Base, 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);
+ procedure Rotate (X, Y : in out Real; Sin, Tau : Real);
+ -- Perform a Givens rotation
+
+ procedure Sort_Eigensystem
+ (Values : in out Real_Vector;
+ Vectors : in out Real_Matrix);
+ -- Sort Values and associated Vectors by decreasing absolute value
+
+ procedure Swap (Left, Right : in out Real);
+ -- Exchange Left and Right.
+
+ function Sqrt (X : Real) return Real;
+ -- Sqrt is implemented locally here, in order to avoid dragging in all of
+ -- the elementary functions. Speed of the square root is not a big concern
+ -- here. This also avoids depending on a specific floating point type.
+
+ ------------
+ -- Rotate --
+ ------------
+
+ procedure Rotate (X, Y : in out Real; Sin, Tau : Real) is
+ Old_X : constant Real := X;
+ Old_Y : constant Real := Y;
+ begin
+ X := Old_X - Sin * (Old_Y + Old_X * Tau);
+ Y := Old_Y + Sin * (Old_X - Old_Y * Tau);
+ end Rotate;
+
+ ----------
+ -- Sqrt --
+ ----------
+
+ function Sqrt (X : Real) return Real is
+ Root, Next : Real;
+
+ begin
+ -- Be defensive: any comparisons with NaN values will yield False.
+
+ if not (X > 0.0) then
+ if X = 0.0 then
+ return X;
+
+ else
+ raise Argument_Error;
+ end if;
+ end if;
+
+ -- Compute an initial estimate based on:
+
+ -- X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0),
+
+ -- where M is the mantissa, R is the radix and E the exponent.
+
+ -- By ignoring the mantissa and ignoring the case of an odd
+ -- exponent, we get a final error that is at most R. In other words,
+ -- the result has about a single bit precision.
+
+ Root := Real (Real'Machine_Radix) ** (Real'Exponent (X) / 2);
+
+ -- Because of the poor initial estimate, use the Babylonian method of
+ -- computing the square root, as it is stable for all inputs. Every step
+ -- will roughly double the precision of the result. Just a few steps
+ -- suffice in most cases. Eight iterations should give about 2**8 bits
+ -- of precision.
+
+ for J in 1 .. 8 loop
+ Next := (Root + X / Root) / 2.0;
+
+ exit when Root = Next;
+
+ Root := Next;
+ end loop;
+
+ return Root;
+ end Sqrt;
+
+ ----------
+ -- Swap --
+ ----------
+
+ procedure Swap (Left, Right : in out Real) is
+ Temp : constant Real := Left;
+ begin
+ Left := Right;
+ Right := Temp;
+ end Swap;
-- Instantiating the following subprograms directly would lead to
-- name clashes, so use a local package.
Right_Vector => Real_Vector,
Matrix => Real_Matrix);
+ function "*" is new
+ Inner_Product
+ (Left_Scalar => Real'Base,
+ Right_Scalar => Real'Base,
+ Result_Scalar => Real'Base,
+ Left_Vector => Real_Vector,
+ Right_Vector => Real_Vector,
+ Zero => 0.0);
+
+ function "*" is new
+ Matrix_Vector_Product
+ (Left_Scalar => Real'Base,
+ Right_Scalar => Real'Base,
+ Result_Scalar => Real'Base,
+ Matrix => Real_Matrix,
+ Right_Vector => Real_Vector,
+ Result_Vector => Real_Vector,
+ Zero => 0.0);
+
+ function "*" is new
+ Vector_Matrix_Product
+ (Left_Scalar => Real'Base,
+ Right_Scalar => Real'Base,
+ Result_Scalar => Real'Base,
+ Left_Vector => Real_Vector,
+ Matrix => Real_Matrix,
+ Result_Vector => Real_Vector,
+ Zero => 0.0);
+
+ function "*" is new
+ Matrix_Matrix_Product
+ (Left_Scalar => Real'Base,
+ Right_Scalar => Real'Base,
+ Result_Scalar => Real'Base,
+ Left_Matrix => Real_Matrix,
+ Right_Matrix => Real_Matrix,
+ Result_Matrix => Real_Matrix,
+ Zero => 0.0);
+
function "/" is new
Vector_Scalar_Elementwise_Operation
(Left_Scalar => Real'Base,
Result_Matrix => Real_Matrix,
Operation => "/");
+ function "abs" is new
+ L2_Norm
+ (Scalar => Real'Base,
+ Vector => Real_Vector,
+ Inner_Product => "*",
+ Sqrt => Sqrt);
+
function "abs" is new
Vector_Elementwise_Operation
(X_Scalar => Real'Base,
-- 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'Base
+ renames Instantiations."*";
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);
+ function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector
+ renames Instantiations."*";
- return R;
- end "*";
+ function "*" (Left : Real_Matrix; Right : Real_Vector) return Real_Vector
+ renames Instantiations."*";
-- 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, Right : Real_Matrix) 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'Base
+ renames Instantiations."abs";
function "abs" (Right : Real_Vector) return Real_Vector
renames Instantiations."abs";
-----------------
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;
+ M : Real_Matrix := A;
+ B : Real_Matrix (A'Range (1), 1 .. 0);
+ R : Real'Base;
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;
+ Forward_Eliminate (M, B, R);
- 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;
+ return R;
end Determinant;
-----------------
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);
+ Jacobi (A, Values, Vectors, Compute_Vectors => True);
+ Sort_Eigensystem (Values, Vectors);
+ end Eigensystem;
- if A /= Vectors then
- raise Argument_Error with "matrix not symmetric";
- end if;
+ -----------------
+ -- Eigenvalues --
+ -----------------
- -- Compute size of additional working space
+ function Eigenvalues (A : Real_Matrix) return Real_Vector is
+ Values : Real_Vector (A'Range (1));
+ Vectors : Real_Matrix (1 .. 0, 1 .. 0);
+ begin
+ Jacobi (A, Values, Vectors, Compute_Vectors => False);
+ Sort_Eigensystem (Values, Vectors);
- 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);
+ return Values;
+ end Eigenvalues;
- declare
- Work : Real_Vector (1 .. Integer'Max (Integer (L_Work (1)), 2 * N));
- pragma Warnings (Off, Work);
+ -------------
+ -- Inverse --
+ -------------
- Comp_Z : aliased constant Character := 'V';
+ function Inverse (A : Real_Matrix) return Real_Matrix is
+ (Solve (A, Unit_Matrix (Length (A))));
- 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;
+ ------------
+ -- Jacobi --
+ ------------
- -- Generate the real orthogonal matrix determined by sytrd
+ procedure Jacobi
+ (A : Real_Matrix;
+ Values : out Real_Vector;
+ Vectors : out Real_Matrix;
+ Compute_Vectors : Boolean := True)
+ is
+ -- This subprogram uses Carl Gustav Jacob Jacobi's iterative method
+ -- for computing eigenvalues and eigenvectors and is based on
+ -- Rutishauser's implementation.
- orgtr (Uplo => Lower'Access,
- N => N,
- A => Vectors,
- Ld_A => N,
- Tau => Tau,
- Work => Work,
- L_Work => Work'Length,
- Info => Info'Access);
+ -- The given real symmetric matrix is transformed iteratively to
+ -- diagonal form through a sequence of appropriately chosen elementary
+ -- orthogonal transformations, called Jacobi rotations here.
- if Info /= 0 then
- raise Program_Error;
- end if;
+ -- The Jacobi method produces a systematic decrease of the sum of the
+ -- squares of off-diagonal elements. Convergence to zero is quadratic,
+ -- both for this implementation, as for the classic method that doesn't
+ -- use row-wise scanning for pivot selection.
- -- Compute all eigenvalues and eigenvectors using QR algorithm
+ -- The numerical stability and accuracy of Jacobi's method make it the
+ -- best choice here, even though for large matrices other methods will
+ -- be significantly more efficient in both time and space.
- steqr (Comp_Z => Comp_Z'Access,
- N => N,
- D => Values,
- E => E,
- Z => Vectors,
- Ld_Z => N,
- Work => Work,
- Info => Info'Access);
+ -- While the eigensystem computations are absolutely foolproof for all
+ -- real symmetric matrices, in presence of invalid values, or similar
+ -- exceptional situations it might not. In such cases the results cannot
+ -- be trusted and Constraint_Error is raised.
- if Info /= 0 then
- raise Constraint_Error with
- "eigensystem computation failed to converge";
- end if;
- end;
- end Eigensystem;
+ -- Note: this implementation needs temporary storage for 2 * N + N**2
+ -- values of type Real.
- -----------------
- -- Eigenvalues --
- -----------------
+ Max_Iterations : constant := 50;
- procedure Eigenvalues
- (A : Real_Matrix;
- Values : out Real_Vector)
- is
- N : constant Natural := Length (A);
- L_Work : Real_Vector (1 .. 1);
- Info : aliased Integer;
+ N : constant Natural := Length (A);
- 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);
+ subtype Square_Matrix is Real_Matrix (1 .. N, 1 .. N);
- begin
- if Values'Length /= N then
- raise Constraint_Error with "wrong length for output vector";
- end if;
+ -- In order to annihilate the M (Row, Col) element, the
+ -- rotation parameters Cos and Sin are computed as
+ -- follows:
- if N = 0 then
- return;
- end if;
+ -- Theta = Cot (2.0 * Phi)
+ -- = (Diag (Col) - Diag (Row)) / (2.0 * M (Row, Col))
- -- Initialize working matrix and check for symmetric input matrix
+ -- Then Tan (Phi) as the smaller root (in modulus) of
- Transpose (A, B);
+ -- T**2 + 2 * T * Theta = 1 (or 0.5 / Theta, if Theta is large)
- if A /= B then
- raise Argument_Error with "matrix not symmetric";
- end if;
+ function Compute_Tan (Theta : Real) return Real is
+ (Real'Copy_Sign (1.0 / (abs Theta + Sqrt (1.0 + Theta**2)), Theta));
- -- Find size of work area
+ function Compute_Tan (P, H : Real) return Real is
+ (if Is_Tiny (P, Compared_To => H) then P / H
+ else Compute_Tan (Theta => H / (2.0 * P)));
- 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);
+ function Sum_Strict_Upper (M : Square_Matrix) return Real;
+ -- Return the sum of all elements in the strict upper triangle of M
- declare
- Work : Real_Vector (1 .. Integer'Min (Integer (L_Work (1)), 4 * N));
- pragma Warnings (Off, Work);
+ ----------------------
+ -- Sum_Strict_Upper --
+ ----------------------
+ function Sum_Strict_Upper (M : Square_Matrix) return Real is
+ Sum : Real := 0.0;
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;
+ for Row in 1 .. N - 1 loop
+ for Col in Row + 1 .. N loop
+ Sum := Sum + abs M (Row, Col);
+ end loop;
+ end loop;
+
+ return Sum;
+ end Sum_Strict_Upper;
+
+ M : Square_Matrix := A; -- Work space for solving eigensystem
+ Threshold : Real;
+ Sum : Real;
+ Diag : Real_Vector (1 .. N);
+ Diag_Adj : Real_Vector (1 .. N);
+
+ -- The vector Diag_Adj indicates the amount of change in each value,
+ -- while Diag tracks the value itself and Values holds the values as
+ -- they were at the beginning. As the changes typically will be small
+ -- compared to the absolute value of Diag, at the end of each iteration
+ -- Diag is computed as Diag + Diag_Adj thus avoiding accumulating
+ -- rounding errors. This technique is due to Rutishauser.
begin
- -- All computations are done using column-major order, but this works
- -- out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
-
- R := A;
+ if Compute_Vectors
+ and then (Vectors'Length (1) /= N or else Vectors'Length (2) /= N)
+ then
+ raise Constraint_Error with "incompatible matrix dimensions";
- -- Compute LU decomposition
+ elsif Values'Length /= N then
+ raise Constraint_Error with "incompatible vector length";
- getrf (M => N,
- N => N,
- A => R,
- Ld_A => N,
- I_Piv => Piv,
- Info => Info'Access);
+ elsif not Is_Symmetric (M) then
+ raise Constraint_Error with "matrix not symmetric";
+ end if;
- if Info /= 0 then
- raise Constraint_Error with "inverting singular matrix";
+ -- Note: Only the locally declared matrix M and vectors (Diag, Diag_Adj)
+ -- have lower bound equal to 1. The Vectors matrix may have
+ -- different bounds, so take care indexing elements. Assignment
+ -- as a whole is fine as sliding is automatic in that case.
+
+ Vectors := (if not Compute_Vectors then (1 .. 0 => (1 .. 0 => 0.0))
+ else Unit_Matrix (Vectors'Length (1), Vectors'Length (2)));
+ Values := Diagonal (M);
+
+ Sweep : for Iteration in 1 .. Max_Iterations loop
+
+ -- The first three iterations, perform rotation for any non-zero
+ -- element. After this, rotate only for those that are not much
+ -- smaller than the average off-diagnal element. After the fifth
+ -- iteration, additionally zero out off-diagonal elements that are
+ -- very small compared to elements on the diagonal with the same
+ -- column or row index.
+
+ Sum := Sum_Strict_Upper (M);
+
+ exit Sweep when Sum = 0.0;
+
+ Threshold := (if Iteration < 4 then 0.2 * Sum / Real (N**2) else 0.0);
+
+ -- Iterate over all off-diagonal elements, rotating any that have
+ -- an absolute value that exceeds the threshold.
+
+ Diag := Values;
+ Diag_Adj := (others => 0.0); -- Accumulates adjustments to Diag
+
+ for Row in 1 .. N - 1 loop
+ for Col in Row + 1 .. N loop
+
+ -- If, before the rotation M (Row, Col) is tiny compared to
+ -- Diag (Row) and Diag (Col), rotation is skipped. This is
+ -- meaningful, as it produces no larger error than would be
+ -- produced anyhow if the rotation had been performed.
+ -- Suppress this optimization in the first four sweeps, so
+ -- that this procedure can be used for computing eigenvectors
+ -- of perturbed diagonal matrices.
+
+ if Iteration > 4
+ and then Is_Tiny (M (Row, Col), Compared_To => Diag (Row))
+ and then Is_Tiny (M (Row, Col), Compared_To => Diag (Col))
+ then
+ M (Row, Col) := 0.0;
+
+ elsif abs M (Row, Col) > Threshold then
+ Perform_Rotation : declare
+ Tan : constant Real := Compute_Tan (M (Row, Col),
+ Diag (Col) - Diag (Row));
+ Cos : constant Real := 1.0 / Sqrt (1.0 + Tan**2);
+ Sin : constant Real := Tan * Cos;
+ Tau : constant Real := Sin / (1.0 + Cos);
+ Adj : constant Real := Tan * M (Row, Col);
+
+ begin
+ Diag_Adj (Row) := Diag_Adj (Row) - Adj;
+ Diag_Adj (Col) := Diag_Adj (Col) + Adj;
+ Diag (Row) := Diag (Row) - Adj;
+ Diag (Col) := Diag (Col) + Adj;
+
+ M (Row, Col) := 0.0;
+
+ for J in 1 .. Row - 1 loop -- 1 <= J < Row
+ Rotate (M (J, Row), M (J, Col), Sin, Tau);
+ end loop;
+
+ for J in Row + 1 .. Col - 1 loop -- Row < J < Col
+ Rotate (M (Row, J), M (J, Col), Sin, Tau);
+ end loop;
+
+ for J in Col + 1 .. N loop -- Col < J <= N
+ Rotate (M (Row, J), M (Col, J), Sin, Tau);
+ end loop;
+
+ for J in Vectors'Range (1) loop
+ Rotate (Vectors (J, Row - 1 + Vectors'First (2)),
+ Vectors (J, Col - 1 + Vectors'First (2)),
+ Sin, Tau);
+ end loop;
+ end Perform_Rotation;
+ end if;
+ end loop;
+ end loop;
+
+ Values := Values + Diag_Adj;
+ end loop Sweep;
+
+ -- All normal matrices with valid values should converge perfectly.
+
+ if Sum /= 0.0 then
+ raise Constraint_Error with "eigensystem solution does not converge";
end if;
+ end Jacobi;
- -- Determine size of work area
+ -----------
+ -- Solve --
+ -----------
- getri (N => N,
- A => R,
- Ld_A => N,
- I_Piv => Piv,
- Work => L_Work,
- L_Work => -1,
- Info => Info'Access);
+ 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;
- if Info /= 0 then
- raise Constraint_Error;
+ begin
+ if X'Length /= N then
+ raise Constraint_Error with "incompatible vector length";
end if;
- declare
- Work : Real_Vector (1 .. Integer (L_Work (1)));
- pragma Warnings (Off, Work);
+ for J in 0 .. MX'Length (1) - 1 loop
+ MX (MX'First (1) + J, 1) := X (X'First + J);
+ end loop;
- 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;
+ Forward_Eliminate (MA, MX, Det);
+ Back_Substitute (MA, MX);
- -- ??? Should iterate with gerfs, based on implementation advice
- end;
- end Inverse;
+ for J in 0 .. R'Length - 1 loop
+ R (R'First + J) := MX (MX'First (1) + J, 1);
+ end loop;
- 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;
+ end Solve;
- -----------
- -- 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;
- 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";
+ if X'Length (1) /= N then
+ raise Constraint_Error with "matrices have unequal number of rows";
end if;
- -- ??? Should solve directly, is faster and more accurate
+ 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);
- B := Inverse (A) * X;
+ return MB;
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;
+ ----------------------
+ -- Sort_Eigensystem --
+ ----------------------
- -- ??? Should solve directly, is faster and more accurate
+ procedure Sort_Eigensystem
+ (Values : in out Real_Vector;
+ Vectors : in out Real_Matrix)
+ is
- B := Inverse (A) * X;
- end Solve;
+ procedure Swap (Left, Right : Integer);
+ -- Swap Values (Left) with Values (Right), and also swap the
+ -- corresponding eigenvectors. Note that lowerbounds may differ.
- 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 Less (Left, Right : Integer) return Boolean is
+ (Values (Left) > Values (Right));
+ -- Sort by decreasing eigenvalue, see RM G.3.1 (76).
+
+ procedure Sort is new Generic_Anonymous_Array_Sort (Integer);
+ -- Sorts eigenvalues and eigenvectors by decreasing value
+
+ procedure Swap (Left, Right : Integer) is
+ begin
+ Swap (Values (Left), Values (Right));
+ Swap_Column (Vectors, Left - Values'First + Vectors'First (2),
+ Right - Values'First + Vectors'First (2));
+ end Swap;
- 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;
+ Sort (Values'First, Values'Last);
+ end Sort_Eigensystem;
---------------
-- Transpose --
Prefix => New_Reference_To (Temp, Loc))),
Typ => T));
end if;
+
+ -- Generate:
+ -- Set_Finalize_Address (<PtrT>FM, <T>FD'Unrestricted_Access);
+
+ -- Since .NET/JVM compilers do not support address arithmetic,
+ -- this call is skipped. The same is done for CodePeer because
+ -- primitive Finalize_Address is never generated.
+
+ if VM_Target = No_VM
+ and then not CodePeer_Mode
+ and then Present (Finalization_Master (PtrT))
+ then
+ Insert_Action (N,
+ Make_Set_Finalize_Address_Call
+ (Loc => Loc,
+ Typ => T,
+ Ptr_Typ => PtrT));
+ end if;
end if;
Rewrite (N, New_Reference_To (Temp, Loc));
if Ekind (PtrT) = E_Anonymous_Access_Type
and then Needs_Finalization (Dtyp)
then
- -- Anonymous access-to-controlled types allocate on the global pool
+ -- Anonymous access-to-controlled types allocate on the global pool.
+ -- Do not set this attribute on .NET/JVM since those targets do not
+ -- support pools.
- if No (Associated_Storage_Pool (PtrT)) then
+ if No (Associated_Storage_Pool (PtrT))
+ and then VM_Target = No_VM
+ then
Set_Associated_Storage_Pool (PtrT,
Get_Global_Pool_For_Access_Type (PtrT));
end if;
(Obj_Ref => New_Copy_Tree (Init_Arg1),
Typ => T));
- -- Special processing for .NET/JVM, the allocated object is
- -- attached to the finalization master. Generate:
+ if Present (Finalization_Master (PtrT)) then
- -- Attach (<PtrT>FC, Root_Controlled_Ptr (Init_Arg1));
+ -- Special processing for .NET/JVM, the allocated object
+ -- is attached to the finalization master. Generate:
- -- Types derived from [Limited_]Controlled are the only
- -- ones considered since they have fields Prev and Next.
+ -- Attach (<PtrT>FM, Root_Controlled_Ptr (Init_Arg1));
- if VM_Target /= No_VM
- and then Present (Finalization_Master (PtrT))
- and then Is_Controlled (T)
- then
- Insert_Action (N,
- Make_Attach_Call
- (Obj_Ref => New_Copy_Tree (Init_Arg1),
- Ptr_Typ => PtrT));
+ -- Types derived from [Limited_]Controlled are the only
+ -- ones considered since they have fields Prev and Next.
+
+ if VM_Target /= No_VM
+ and then Is_Controlled (T)
+ then
+ Insert_Action (N,
+ Make_Attach_Call
+ (Obj_Ref => New_Copy_Tree (Init_Arg1),
+ Ptr_Typ => PtrT));
+
+ -- Default case, generate:
+
+ -- Set_Finalize_Address
+ -- (<PtrT>FM, <T>FD'Unrestricted_Access);
+
+ -- Do not generate the above for CodePeer compilations
+ -- because primitive Finalize_Address is never built.
+
+ elsif not CodePeer_Mode then
+ Insert_Action (N,
+ Make_Set_Finalize_Address_Call
+ (Loc => Loc,
+ Typ => T,
+ Ptr_Typ => PtrT));
+ end if;
end if;
end if;
(Func_Call, Function_Id, Return_Object => Empty);
end if;
+ -- If the build-in-place function call returns a controlled object,
+ -- the finalization master will require a reference to routine
+ -- Finalize_Address of the designated type. Setting this attribute
+ -- is done in the same manner to expansion of allocators.
+
+ if Needs_Finalization (Result_Subt) then
+
+ -- Controlled types with supressed finalization do not need to
+ -- associate the address of their Finalize_Address primitives with
+ -- a master since they do not need a master to begin with.
+
+ if Is_Library_Level_Entity (Acc_Type)
+ and then Finalize_Storage_Only (Result_Subt)
+ then
+ null;
+
+ -- Do not generate the call to Make_Set_Finalize_Address for
+ -- CodePeer compilations because Finalize_Address is never built.
+
+ elsif not CodePeer_Mode then
+ Insert_Action (Allocator,
+ Make_Set_Finalize_Address_Call (Loc,
+ Typ => Etype (Function_Id),
+ Ptr_Typ => Acc_Type));
+ end if;
+ end if;
+
-- Finally, replace the allocator node with a reference to the result
-- of the function call itself (which will effectively be an access
-- to the object created by the allocator).
Statements => Make_Deep_Record_Body (Finalize_Case, Typ, True)));
end Make_Local_Deep_Finalize;
+ ------------------------------------
+ -- Make_Set_Finalize_Address_Call --
+ ------------------------------------
+
+ function Make_Set_Finalize_Address_Call
+ (Loc : Source_Ptr;
+ Typ : Entity_Id;
+ Ptr_Typ : Entity_Id) return Node_Id
+ is
+ Desig_Typ : constant Entity_Id :=
+ Available_View (Designated_Type (Ptr_Typ));
+ Utyp : Entity_Id;
+
+ begin
+ -- If the context is a class-wide allocator, we use the class-wide type
+ -- to obtain the proper Finalize_Address routine.
+
+ if Is_Class_Wide_Type (Desig_Typ) then
+ Utyp := Desig_Typ;
+
+ else
+ Utyp := Typ;
+
+ if Is_Private_Type (Utyp) and then Present (Full_View (Utyp)) then
+ Utyp := Full_View (Utyp);
+ end if;
+
+ if Is_Concurrent_Type (Utyp) then
+ Utyp := Corresponding_Record_Type (Utyp);
+ end if;
+ end if;
+
+ Utyp := Underlying_Type (Base_Type (Utyp));
+
+ -- Deal with non-tagged derivation of private views. If the parent is
+ -- now known to be protected, the finalization routine is the one
+ -- defined on the corresponding record of the ancestor (corresponding
+ -- records do not automatically inherit operations, but maybe they
+ -- should???)
+
+ if Is_Untagged_Derivation (Typ) then
+ if Is_Protected_Type (Typ) then
+ Utyp := Corresponding_Record_Type (Root_Type (Base_Type (Typ)));
+ else
+ Utyp := Underlying_Type (Root_Type (Base_Type (Typ)));
+
+ if Is_Protected_Type (Utyp) then
+ Utyp := Corresponding_Record_Type (Utyp);
+ end if;
+ end if;
+ end if;
+
+ -- If the underlying_type is a subtype, we are dealing with the
+ -- completion of a private type. We need to access the base type and
+ -- generate a conversion to it.
+
+ if Utyp /= Base_Type (Utyp) then
+ pragma Assert (Is_Private_Type (Typ));
+
+ Utyp := Base_Type (Utyp);
+ end if;
+
+ -- Generate:
+ -- Set_Finalize_Address (<Ptr_Typ>FM, <Utyp>FD'Unrestricted_Access);
+
+ return
+ Make_Procedure_Call_Statement (Loc,
+ Name =>
+ New_Reference_To (RTE (RE_Set_Finalize_Address), Loc),
+ Parameter_Associations => New_List (
+ New_Reference_To (Finalization_Master (Ptr_Typ), Loc),
+ Make_Attribute_Reference (Loc,
+ Prefix =>
+ New_Reference_To (TSS (Utyp, TSS_Finalize_Address), Loc),
+ Attribute_Name => Name_Unrestricted_Access)));
+ end Make_Set_Finalize_Address_Call;
+
--------------------------
-- Make_Transient_Block --
--------------------------
-- Create a special version of Deep_Finalize with identifier Nam. The
-- routine has state information and can parform partial finalization.
+ function Make_Set_Finalize_Address_Call
+ (Loc : Source_Ptr;
+ Typ : Entity_Id;
+ Ptr_Typ : Entity_Id) return Node_Id;
+ -- Generate the following call:
+ --
+ -- Set_Finalize_Address (<Ptr_Typ>FM, <Typ>FD'Unrestricted_Access);
+ --
+ -- where Finalize_Address is the corresponding TSS primitive of type Typ
+ -- and Ptr_Typ is the access type of the related allocation. Loc is the
+ -- source location of the related allocator.
+
--------------------------------------------
-- Task and Protected Object finalization --
--------------------------------------------
@item @b{ppc-aix}
@item @code{@ @ }@i{rts-native (default)}
@item @code{@ @ @ @ }Tasking @tab native AIX threads
+@item @code{@ @ @ @ }Exceptions @tab ZCX
+@*
+@item @code{@ @ }@i{rts-sjlj}
+@item @code{@ @ @ @ }Tasking @tab native AIX threads
@item @code{@ @ @ @ }Exceptions @tab SJLJ
@*
@item @b{ppc-darwin}
@item @b{x86-solaris}
@item @code{@ @ }@i{rts-native (default)}
@item @code{@ @ @ @ }Tasking @tab native Solaris threads
+@item @code{@ @ @ @ }Exceptions @tab ZCX
+@*
+@item @code{@ @ }@i{rts-sjlj}
+@item @code{@ @ @ @ }Tasking @tab native Solaris threads library
@item @code{@ @ @ @ }Exceptions @tab SJLJ
@*
@item @b{x86-windows}
pragma Assert (SCO_Table.Table (Start).C1 = 's');
end loop;
- Write_Info_Terminate;
+ if Ctr > 0 then
+ Write_Info_Terminate;
+ end if;
-- Statement continuations should not occur since they
-- are supposed to have been handled in the loop above.
typedef struct
{
_Unwind_Action phase;
- char * description;
+ const char * description;
} phase_descriptor;
static const phase_descriptor phase_descriptors[]
} region_descriptor;
-static void
-db_region_for (region_descriptor *region, _Unwind_Context *uw_context)
+/* Extract and adjust the IP (instruction pointer) from an exception
+ context. */
+
+static _Unwind_Ptr
+get_ip_from_context (_Unwind_Context *uw_context)
{
int ip_before_insn = 0;
#ifdef HAVE_GETIPINFO
#else
_Unwind_Ptr ip = _Unwind_GetIP (uw_context);
#endif
+ /* Subtract 1 if necessary because GetIPInfo yields a call return address
+ in this case, while we are interested in information for the call point.
+ This does not always yield the exact call instruction address but always
+ brings the IP back within the corresponding region. */
if (!ip_before_insn)
ip--;
+ return ip;
+}
+
+static void
+db_region_for (region_descriptor *region, _Unwind_Context *uw_context)
+{
+ _Unwind_Ptr ip;
+
if (! (db_accepted_codes () & DB_REGIONS))
return;
+ ip = get_ip_from_context (uw_context);
+
db (DB_REGIONS, "For ip @ 0x%08x => ", ip);
if (region->lsda)
static void
db_action_for (action_descriptor *action, _Unwind_Context *uw_context)
{
- int ip_before_insn = 0;
-#ifdef HAVE_GETIPINFO
- _Unwind_Ptr ip = _Unwind_GetIPInfo (uw_context, &ip_before_insn);
-#else
- _Unwind_Ptr ip = _Unwind_GetIP (uw_context);
-#endif
- if (!ip_before_insn)
- ip--;
+ _Unwind_Ptr ip = get_ip_from_context (uw_context);
db (DB_ACTIONS, "For ip @ 0x%08x => ", ip);
region_descriptor *region,
action_descriptor *action)
{
- int ip_before_insn = 0;
-#ifdef HAVE_GETIPINFO
- _Unwind_Ptr call_site = _Unwind_GetIPInfo (uw_context, &ip_before_insn);
-#else
- _Unwind_Ptr call_site = _Unwind_GetIP (uw_context);
-#endif
- /* Subtract 1 if necessary because GetIPInfo returns the actual call site
- value + 1 in this case. */
- if (!ip_before_insn)
- call_site--;
+ _Unwind_Ptr call_site = get_ip_from_context (uw_context);
/* call_site is a direct index into the call-site table, with two special
values : -1 for no-action and 0 for "terminate". The latter should never
action_descriptor *action)
{
const unsigned char *p = region->call_site_table;
- int ip_before_insn = 0;
-#ifdef HAVE_GETIPINFO
- _Unwind_Ptr ip = _Unwind_GetIPInfo (uw_context, &ip_before_insn);
-#else
- _Unwind_Ptr ip = _Unwind_GetIP (uw_context);
-#endif
- /* Subtract 1 if necessary because GetIPInfo yields a call return address
- in this case, while we are interested in information for the call point.
- This does not always yield the exact call instruction address but always
- brings the IP back within the corresponding region. */
- if (!ip_before_insn)
- ip--;
+ _Unwind_Ptr ip = get_ip_from_context (uw_context);
/* Unless we are able to determine otherwise... */
action->kind = nothing;
RE_Finalization_Master, -- System.Finalization_Masters
RE_Finalization_Master_Ptr, -- System.Finalization_Masters
RE_Set_Base_Pool, -- System.Finalization_Masters
+ RE_Set_Finalize_Address, -- System.Finalization_Masters
RE_Root_Controlled, -- System.Finalization_Root
RE_Root_Controlled_Ptr, -- System.Finalization_Root
RE_Finalization_Master => System_Finalization_Masters,
RE_Finalization_Master_Ptr => System_Finalization_Masters,
RE_Set_Base_Pool => System_Finalization_Masters,
+ RE_Set_Finalize_Address => System_Finalization_Masters,
RE_Root_Controlled => System_Finalization_Root,
RE_Root_Controlled_Ptr => System_Finalization_Root,
------------------------------------------------------------------------------
with Ada.Exceptions; use Ada.Exceptions;
-
+with System.Address_Image;
+with System.IO; use System.IO;
with System.Soft_Links; use System.Soft_Links;
with System.Storage_Elements; use System.Storage_Elements;
procedure Detach (N : not null FM_Node_Ptr) is
begin
- -- N must be attached to some list
-
- pragma Assert (N.Next /= null and then N.Prev /= null);
-
- Lock_Task.all;
+ if N.Prev /= null and then N.Next /= null then
+ Lock_Task.all;
- N.Prev.Next := N.Next;
- N.Next.Prev := N.Prev;
+ N.Prev.Next := N.Next;
+ N.Next.Prev := N.Prev;
+ N.Prev := null;
+ N.Next := null;
- Unlock_Task.all;
+ Unlock_Task.all;
+ end if;
-- Note: No need to unlock in case of an exception because the above
-- code can never raise one.
Obj_Addr : Address;
Raised : Boolean := False;
+ function Is_Empty_List (L : not null FM_Node_Ptr) return Boolean;
+ -- Determine whether a list contains only one element, the dummy head
+
+ -------------------
+ -- Is_Empty_List --
+ -------------------
+
+ function Is_Empty_List (L : not null FM_Node_Ptr) return Boolean is
+ begin
+ return L.Next = L and then L.Prev = L;
+ end Is_Empty_List;
+
+ -- Start of processing for Finalize
+
begin
-- It is possible for multiple tasks to cause the finalization of the
-- same master. Let only one task finalize the objects.
Master.Finalization_Started := True;
- -- Skip the dummy head
+ while not Is_Empty_List (Master.Objects'Unchecked_Access) loop
+ Curr_Ptr := Master.Objects.Next;
- Curr_Ptr := Master.Objects.Next;
- while Curr_Ptr /= Master.Objects'Unchecked_Access loop
+ Detach (Curr_Ptr);
- -- If primitive Finalize_Address is not set, then the expansion of
- -- the designated type or that of the allocator failed. This is a
- -- serious error.
-
- if Master.Finalize_Address = null then
- raise Program_Error
- with "primitive Finalize_Address not available";
- end if;
+ if Master.Finalize_Address /= null then
- -- Skip the list header in order to offer proper object layout for
- -- finalization and call Finalize_Address.
+ -- Skip the list header in order to offer proper object layout for
+ -- finalization and call Finalize_Address.
- Obj_Addr := Curr_Ptr.all'Address + Header_Offset;
+ Obj_Addr := Curr_Ptr.all'Address + Header_Offset;
- begin
- Master.Finalize_Address (Obj_Addr);
+ begin
+ Master.Finalize_Address (Obj_Addr);
- exception
- when Fin_Occur : others =>
- if not Raised then
- Raised := True;
- Save_Occurrence (Ex_Occur, Fin_Occur);
- end if;
- end;
-
- Curr_Ptr := Curr_Ptr.Next;
+ exception
+ when Fin_Occur : others =>
+ if not Raised then
+ Raised := True;
+ Save_Occurrence (Ex_Occur, Fin_Occur);
+ end if;
+ end;
+ end if;
end loop;
-- If the finalization of a particular object failed or Finalize_Address
Master.Objects.Prev := Master.Objects'Unchecked_Access;
end Initialize;
+ --------
+ -- pm --
+ --------
+
+ procedure pm (Master : Finalization_Master) is
+ Head : constant FM_Node_Ptr := Master.Objects'Unrestricted_Access;
+ Head_Seen : Boolean := False;
+ N_Ptr : FM_Node_Ptr;
+
+ begin
+ -- Output the basic contents of a master
+
+ -- Master : 0x123456789
+ -- Base_Pool: null <or> 0x123456789
+ -- Fin_Addr : null <or> 0x123456789
+ -- Fin_Start: TRUE <or> FALSE
+
+ Put ("Master : ");
+ Put_Line (Address_Image (Master'Address));
+
+ Put ("Base_Pool: ");
+
+ if Master.Base_Pool = null then
+ Put_Line (" null");
+ else
+ Put_Line (Address_Image (Master.Base_Pool'Address));
+ end if;
+
+ Put ("Fin_Addr : ");
+
+ if Master.Finalize_Address = null then
+ Put_Line ("null");
+ else
+ Put_Line (Address_Image (Master.Finalize_Address'Address));
+ end if;
+
+ Put ("Fin_Start: ");
+ Put_Line (Master.Finalization_Started'Img);
+
+ -- Output all chained elements. The format is the following:
+
+ -- ^ <or> ? <or> null
+ -- |Header: 0x123456789 (dummy head)
+ -- | Prev: 0x123456789
+ -- | Next: 0x123456789
+ -- V
+
+ -- ^ - the current element points back to the correct element
+ -- ? - the current element points back to an erroneous element
+ -- n - the current element points back to null
+
+ -- Header - the address of the list header
+ -- Prev - the address of the list header which the current element
+ -- - points back to
+ -- Next - the address of the list header which the current element
+ -- - points to
+ -- (dummy head) - present if dummy head
+
+ N_Ptr := Head;
+ while N_Ptr /= null loop -- Should never be null; we being defensive
+ Put_Line ("V");
+
+ -- We see the head initially; we want to exit when we see the head a
+ -- SECOND time.
+
+ if N_Ptr = Head then
+ exit when Head_Seen;
+
+ Head_Seen := True;
+ end if;
+
+ -- The current element is null. This should never happen since the
+ -- list is circular.
+
+ if N_Ptr.Prev = null then
+ Put_Line ("null (ERROR)");
+
+ -- The current element points back to the correct element
+
+ elsif N_Ptr.Prev.Next = N_Ptr then
+ Put_Line ("^");
+
+ -- The current element points to an erroneous element
+
+ else
+ Put_Line ("? (ERROR)");
+ end if;
+
+ -- Output the header and fields
+
+ Put ("|Header: ");
+ Put (Address_Image (N_Ptr.all'Address));
+
+ -- Detect the dummy head
+
+ if N_Ptr = Head then
+ Put_Line (" (dummy head)");
+ else
+ Put_Line ("");
+ end if;
+
+ Put ("| Prev: ");
+
+ if N_Ptr.Prev = null then
+ Put_Line ("null");
+ else
+ Put_Line (Address_Image (N_Ptr.Prev.all'Address));
+ end if;
+
+ Put ("| Next: ");
+
+ if N_Ptr.Next = null then
+ Put_Line ("null");
+ else
+ Put_Line (Address_Image (N_Ptr.Next.all'Address));
+ end if;
+
+ N_Ptr := N_Ptr.Next;
+ end loop;
+ end pm;
+
-------------------
-- Set_Base_Pool --
-------------------
Master.Base_Pool := Pool_Ptr;
end Set_Base_Pool;
+ --------------------------
+ -- Set_Finalize_Address --
+ --------------------------
+
+ procedure Set_Finalize_Address
+ (Master : in out Finalization_Master;
+ Fin_Addr_Ptr : Finalize_Address_Ptr)
+ is
+ begin
+ Master.Finalize_Address := Fin_Addr_Ptr;
+ end Set_Finalize_Address;
+
end System.Finalization_Masters;
overriding procedure Initialize (Master : in out Finalization_Master);
-- Initialize the dummy head of a finalization master
+ procedure pm (Master : Finalization_Master);
+ -- Debug routine, outputs the contents of a master
+
procedure Set_Base_Pool
(Master : in out Finalization_Master;
Pool_Ptr : Any_Storage_Pool_Ptr);
-- Set the underlying pool of a finalization master
+ procedure Set_Finalize_Address
+ (Master : in out Finalization_Master;
+ Fin_Addr_Ptr : Finalize_Address_Ptr);
+ -- Set the clean up routine of a finalization master
+
end System.Finalization_Masters;
-- --
-- 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 --
-- terms of the GNU General Public License as published by the Free Soft- --
First : Integer) return Integer;
pragma Inline_Always (Check_Unit_Last);
+ --------------
+ -- Diagonal --
+ --------------
+
+ function Diagonal (A : Matrix) return Vector is
+
+ N : constant Natural := Natural'Min (A'Length (1), A'Length (2));
+ R : Vector (A'First (1) .. A'First (1) + N - 1);
+
+ begin
+ for J in 0 .. N - 1 loop
+ R (R'First + J) := A (A'First (1) + J, A'First (2) + J);
+ end loop;
+
+ return R;
+ end Diagonal;
+
+ --------------------------
+ -- Square_Matrix_Length --
+ --------------------------
+
function Square_Matrix_Length (A : Matrix) return Natural is
begin
if A'Length (1) /= A'Length (2) then
return First + (Order - 1);
end Check_Unit_Last;
+ ---------------------
+ -- Back_Substitute --
+ ---------------------
+
+ procedure Back_Substitute (M, N : in out Matrix) is
+ pragma Assert (M'First (1) = N'First (1) and then
+ M'Last (1) = N'Last (1));
+ Max_Col : Integer := M'Last (2);
+
+ procedure Sub_Row
+ (M : in out Matrix;
+ Target : Integer;
+ Source : Integer;
+ Factor : Scalar);
+
+ procedure Sub_Row
+ (M : in out Matrix;
+ Target : Integer;
+ Source : Integer;
+ Factor : Scalar) is
+ begin
+ for J in M'Range (2) loop
+ M (Target, J) := M (Target, J) - Factor * M (Source, J);
+ end loop;
+ end Sub_Row;
+
+ begin
+ for Row in reverse M'Range (1) loop
+ Find_Non_Zero : for Col in M'First (2) .. Max_Col loop
+ if Is_Non_Zero (M (Row, Col)) then
+ -- Found first non-zero element, so subtract a multiple
+ -- of this row from all higher rows, to reduce all other
+ -- elements in this column to zero.
+
+ for J in M'First (1) .. Row - 1 loop
+ Sub_Row (N, J, Row, (M (J, Col) / M (Row, Col)));
+ Sub_Row (M, J, Row, (M (J, Col) / M (Row, Col)));
+ end loop;
+
+ Max_Col := Col - 1;
+ exit Find_Non_Zero;
+ end if;
+ end loop Find_Non_Zero;
+ end loop;
+ end Back_Substitute;
+
+ -----------------------
+ -- Forward_Eliminate --
+ -----------------------
+
+ procedure Forward_Eliminate
+ (M : in out Matrix;
+ N : in out Matrix;
+ Det : out Scalar)
+ is
+ pragma Assert (M'First (1) = N'First (1) and then
+ M'Last (1) = N'Last (1));
+
+ function "abs" (X : Scalar) return Scalar is
+ (if X < Zero then Zero - X else X);
+
+ procedure Sub_Row
+ (M : in out Matrix;
+ Target : Integer;
+ Source : Integer;
+ Factor : Scalar);
+
+ procedure Divide_Row
+ (M, N : in out Matrix;
+ Row : Integer;
+ Scale : Scalar);
+
+ procedure Switch_Row
+ (M, N : in out Matrix;
+ Row_1 : Integer;
+ Row_2 : Integer);
+
+ -------------
+ -- Sub_Row --
+ -------------
+
+ procedure Sub_Row
+ (M : in out Matrix;
+ Target : Integer;
+ Source : Integer;
+ Factor : Scalar) is
+ begin
+ for J in M'Range (2) loop
+ M (Target, J) := M (Target, J) - Factor * M (Source, J);
+ end loop;
+ end Sub_Row;
+
+ ----------------
+ -- Divide_Row --
+ ----------------
+
+ procedure Divide_Row
+ (M, N : in out Matrix;
+ Row : Integer;
+ Scale : Scalar)
+ is
+ begin
+ Det := Det * Scale;
+
+ for J in M'Range (2) loop
+ M (Row, J) := M (Row, J) / Scale;
+ end loop;
+
+ for J in N'Range (2) loop
+ N (Row - M'First (1) + N'First (1), J)
+ := N (Row - M'First (1) + N'First (1), J) / Scale;
+ end loop;
+ end Divide_Row;
+
+ ----------------
+ -- Switch_Row --
+ ----------------
+
+ procedure Switch_Row
+ (M, N : in out Matrix;
+ Row_1 : Integer;
+ Row_2 : Integer)
+ is
+ procedure Swap (X, Y : in out Scalar);
+ -- Exchange the values of X and Y
+
+ procedure Swap (X, Y : in out Scalar) is
+ T : constant Scalar := X;
+ begin
+ X := Y;
+ Y := T;
+ end Swap;
+
+ begin
+ if Row_1 /= Row_2 then
+ Det := Zero - Det;
+
+ for J in M'Range (2) loop
+ Swap (M (Row_1, J), M (Row_2, J));
+ end loop;
+
+ for J in N'Range (2) loop
+ Swap (N (Row_1 - M'First (1) + N'First (1), J),
+ N (Row_2 - M'First (1) + N'First (1), J));
+ end loop;
+ end if;
+ end Switch_Row;
+
+ I : Integer := M'First (1);
+
+ begin -- Forward_Eliminate
+ Det := One;
+
+ for J in M'Range (2) loop
+ declare
+ Max_I : Integer := I;
+ Max_Abs : Scalar := Zero;
+ begin
+ -- Find best pivot in column J, starting in row I.
+ for K in I .. M'Last (1) loop
+ declare
+ New_Abs : constant Scalar := abs M (K, J);
+ begin
+ if Max_Abs < New_Abs then
+ Max_Abs := New_Abs;
+ Max_I := K;
+ end if;
+ end;
+ end loop;
+
+ if Zero < Max_Abs then
+ Switch_Row (M, N, I, Max_I);
+ Divide_Row (M, N, I, M (I, J));
+
+ for U in I + 1 .. M'Last (1) loop
+ Sub_Row (N, U, I, M (U, J));
+ Sub_Row (M, U, I, M (U, J));
+ end loop;
+
+ exit when I >= M'Last (1);
+
+ I := I + 1;
+
+ else
+ Det := Zero; -- Zero, but we don't have literals
+ end if;
+ end;
+ end loop;
+ end Forward_Eliminate;
+
-------------------
-- Inner_Product --
-------------------
return R;
end Inner_Product;
+ -------------
+ -- L2_Norm --
+ -------------
+
+ function L2_Norm (X : Vector) return Scalar is
+ begin
+ return Sqrt (Inner_Product (X, X));
+ end L2_Norm;
+
----------------------------------
-- Matrix_Elementwise_Operation --
----------------------------------
return R;
end Outer_Product;
+ -----------------
+ -- Swap_Column --
+ -----------------
+
+ procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is
+ Temp : Scalar;
+ begin
+ for J in A'Range (1) loop
+ Temp := A (J, Left);
+ A (J, Left) := A (J, Right);
+ A (J, Right) := Temp;
+ end loop;
+ end Swap_Column;
+
---------------
-- Transpose --
---------------
-- --
-- S p e c --
-- --
--- 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 --
-- terms of the GNU General Public License as published by the Free Soft- --
package System.Generic_Array_Operations is
pragma Pure (Generic_Array_Operations);
+ ---------------------
+ -- Back_Substitute --
+ ---------------------
+
+ generic
+ type Scalar is private;
+ type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+ with function "-" (Left, Right : Scalar) return Scalar is <>;
+ with function "*" (Left, Right : Scalar) return Scalar is <>;
+ with function "/" (Left, Right : Scalar) return Scalar is <>;
+ with function Is_Non_Zero (X : Scalar) return Boolean is <>;
+ procedure Back_Substitute (M, N : in out Matrix);
+
+ --------------
+ -- Diagonal --
+ --------------
+
+ generic
+ type Scalar is private;
+ type Vector is array (Integer range <>) of Scalar;
+ type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+ function Diagonal (A : Matrix) return Vector;
+
+ -----------------------
+ -- Forward_Eliminate --
+ -----------------------
+
+ -- Use elementary row operations to put square matrix M in row echolon
+ -- form. Identical row operations are performed on matrix N, must have the
+ -- same number of rows as M.
+
+ generic
+ type Scalar is private;
+ type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+ with function "-" (Left, Right : Scalar) return Scalar is <>;
+ with function "*" (Left, Right : Scalar) return Scalar is <>;
+ with function "/" (Left, Right : Scalar) return Scalar is <>;
+ with function "<" (Left, Right : Scalar) return Boolean is <>;
+ Zero, One : Scalar;
+ procedure Forward_Eliminate
+ (M : in out Matrix;
+ N : in out Matrix;
+ Det : out Scalar);
+
--------------------------
-- Square_Matrix_Length --
--------------------------
(Left : Left_Vector;
Right : Right_Vector) return Result_Scalar;
+ -------------
+ -- L2_Norm --
+ -------------
+
+ generic
+ type Scalar is private;
+ type Vector is array (Integer range <>) of Scalar;
+ with function Inner_Product (Left, Right : Vector) return Scalar is <>;
+ with function Sqrt (X : Scalar) return Scalar is <>;
+ function L2_Norm (X : Vector) return Scalar;
+
-------------------
-- Outer_Product --
-------------------
(Left : Left_Matrix;
Right : Right_Matrix) return Result_Matrix;
+ -----------------
+ -- Swap_Column --
+ -----------------
+
+ generic
+ type Scalar is private;
+ type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+ procedure Swap_Column (A : in out Matrix; Left, Right : Integer);
+
---------------
-- Transpose --
---------------
package SST renames System.Secondary_Stack;
- NT_Exc_Stack : array (0 .. 8192) of aliased Character;
- for NT_Exc_Stack'Alignment use Standard'Maximum_Alignment;
- -- Allocate an exception stack for the main program to use.
- -- This is currently only used under VMS.
-
NT_TSD : TSD;
-- Note: we rely on the default initialization of NT_TSD
return NT_TSD.Current_Excep'Access;
end Get_Current_Excep_NT;
- ---------------------------
- -- Get_Exc_Stack_Addr_NT --
- ---------------------------
-
- function Get_Exc_Stack_Addr_NT return Address is
- begin
- return NT_Exc_Stack (NT_Exc_Stack'Last)'Address;
- end Get_Exc_Stack_Addr_NT;
-
- -----------------------------
- -- Get_Exc_Stack_Addr_Soft --
- -----------------------------
-
- function Get_Exc_Stack_Addr_Soft return Address is
- begin
- return Get_Exc_Stack_Addr.all;
- end Get_Exc_Stack_Addr_Soft;
-
------------------------
-- Get_GNAT_Exception --
------------------------
Get_Sec_Stack_Addr : Get_Address_Call := Get_Sec_Stack_Addr_NT'Access;
Set_Sec_Stack_Addr : Set_Address_Call := Set_Sec_Stack_Addr_NT'Access;
- function Get_Exc_Stack_Addr_NT return Address;
- Get_Exc_Stack_Addr : Get_Address_Call := Get_Exc_Stack_Addr_NT'Access;
-
function Get_Current_Excep_NT return EOA;
Get_Current_Excep : Get_EOA_Call := Get_Current_Excep_NT'Access;
pragma Inline (Get_Sec_Stack_Addr_Soft);
pragma Inline (Set_Sec_Stack_Addr_Soft);
- function Get_Exc_Stack_Addr_Soft return Address;
-
-- The following is a dummy record designed to mimic Communication_Block as
-- defined in s-tpobop.ads:
N_Ptr := Address_To_FM_Node_Ptr
(N_Addr + Header_And_Padding - Header_Offset);
- if Master.Finalize_Address = null then
- Master.Finalize_Address := Fin_Address;
- end if;
-
-- Prepend the allocated object to the finalization master
Attach (N_Ptr, Master.Objects'Unchecked_Access);
+ if Master.Finalize_Address = null then
+ Master.Finalize_Address := Fin_Address;
+ end if;
+
-- Move the address from the hidden list header to the start of the
-- object. This operation effectively hides the list header.
new Ada.Unchecked_Conversion
(Task_Id, System.Task_Primitives.Task_Address);
- function Get_Exc_Stack_Addr return Address;
- -- Replace System.Soft_Links.Get_Exc_Stack_Addr_NT
-
procedure Timer_Sleep_AST (ID : Address);
pragma Convention (C, Timer_Sleep_AST);
-- Signal the condition variable when AST fires
if Result = 0 then
Succeeded := True;
- Self_ID.Common.LL.Exc_Stack_Ptr := new Exc_Stack_T;
else
if not Single_Lock then
pragma Assert (Result = 0);
end Initialize_TCB;
- ------------------------
- -- Get_Exc_Stack_Addr --
- ------------------------
-
- function Get_Exc_Stack_Addr return Address is
- begin
- return Self.Common.LL.Exc_Stack_Ptr (Exc_Stack_T'Last)'Address;
- end Get_Exc_Stack_Addr;
-
-----------------
-- Create_Task --
-----------------
procedure Free is new
Ada.Unchecked_Deallocation (Ada_Task_Control_Block, Task_Id);
- procedure Free is new Ada.Unchecked_Deallocation
- (Exc_Stack_T, Exc_Stack_Ptr_T);
-
begin
if not Single_Lock then
Result := pthread_mutex_destroy (T.Common.LL.L'Access);
Known_Tasks (T.Known_Tasks_Index) := null;
end if;
- Free (T.Common.LL.Exc_Stack_Ptr);
Free (Tmp);
if Is_Self then
begin
Environment_Task_Id := Environment_Task;
- SSL.Get_Exc_Stack_Addr := Get_Exc_Stack_Addr'Access;
-
-- Initialize the lock used to synchronize chain of all ATCBs
Initialize_Lock (Single_RTS_Lock'Access, RTS_Lock_Level);
-- --
-- S p e c --
-- --
--- Copyright (C) 1991-2009, Free Software Foundation, Inc. --
+-- Copyright (C) 1991-2011, Free Software Foundation, Inc. --
-- --
-- GNARL 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- --
private
- type Exc_Stack_T is array (0 .. 8192) of aliased Character;
- for Exc_Stack_T'Alignment use Standard'Maximum_Alignment;
- type Exc_Stack_Ptr_T is access all Exc_Stack_T;
-
type Lock is record
L : aliased System.OS_Interface.pthread_mutex_t;
Prio : Interfaces.C.int;
L : aliased RTS_Lock;
-- Protection for all components is lock L
- Exc_Stack_Ptr : Exc_Stack_Ptr_T;
- -- ??? This needs comments
-
AST_Pending : Boolean;
-- Used to detect delay and sleep timeouts
git://gcc.gnu.org / gcc.git / commitdiff