Merge dixonsparseelim & dixondenseLU into dixonsolve

Add solveSingular & solveNonsingular to the interface of DixonSolver<..., Method::SparseElimination>

Careful : I had to remove the last argument m.trialsBeforeFailure of solve to have a unified interface for DixonSolver<..., Method::SparseElimination> and DixonSolver<...,Method::DenseElimination>

Author: Romain Lebreton

examples/Makefile.am

@@ -29,7 +29,7 @@ SUBDIRS=data
 
 EXAMPLES=rank det minpoly valence solve dot-product echelon sparseelimdet \
 sparseelimrank checksolve doubledet smithvalence charpoly blassolve \
-dixonsparseelim dixondenseLU sparsesolverat densesolverat \
+dixonsolve sparsesolverat densesolverat \
 ratdet poweroftwo_ranks power_rank genprime smithsparse nullspacebasis_rat nullspacebasis matrices
 #polysmith bench-fft bench-matpoly-mult
 # EXAMPLES+=nulp yabla 
@@ -66,8 +66,7 @@ echelon_SOURCES        = echelon.C
 smithvalence_SOURCES   = smithvalence.C
 sparseelimdet_SOURCES  = sparseelimdet.C
 sparseelimrank_SOURCES = sparseelimrank.C
-dixonsparseelim_SOURCES  = dixonsparseelim.C
-dixondenseLU_SOURCES = dixondenseLU.C
+dixonsolve_SOURCES  = dixonsolve.C
 ratdet_SOURCES       = ratdet.C
 sparsesolverat_SOURCES = sparsesolverat.C
 densesolverat_SOURCES = densesolverat.C

examples/dixondenseLU.C renamed to examples/dixonsolve.C

@@ -18,20 +18,28 @@
  * ========LICENCE========
  */
 
-/**\file examples/dixondenseelim.C
- @example examples/dixondenseelim.C
+/**\file examples/dixonsolve.C
+ @example examples/dixonsolve.C
  @author [email protected]
- * \brief Dixon System Solving Lifting using dense LU
+ @author [email protected]
+ * \brief Dixon System Solving via Lifting using dense LU or sparse LU
  * \ingroup examples
  */
 #include <iostream>
 #include <omp.h>
-#include "linbox/matrix/dense-matrix.h"
-#include "linbox/solutions/solve.h"
-#include "linbox/util/matrix-stream.h"
+
+#include "linbox/algorithms/rational-solver.h"
 #include "linbox/solutions/methods.h"
+#include "linbox/solutions/solve.h"
 #include "linbox/util/args-parser.h"
-#include "linbox/algorithms/rational-solver.h"
+#include "linbox/util/error.h"
+#include "linbox/util/matrix-stream.h"
+
+// DenseLU
+#include "linbox/matrix/dense-matrix.h"
+
+// SparseElim
+#include "linbox/matrix/sparse-matrix.h"
 
 using namespace LinBox;
 typedef Givaro::ZRing<Givaro::Integer> Ints;
@@ -49,65 +57,45 @@ struct FixPrime {
     template<class _ModField> void setBitsField() { }
 };
 
+template<typename _Matrix, typename _EliminationMethod>
+int test(_Matrix A, std::string vector_file, bool inv, bool pp, bool sparse_elim) {
+    // TODO : git rm dixondenseLU dixonsparseelim
 
-int main (int argc, char **argv) {
-
-    // TODO : seed ?
-    std::string matrix_file = "";
-    std::string vector_file = "";
-    bool inv = false;
-    bool pp = false;
+    std::cout << "A is " << A.rowdim() << " by " << A.coldim() << std::endl;
 
-    Argument as[] = {
-        { 'm', "-m FILE", "Set the input file for the matrix.",  TYPE_STR , &matrix_file },
-        { 'v', "-v FILE", "Set the input file for the vector.",  TYPE_STR , &vector_file },
-        { 'i', "-i"     , "whether the matrix is invertible.",  TYPE_BOOL , &inv },
-        { 'p', "-p"     , "whether you want to pretty print the matrix.",  TYPE_BOOL , &pp },
-        END_OF_ARGUMENTS
-    };
+    if (pp)
+    {
+            // Print Matrix
 
-    FFLAS::parseArguments(argc,argv,as);
+            // Matrix Market
+            // std::cout << "A is " << A << std::endl;
 
-    // Matrix File
-    if (matrix_file.empty()) {
-        std::cerr << "You must specify an input file for the matrix with -m" << std::endl;
-        exit(-1);
+            // Maple
+        A.write(std::cout << "A:=", Tag::FileFormat::Maple) << ';' << std::endl;
     }
-    std::ifstream input (matrix_file);
-    if (!input) { std::cerr << "Error opening matrix file " << argv[1] << std::endl; return -1; }
 
-    // Vector File 
+        // Vector File
+    Ints ZZ;
     std::ifstream invect;
+
+    ZVector B(ZZ, A.rowdim());
+
     bool createB = vector_file.empty();
     if (!createB) {
         invect.open (vector_file, std::ifstream::in);
         if (!invect) {
             createB = true;
+        } else {
+            for(ZVector::iterator it=B.begin(); it != B.end(); ++it)
+                invect >> *it;
         }
     }
 
-        // Read Integral matrix from File
-    Ints ZZ;
-    MatrixStream< Ints > ms( ZZ, input );
-    DenseMatrix<Ints> A (ms);
-    Ints::Element d;
-    std::cout << "A is " << A.rowdim() << " by " << A.coldim() << std::endl;
-
-    if (pp)
-    {
-            // Print Matrix
-
-            // Matrix Market
-            // std::cout << "A is " << A << std::endl;
-
-            // Maple
-        A.write(std::cout << "A:=", Tag::FileFormat::Maple) << ';' << std::endl;
-    }
-
         // Vectors
-    ZVector X(ZZ, A.coldim()),U(ZZ, A.coldim()),B(ZZ, A.rowdim());
+    ZVector X(ZZ, A.coldim());
 
     if (createB) {
+        ZVector U(ZZ, A.coldim());
         if (inv) {
             std::cerr << "Creating a random {-1,1} vector " << std::endl;
             srand48( BaseTimer::seed() );
@@ -119,64 +107,84 @@ int main (int argc, char **argv) {
         } else {
             std::cerr << "Creating a random consistant {-1,1} vector " << std::endl;
             srand48( BaseTimer::seed() );
-            for(auto& it:U)
+            for(FFPACK::rns_double::integer& it:U)
                 if (drand48() <0.5)
                     it = -1;
                 else
                     it = 1;
-            
+
             // B = A U
             A.apply(B,U);
         }
-    } else {
-        for(ZVector::iterator it=B.begin();
-            it != B.end(); ++it)
-            invect >> *it;
     }
 
     if(pp)
     {
             // Print RHS
         B.write(std::cout << "B:=", Tag::FileFormat::Maple) << ';' << std::endl;
     }
-    
+
     std::cout << "B is " << B.size() << "x1" << std::endl;
 
     Timer chrono;
 
         // BlasElimination
-    Method::DenseElimination M;
+    // TODO : à vérifier si cela marche avec Method::DenseElimination
+    _EliminationMethod M;
+    if (inv){
+        M.singularity = Singularity::NonSingular;
+    }
+
 
         //====================================================
         // BEGIN Replacement solve with fixed prime
+
+    Ints::Element d;
     Method::Dixon m(M);
     typedef Givaro::Modular<double> Field;
 		// 0.7213475205 is an upper approximation of 1/(2log(2))
     size_t bitsize((size_t)( 26-(int)ceil(log((double)A.rowdim())*0.7213475205)));
     Givaro::Integer randomPrime( *(PrimeIterator<>(bitsize)) );
 
     FixPrime fixedprime( randomPrime );
-    DixonSolver<Ints, Field, FixPrime, Method::DenseElimination> rsolve(A.field(), fixedprime);
+    DixonSolver<Ints, Field, FixPrime, _EliminationMethod> rsolve(A.field(), fixedprime);
     std::cout << "Using: " << *fixedprime << " as the fixed p-adic." << std::endl;
 
     chrono.start();
-    if (inv)
-    {
-        SolverReturnStatus ss = rsolve.solveNonsingular(X, d, A, B, false,(int)m.trialsBeforeFailure);
-        if (ss != SS_OK) {
-            std::cerr << "Error during solveNonsingular (possibly singular matrix or p-adic precision too small)" << std::endl;
-            exit(-1);
+    if (!sparse_elim){
+            // Dense Elimination
+        if (inv)
+        {
+            std::cout << "Solving using Dense Elimination for non singular system" << std::endl;
+            SolverReturnStatus ss = rsolve.solveNonsingular(X, d, A, B);
+            if (ss != SS_OK) {
+                std::cerr << "Error during solveNonsingular (possibly singular matrix or p-adic precision too small)" << std::endl;
+                exit(-1);
+            }
         }
-    }
-    else
-    {
-        SolverReturnStatus ss = rsolve.solve(X, d, A, B, false,(int)m.trialsBeforeFailure);
-        if (ss == SS_FAILED){
-            std::cerr << "Error during solve (all primes used were bad)" << std::endl;
-            exit(-1);
+        else
+        {
+            std::cout << "Solving using Dense Elimination for any system" << std::endl;
+            SolverReturnStatus ss = rsolve.solve(X, d, A, B);
+            if (ss == SS_FAILED){
+                std::cerr << "Error during solve (all primes used were bad)" << std::endl;
+                exit(-1);
+            }
+            if (ss == SS_INCONSISTENT){
+                std::cerr << "Error: system appeared inconsistent" << std::endl;
+                exit(-1);
+            }
         }
-        if (ss == SS_INCONSISTENT){
-            std::cerr << "Error: system appeared inconsistent" << std::endl;
+    } else {
+            // Sparse Elimination
+        try
+        {
+            std::cout << "Solving using Sparse Elimination for any system" << std::endl;
+            rsolve.solve(X, d, A, B);
+        }
+        catch(LinboxError& e)
+        {
+            std::cerr << e << '\n';
             exit(-1);
         }
     }
@@ -215,15 +223,59 @@ int main (int argc, char **argv) {
     {
             // Print Solution
 
-        std::cout << "(DenseElimination) Solution is [";
+        std::cout << "Solution is [";
         for(auto it:X) ZZ.write(std::cout, it) << " ";
         std::cout << "] / ";
         ZZ.write(std::cout, d)<< std::endl;
     }
-
     return 0;
 }
 
+int main (int argc, char **argv) {
+
+    // TODO : seed ?
+    std::string matrix_file = "";
+    std::string vector_file = "";
+    bool inv = false;
+    bool pp = false;
+    bool sparse_elim = false;
+
+
+    Argument as[] = {
+        { 'm', "-m FILE", "Set the input file for the matrix.",  TYPE_STR , &matrix_file },
+        { 'v', "-v FILE", "Set the input file for the vector.",  TYPE_STR , &vector_file },
+        { 'i', "-i"     , "whether the matrix is invertible.",  TYPE_BOOL , &inv },
+        { 'p', "-p"     , "whether you want to pretty print the matrix.",  TYPE_BOOL , &pp },
+        { 's', "-s"     , "whether to use sparse elimination.",  TYPE_BOOL , &sparse_elim },
+        END_OF_ARGUMENTS
+    };
+
+    FFLAS::parseArguments(argc,argv,as);
+
+    // Matrix File
+    if (matrix_file.empty()) {
+        std::cerr << "You must specify an input file for the matrix with -m" << std::endl;
+        exit(-1);
+    }
+    std::ifstream input (matrix_file);
+    if (!input) { std::cerr << "Error opening matrix file " << argv[1] << std::endl; exit(-1); }
+
+
+        // Read Integral matrix from File
+    Ints ZZ;
+    MatrixStream< Ints > ms( ZZ, input );
+
+    if (sparse_elim){
+        SparseMatrix<Ints> A(ms);
+        return test<SparseMatrix<Ints>,Method::SparseElimination>
+            (A, vector_file, inv , pp, sparse_elim);
+    } else {
+        DenseMatrix<Ints> A(ms);
+        return test<DenseMatrix<Ints>,Method::DenseElimination>
+            (A, vector_file, inv , pp, sparse_elim);
+    }
+}
+
 
 // Local Variables:
 // mode: C++

linbox/algorithms/rational-solver.h

@@ -716,11 +716,23 @@ namespace LinBox {
 		{}
 
 
-		// solve non singular system
+		// solve any system
 		template<class IMatrix, class Vector1, class Vector2>
 		SolverReturnStatus solve(Vector1& num, Integer& den,
                                  const IMatrix& A, const Vector2& b,
                                  int maxPrimes = DEFAULT_MAXPRIMES) const;
+
+		// solve non singular system
+		template<class IMatrix, class Vector1, class Vector2>
+		SolverReturnStatus solveNonsingular(Vector1& num, Integer& den,
+                                 const IMatrix& A, const Vector2& b,
+                                 int maxPrimes = DEFAULT_MAXPRIMES) const;
+
+		// solve singular system
+		template<class IMatrix, class Vector1, class Vector2>
+		SolverReturnStatus solveSingular(Vector1& num, Integer& den,
+                                 const IMatrix& A, const Vector2& b,
+                                 int maxPrimes = DEFAULT_MAXPRIMES) const;
 	};
 
 }

linbox/algorithms/rational-solver.inl

@@ -636,6 +636,30 @@ namespace LinBox
 			return SS_OK;
 	}
 
+	template <class Ring, class Field, class RandomPrime>
+	template <class IMatrix, class Vector1, class Vector2>
+	SolverReturnStatus
+	DixonSolver<Ring,Field,RandomPrime,Method::SparseElimination>::solveNonsingular(Vector1& num,
+                                                                          Integer& den,
+                                                                          const IMatrix& A,
+                                                                          const Vector2& b,
+                                                                          int maxPrimes) const
+	{
+		return solve(num, den, A, b, maxPrimes);
+	}
+
+	template <class Ring, class Field, class RandomPrime>
+	template <class IMatrix, class Vector1, class Vector2>
+	SolverReturnStatus
+	DixonSolver<Ring,Field,RandomPrime,Method::SparseElimination>::solveSingular(Vector1& num,
+                                                                          Integer& den,
+                                                                          const IMatrix& A,
+                                                                          const Vector2& b,
+                                                                          int maxPrimes) const
+	{
+		return solve(num, den, A, b, maxPrimes);
+	}
+
 } //end of namespace LinBox
 
 //BB : moved the following "guarded" code in a new file, verbatim :

linbox/vector/vector-domain.h

@@ -82,7 +82,7 @@ namespace LinBox
 		VectorDomainBase () : _faxpy(nullptr) {std::cout<<"CST DEFAULT VDB"<<std::endl;}
 		VectorDomainBase (const Field &F) :  _faxpy(new FieldAXPY<Field>(F))
                 { /*std::cerr << "VDD cstor " << this << std::endl;*/}
-                ~VectorDomainBase()	{ if (_faxpy != nullptr) delete _faxpy; }
+		~VectorDomainBase()	{ if (_faxpy != nullptr) delete _faxpy; }
 
 		void init(const Field &F) { if (_faxpy != nullptr) delete _faxpy; _faxpy = new FieldAXPY<Field>(F); }
 		inline const Field & field() const { return _faxpy->field(); }