From patchwork Wed Apr 27 13:05:55 2016 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Arnaud Charlet X-Patchwork-Id: 615631 Return-Path: X-Original-To: incoming@patchwork.ozlabs.org Delivered-To: patchwork-incoming@bilbo.ozlabs.org Received: from sourceware.org (server1.sourceware.org [209.132.180.131]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (No client certificate requested) by ozlabs.org (Postfix) with ESMTPS id 3qw0Zr0kqTz9t4c for ; Wed, 27 Apr 2016 23:06:09 +1000 (AEST) Authentication-Results: ozlabs.org; dkim=pass (1024-bit key; unprotected) header.d=gcc.gnu.org header.i=@gcc.gnu.org header.b=cVTsRlHZ; dkim-atps=neutral DomainKey-Signature: a=rsa-sha1; c=nofws; d=gcc.gnu.org; h=list-id :list-unsubscribe:list-archive:list-post:list-help:sender:date :from:to:cc:subject:message-id:mime-version:content-type; q=dns; s=default; b=JhYo6bycaKnxRDNUjGLwMWsRyFXVI/sq0WzQB8duzBeeRAJWVc PIY8N3OIDBVw8tv/yX6Lw1wcrQUSyoG6jWQ+oxPEVeTC1Ufb8e2hec3BDdQBawEU dR4E1rQ3CP47NjVMqa0EGmE92CR6kvJxB5dOE/C4YSpBtK5q6HkC8yeOM= DKIM-Signature: v=1; a=rsa-sha1; c=relaxed; d=gcc.gnu.org; h=list-id :list-unsubscribe:list-archive:list-post:list-help:sender:date :from:to:cc:subject:message-id:mime-version:content-type; s= default; bh=gLc4EPodgBBQzllDtRqKEHoLVgw=; b=cVTsRlHZW74nXJWDH+uM k5F9IH6CUmq9RvIOBsJmLNf+GOQDBLGpEPsheIDJvEJ2VWpof3MbXrOPWp3nUBpU 9PJxnFFOC710TDuKGKN3degeUqnsZaCBKt1gGZhZC9u8sfDCTf5pPPHe8EqEKVjN 1YPQohrmf0H6HpJSJtYNtH8= Received: (qmail 8645 invoked by alias); 27 Apr 2016 13:06:01 -0000 Mailing-List: contact gcc-patches-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Archive: List-Post: List-Help: Sender: gcc-patches-owner@gcc.gnu.org Delivered-To: mailing list gcc-patches@gcc.gnu.org Received: (qmail 8623 invoked by uid 89); 27 Apr 2016 13:06:01 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=2.2 required=5.0 tests=AWL, BAYES_50, KAM_ASCII_DIVIDERS, KAM_LAZY_DOMAIN_SECURITY, RCVD_IN_DNSWL_NONE autolearn=no version=3.3.2 spammy=Solution, 235481, Schonberg, schonbergadacorecom X-HELO: rock.gnat.com Received: from rock.gnat.com (HELO rock.gnat.com) (205.232.38.15) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with (AES256-SHA encrypted) ESMTPS; Wed, 27 Apr 2016 13:05:57 +0000 Received: from localhost (localhost.localdomain [127.0.0.1]) by filtered-rock.gnat.com (Postfix) with ESMTP id 7FD1E116CCA; Wed, 27 Apr 2016 09:05:55 -0400 (EDT) Received: from rock.gnat.com ([127.0.0.1]) by localhost (rock.gnat.com [127.0.0.1]) (amavisd-new, port 10024) with LMTP id YTRRzYhjahW4; Wed, 27 Apr 2016 09:05:55 -0400 (EDT) Received: from tron.gnat.com (tron.gnat.com [205.232.38.10]) by rock.gnat.com (Postfix) with ESMTP id 6F3F1116684; Wed, 27 Apr 2016 09:05:55 -0400 (EDT) Received: by tron.gnat.com (Postfix, from userid 4192) id 6C956370; Wed, 27 Apr 2016 09:05:55 -0400 (EDT) Date: Wed, 27 Apr 2016 09:05:55 -0400 From: Arnaud Charlet To: gcc-patches@gcc.gnu.org Cc: Ed Schonberg Subject: [Ada] Detect singular matrices in Solve primitives for vectors and matrices. Message-ID: <20160427130555.GA79191@adacore.com> MIME-Version: 1.0 Content-Disposition: inline User-Agent: Mutt/1.5.23 (2014-03-12) This patch detects cases when the determinant of a system of equations is exactly Zero, and raises the proper exception, as mandated by RN G.3.2 (68/2) and G.3.2 (89/2). Executing the following: gnatmake -q debug_solve.adb debug_solve must yield: System is: 1.00000E+00 x_1 + 0.00000E+00 x_2 + 1.00000E+00 x_3 = 1.00000E+00 2.00000E+00 x_1 + 0.00000E+00 x_2 + 2.00000E+00 x_3 = 2.00000E+00 3.00000E+00 x_1 + 0.00000E+00 x_2 + 1.00000E+00 x_3 = 2.00000E+00 System is: raised CONSTRAINT_ERROR : debug_solve.float_arrays.Instantiations.Solve: matrix is singular --- with Ada.Numerics.Generic_Real_Arrays; with Ada.Text_IO; procedure debug_solve is package float_arrays is new Ada.Numerics.Generic_Real_Arrays (Real => Float); matrix : constant float_arrays.Real_Matrix (1 .. 3, 1 .. 3) := ((1.0, 0.0, 1.0), (2.0, 0.0, 2.0), (3.0, 0.0, 1.0)); rhs : constant float_arrays.Real_Vector (1 .. 3) := (1.0, 2.0, 2.0); solution : float_arrays.Real_Vector (1 .. 3); begin Ada.Text_IO.Put_Line ("System is:"); Ada.Text_IO.Put_Line (Float'Image (matrix (1, 1)) & " x_1 + " & Float'Image (matrix (1, 2)) & " x_2 + " & Float'Image (matrix (1, 3)) & " x_3 = " & Float'Image(rhs(1))); Ada.Text_IO.Put_Line (Float'Image (matrix (2, 1)) & " x_1 + " & Float'Image (matrix (2, 2)) & " x_2 + " & Float'Image (matrix (2, 3)) & " x_3 = " & Float'Image(rhs(2))); Ada.Text_IO.Put_Line (Float'Image (matrix (3, 1)) & " x_1 + " & Float'Image (matrix (3, 2)) & " x_2 + " & Float'Image (matrix (3, 3)) & " x_3 = " & Float'Image(rhs(3))); Ada.Text_IO.Put_Line ("System is:"); solution := float_arrays.Solve (A => matrix, X => rhs); Ada.Text_IO.Put_Line ("Solution is x = (" & Float'Image (solution (1)) & ", " & Float'Image (solution (2)) & ", " & Float'Image (solution (3)) & ")"); end; Tested on x86_64-pc-linux-gnu, committed on trunk 2016-04-27 Ed Schonberg * s-gearop.ads (Matrix_Vector_Solution, Matrix_Matrix_Solution): Add scalar formal object Zero, to allow detection and report when the matrix is singular. * s-gearop.adb (Matrix_Vector_Solution, Matrix_Matrix_Solution): Raise Constraint_Error if the Forward_Eliminate pass has determined that determinant is Zero.o * s-ngrear.adb (Solve): Add actual for Zero in corresponding instantiations. * s-ngcoar.adb (Solve): Ditto. Index: a-ngrear.adb =================================================================== --- a-ngrear.adb (revision 235481) +++ a-ngrear.adb (working copy) @@ -6,7 +6,7 @@ -- -- -- B o d y -- -- -- --- Copyright (C) 2006-2012, Free Software Foundation, Inc. -- +-- Copyright (C) 2006-2016, 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- -- @@ -337,10 +337,11 @@ Result_Matrix => Real_Matrix, Operation => "abs"); - function Solve is - new Matrix_Vector_Solution (Real'Base, Real_Vector, Real_Matrix); + function Solve is new + Matrix_Vector_Solution (Real'Base, 0.0, Real_Vector, Real_Matrix); - function Solve is new Matrix_Matrix_Solution (Real'Base, Real_Matrix); + function Solve is new + Matrix_Matrix_Solution (Real'Base, 0.0, Real_Matrix); function Unit_Matrix is new Generic_Array_Operations.Unit_Matrix Index: s-gearop.adb =================================================================== --- s-gearop.adb (revision 235481) +++ s-gearop.adb (working copy) @@ -6,7 +6,7 @@ -- -- -- B o d y -- -- -- --- Copyright (C) 2006-2012, Free Software Foundation, Inc. -- +-- Copyright (C) 2006-2016, Free Software Foundation, Inc. -- -- -- -- GNAT is free software; you can redistribute it and/or modify it under -- -- terms of the GNU General Public License as published by the Free Soft- -- @@ -30,9 +30,7 @@ ------------------------------------------------------------------------------ with Ada.Numerics; use Ada.Numerics; - package body System.Generic_Array_Operations is - function Check_Unit_Last (Index : Integer; Order : Positive; @@ -696,6 +694,11 @@ end loop; Forward_Eliminate (MA, MX, Det); + + if Det = Zero then + raise Constraint_Error with "matrix is singular"; + end if; + Back_Substitute (MA, MX); for J in 0 .. R'Length - 1 loop @@ -735,6 +738,11 @@ end loop; Forward_Eliminate (MA, MB, Det); + + if Det = Zero then + raise Constraint_Error with "matrix is singular"; + end if; + Back_Substitute (MA, MB); return MB; Index: s-gearop.ads =================================================================== --- s-gearop.ads (revision 235481) +++ s-gearop.ads (working copy) @@ -6,7 +6,7 @@ -- -- -- S p e c -- -- -- --- Copyright (C) 2006-2011, Free Software Foundation, Inc. -- +-- Copyright (C) 2006-2016, 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- -- @@ -396,6 +396,7 @@ generic type Scalar is private; + Zero : Scalar; type Vector is array (Integer range <>) of Scalar; type Matrix is array (Integer range <>, Integer range <>) of Scalar; with procedure Back_Substitute (M, N : in out Matrix) is <>; @@ -411,6 +412,7 @@ generic type Scalar is private; + Zero : Scalar; type Matrix is array (Integer range <>, Integer range <>) of Scalar; with procedure Back_Substitute (M, N : in out Matrix) is <>; with procedure Forward_Eliminate Index: a-ngcoar.adb =================================================================== --- a-ngcoar.adb (revision 235481) +++ a-ngcoar.adb (working copy) @@ -6,7 +6,7 @@ -- -- -- B o d y -- -- -- --- Copyright (C) 2006-2012, Free Software Foundation, Inc. -- +-- Copyright (C) 2006-2016, Free Software Foundation, Inc. -- -- -- -- GNAT is free software; you can redistribute it and/or modify it under -- -- terms of the GNU General Public License as published by the Free Soft- -- @@ -30,7 +30,6 @@ ------------------------------------------------------------------------------ with System.Generic_Array_Operations; use System.Generic_Array_Operations; -with Ada.Numerics; use Ada.Numerics; package body Ada.Numerics.Generic_Complex_Arrays is @@ -694,11 +693,11 @@ -- Solve -- ----------- - function Solve is - new Matrix_Vector_Solution (Complex, Complex_Vector, Complex_Matrix); + function Solve is new Matrix_Vector_Solution + (Complex, (0.0, 0.0), Complex_Vector, Complex_Matrix); - function Solve is - new Matrix_Matrix_Solution (Complex, Complex_Matrix); + function Solve is new Matrix_Matrix_Solution + (Complex, (0.0, 0.0), Complex_Matrix); ----------------- -- Unit_Matrix --