[WIP] Multi modular lifting container

benchmarks/benchmark-dense-solve.C

@@ -54,8 +54,10 @@ namespace {
         int n = 500;
         int bits = 10;
         int seed = -1;
+        int primesCount = -1;
         std::string dispatchString = "Auto";
         std::string methodString = "Auto";
+        std::string rnsFgemmString = "ParallelRnsOnly";
     };
 
     template <typename Vector>
@@ -145,6 +147,8 @@ int main(int argc, char** argv)
                      {'b', "-b", "bit size", TYPE_INT, &args.bits},
                      {'s', "-s", "Seed for randomness.", TYPE_INT, &args.seed},
                      {'d', "-d", "Dispatch mode (any of: Auto, Sequential, SMP, Distributed).", TYPE_STR, &args.dispatchString},
+                     {'r', "-r", "RNS-FGEMM type (either BothParallel, BothSequential, ParallelRnsOnly or ParallelFgemmOnly).", TYPE_STR, &args.rnsFgemmString},
+                     {'p', "-p", "Enable multi-modular method, and tells how many primes to use.", TYPE_INT, &args.primesCount},
 		             {'t', "-t", "Number of threads.", TYPE_INT, &numThreads },
                      {'M', "-M",
                       "Choose the solve method (any of: Auto, Elimination, DenseElimination, SparseElimination, "
@@ -171,11 +175,26 @@ int main(int argc, char** argv)
 
     MethodBase method;
     method.pCommunicator = &communicator;
+    if (args.primesCount > 0) {
+        method.multiModularLifting = true;
+        method.primesCount = args.primesCount;
+    } else {
+        method.multiModularLifting = false;
+    }
     if (args.dispatchString == "Sequential")        method.dispatch = Dispatch::Sequential;
     else if (args.dispatchString == "SMP")          method.dispatch = Dispatch::SMP;
     else if (args.dispatchString == "Distributed")  method.dispatch = Dispatch::Distributed;
     else                                            method.dispatch = Dispatch::Auto;
 
+    if (args.rnsFgemmString == "BothParallel")              method.rnsFgemmType = RnsFgemmType::BothParallel;
+    else if (args.rnsFgemmString == "BothSequential")       method.rnsFgemmType = RnsFgemmType::BothSequential;
+    else if (args.rnsFgemmString == "ParallelRnsOnly")      method.rnsFgemmType = RnsFgemmType::ParallelRnsOnly;
+    else if (args.rnsFgemmString == "ParallelFgemmOnly")    method.rnsFgemmType = RnsFgemmType::ParallelFgemmOnly;
+    else {
+        std::cerr << "-r RNS-FGEMM type should be either BothParallel, BothSequential, ParallelRnsOnly or ParallelFgemmOnly" << std::endl;
+        return EXIT_FAILURE;
+    }
+
     // Real benchmark
 
     bool isModular = false;

linbox/algorithms/dixon-solver/dixon-solver-dense.h

@@ -89,6 +89,7 @@ namespace LinBox {
         mutable Prime _prime;
         Ring _ring;
         Field _field;
+        Method::Dixon _method;
 
         BlasMatrixDomain<Field> _bmdf;
 
@@ -113,10 +114,11 @@ namespace LinBox {
          * @param r   a Ring, set by default
          * @param rp  a RandomPrime generator, set by default
          */
-        DixonSolver(const Ring& r = Ring(), const RandomPrime& rp = RandomPrime())
+        DixonSolver(const Ring& r = Ring(), const RandomPrime& rp = RandomPrime(), const Method::Dixon& method = Method::Dixon())
             : lastCertificate(r, 0)
             , _genprime(rp)
             , _ring(r)
+            , _method(method)
         {
             _genprime.setBits(FieldTraits<Field>::bestBitSize());
             _prime = *_genprime;

linbox/algorithms/dixon-solver/dixon-solver-dense.inl

@@ -24,8 +24,10 @@
 #include "linbox/util/debug.h"
 
 #include "linbox/algorithms/lifting-container.h"
+#include "linbox/algorithms/multi-mod-lifting-container.h"
 #include "linbox/algorithms/matrix-inverse.h"
 #include "linbox/algorithms/rational-reconstruction.h"
+#include "linbox/algorithms/multi-mod-rational-reconstruction.h"
 
 namespace LinBox {
 
@@ -128,13 +130,24 @@ namespace LinBox {
             }
         } while (notfr);
 
-        typedef DixonLiftingContainer<Ring, Field, IMatrix, BlasMatrix<Field>> LiftingContainer;
-        LiftingContainer lc(_ring, *F, A, *FMP, b, _prime);
-        RationalReconstruction<LiftingContainer> re(lc);
-        if (!re.getRational(num, den, 0)) {
-            delete FMP;
-            return SS_FAILED;
+        if (_method.multiModularLifting) {
+            using LiftingContainer = MultiModLiftingContainer<Field, Ring, RandomPrime>;
+            LiftingContainer lc(_ring, _genprime, A, b, _method);
+            MultiModRationalReconstruction<LiftingContainer> re(lc);
+            if (!re.getRational(num, den)) {
+                delete FMP;
+                return SS_FAILED;
+            }
+        } else {
+            using LiftingContainer = DixonLiftingContainer<Ring, Field, IMatrix, BlasMatrix<Field>>;
+            LiftingContainer lc(_ring, *F, A, *FMP, b, _prime);
+            RationalReconstruction<LiftingContainer> re(lc);
+            if (!re.getRational(num, den, 0)) {
+                delete FMP;
+                return SS_FAILED;
+            }
         }
+
 #ifdef RSTIMING
         ttNonsingularSolve.update(re, lc);
 #endif
@@ -271,8 +284,6 @@ namespace LinBox {
     SolverReturnStatus DixonSolver<Ring, Field, RandomPrime, Method::DenseElimination>::solveApparentlyInconsistent(
         const BlasMatrix<Ring>& A, TAS& tas, BlasMatrix<Field>* Atp_minor_inv, size_t rank, const MethodBase& method)
     {
-        using LiftingContainer = DixonLiftingContainer<Ring, Field, BlasMatrix<Ring>, BlasMatrix<Field>>;
-
         if (!method.certifyInconsistency) return SS_INCONSISTENT;
 
         // @fixme Put these as class members!
@@ -295,15 +306,24 @@ namespace LinBox {
         ttCheckConsistency += tCheckConsistency;
 #endif
 
-        LiftingContainer lc(_ring, _field, At_minor, *Atp_minor_inv, zt, _prime);
-        RationalReconstruction<LiftingContainer> re(lc);
-
         BlasVector<Ring> shortNum(A.field(), rank);
         Integer shortDen;
 
-        // Dirty, but should not be called under normal circumstances
-        if (!re.getRational(shortNum, shortDen, 0)) {
-            return SS_FAILED;
+        if (_method.multiModularLifting) {
+            using LiftingContainer = MultiModLiftingContainer<Field, Ring, RandomPrime>;
+            LiftingContainer lc(_ring, _genprime, At_minor, zt, _method);
+            MultiModRationalReconstruction<LiftingContainer> re(lc);
+            if (!re.getRational(shortNum, shortDen)) {
+                return SS_FAILED;
+            }
+        }
+        else {
+            using LiftingContainer = DixonLiftingContainer<Ring, Field, BlasMatrix<Ring>, BlasMatrix<Field>>;
+            LiftingContainer lc(_ring, _field, At_minor, *Atp_minor_inv, zt, _prime);
+            RationalReconstruction<LiftingContainer> re(lc);
+            if (!re.getRational(shortNum, shortDen, 0)) {
+                return SS_FAILED;
+            }
         }
 
 #ifdef RSTIMING
@@ -699,18 +719,26 @@ namespace LinBox {
             BMDI.mulin_right(tas.Q, newb);
             newb.resize(rank);
 
-            // ----- Do lifting on sub matrix
+            // ----- Do lifting on sub matrix and reconstruct
 
             BlasMatrix<Ring> BBA_minor(A_minor);
-            LiftingContainer lc(_ring, _field, BBA_minor, *Ap_minor_inv, newb, _prime);
-
-            // ----- Reconstruct rational
-
-            RationalReconstruction<LiftingContainer> re(lc);
             VectorFraction<Ring> resultVF(_ring, rank);
-            if (!re.getRational(resultVF.numer, resultVF.denom, 0)) {
-                // dirty, but should not be called
-                return SS_FAILED;
+
+            if (_method.multiModularLifting) {
+                using LiftingContainer = MultiModLiftingContainer<Field, Ring, RandomPrime>;
+                LiftingContainer lc(_ring, _genprime, BBA_minor, newb, _method);
+                MultiModRationalReconstruction<LiftingContainer> re(lc);
+                if (!re.getRational(resultVF.numer, resultVF.denom)) {
+                    return SS_FAILED;
+                }
+            }
+            else {
+                using LiftingContainer = DixonLiftingContainer<Ring, Field, BlasMatrix<Ring>, BlasMatrix<Field>>;
+                LiftingContainer lc(_ring, _field, BBA_minor, *Ap_minor_inv, newb, _prime);
+                RationalReconstruction<LiftingContainer> re(lc);
+                if (!re.getRational(resultVF.numer, resultVF.denom, 0)) {
+                    return SS_FAILED;
+                }
             }
 
 #ifdef RSTIMING

linbox/algorithms/last-invariant-factor.h

@@ -51,7 +51,7 @@ namespace LinBox
 
 	protected:
 
-		typedef BlasVector<Ring>         DVect;
+		typedef DenseVector<Ring>         DVect;
 		Ring        r;
         mutable typename Ring::RandIter _gen;
 		Solver solver;
@@ -110,24 +110,23 @@ namespace LinBox
 			Integer r_den;
 			//std::vector<std::pair<Integer, Integer> > result (A.coldim());
 			//typename std::vector<std::pair<Integer, Integer> >::iterator result_p;
-			// vector b, RHS, 32-bit int is good enough
-			std::vector<int> b(A.rowdim());
-			typename std::vector<int>::iterator b_p;
-			typename Vector::const_iterator Prime_p;
+			DenseVector<Ring> b(r, A.rowdim());
 
 			Integer pri, quo, rem, itmp;
 
 			for (; count < threshold; ++ count) {
 				// assign b to be a random vector
-				for (b_p = b.begin(); b_p != b.end(); ++ b_p) {
-//					* b_p = rand() % 268435456 - 134217728; // may need to change to use ring's random gen.
-//					// dpritcha, 2004-07-26
-                    _gen( itmp );
-                    * b_p = (int)itmp;
+				for (auto b_p = b.begin(); b_p != b.end(); ++ b_p) {
+                   _gen( itmp );
+                        //@enhancement vector b, RHS, 32-bit is good enough
+                   * b_p = Integer((int32_t)itmp);
 				}
 
 				// try to solve Ax = b over Ring
 				tmp = solver.solveNonsingular(r_num, r_den, A, b);
+
+                // std::clog << "r_den: " << r_den << std::endl;
+
 				// If no solution found
 				if (tmp != SS_OK) {
 					r.assign (lif, r.zero);
@@ -138,7 +137,7 @@ namespace LinBox
 			}
 			// filter out primes in PRIMEL from lif.
 			if (!r. isZero (lif))
-				for ( Prime_p = PrimeL.begin();
+				for ( auto Prime_p = PrimeL.begin();
 				      Prime_p != PrimeL.end();
 				      ++ Prime_p) {
 					r.init (pri, *Prime_p);
@@ -172,21 +171,17 @@ namespace LinBox
 			Integer r1_den, r2_den;
 			//std::vector<std::pair<Integer, Integer> > result (A.coldim());
 			//typename std::vector<std::pair<Integer, Integer> >::iterator result_p;
-			// vector b, RHS, 32-bit int is good enough
-			std::vector<int> b1(A. rowdim()), b2(A. rowdim());
-			typename std::vector<int>::iterator b_p;
-			typename Vector::const_iterator Prime_p;
+			//@enhancement vector b, RHS, 32-bit instead would be good enough
+            DenseVector<Ring> b1(r, A. rowdim()), b2(r, A. rowdim());
 			Integer pri, quo, rem;
 
 			for (; count < (threshold + 1) / 2; ++ count) {
 				// assign b to be a random vector
-				for (b_p = b1. begin(); b_p != b1. end(); ++ b_p) {
-//					* b_p = rand();
-                                        *b_p = _gen.random();//(* b_p);
+				for (auto b_p = b1. begin(); b_p != b1. end(); ++ b_p) {
+                    _gen.random(*b_p);
 				}
-				for (b_p = b2. begin(); b_p != b2. end(); ++ b_p) {
-//					* b_p = rand();
-                                        *b_p = _gen.random();//(* b_p);
+				for (auto b_p = b2. begin(); b_p != b2. end(); ++ b_p) {
+                    _gen.random(*b_p);
 				}
 				// try to solve Ax = b1, b2 over Ring
 				tmp1 = solver. solveNonsingular(r1_num, r1_den, A, b1);
@@ -243,7 +238,7 @@ namespace LinBox
 
 			// filter out primes in PRIMEL from lif.
 			if (!r. isZero (lif))
-				for ( Prime_p = PrimeL.begin(); Prime_p != PrimeL.end(); ++ Prime_p) {
+				for ( auto Prime_p = PrimeL.begin(); Prime_p != PrimeL.end(); ++ Prime_p) {
 					r.init (pri, *Prime_p);
 					do {
 						r.quoRem(quo,rem,lif,pri);
@@ -253,7 +248,7 @@ namespace LinBox
 				}
 			r. gcdin (Bonus, lif);
 			if (!r. isZero (Bonus))
-				for ( Prime_p = PrimeL.begin(); Prime_p != PrimeL.end(); ++ Prime_p) {
+				for ( auto Prime_p = PrimeL.begin(); Prime_p != PrimeL.end(); ++ Prime_p) {
 					r.init (pri, *Prime_p);
 					do {
 						r.quoRem(quo,rem,Bonus,pri);
@@ -275,13 +270,11 @@ namespace LinBox
 			if (r_num.size()!=A. coldim()) return lif=0;
 			Integer r_den;
 			DVect b(r,A.rowdim());
-			typename DVect::iterator b_p;
-			//typename Vector::const_iterator Prime_p;
 
 			Integer pri, quo, rem;
 
 			// assign b to be a random vector
-			for (b_p = b.begin(); b_p != b.end(); ++ b_p) {
+			for (auto b_p = b.begin(); b_p != b.end(); ++ b_p) {
 //				* b_p = rand() % 268435456 - 134217728; // may need to change to use ring's random gen.
 //				// dpritcha, 2004-07-26
                 _gen( * b_p );

linbox/algorithms/lifting-container.h

@@ -155,14 +155,15 @@ namespace LinBox
             this->_intRing.convert(Prime,_p);
 
             auto hb = RationalSolveHadamardBound(A, b);
-            N = Integer(1) << static_cast<uint64_t>(std::ceil(hb.numLogBound));
-            D = Integer(1) << static_cast<uint64_t>(std::ceil(hb.denLogBound));
+            N = hb.numBound;
+            D = hb.denBound;
 
             // L = N * D * 2
             // _length = logp(L, Prime) = log2(L) * ln(2) / ln(Prime)
             double primeLog2 = Givaro::logtwo(Prime);
-            _length = std::ceil((1 + hb.numLogBound + hb.denLogBound) / primeLog2); // round up instead of down
+            _length = std::ceil(hb.solutionLogBound / primeLog2); // round up instead of down
 #ifdef DEBUG_LC
+			std::cout << "_length "<< _length << std::endl;
 			std::cout<<" norms computed, p = "<<_p<<"\n";
 			std::cout<<" N = "<<N<<", D = "<<D<<", length = "<<_length<<"\n";
 			_matA.write(std::cout<<"A:=", Tag::FileFormat::Maple) << std::endl;
@@ -226,7 +227,7 @@ namespace LinBox
 
 				// compute v2 = _matA * digit
 				IVector v2 (_lc.ring(),_lc._matA.rowdim());
-				_lc._MAD.applyV(v2,digit, _res);
+				_lc._MAD.applyV(v2,digit, _res); // @fixme This third parameter makes no sense!
 
 #ifdef DEBUG_LC
 
@@ -285,34 +286,6 @@ namespace LinBox
 
 		};
 
-		/*- @brief Bit manipulation function for possible use in optimization.
-		 * efficiently pulls out continuous blocks of bits, from lsb to msb inclusive
-		 * least significant bits start at index 0, so msb >= lsb
-		 * if any bits with index >= 8*numBytes are asked for they will be zeroes
-		 */
-#if 0
-		static long long bytesToBits(unsigned char * byteArray, size_t numBytes, size_t lsb, size_t msb) {
-			linbox_check(msb >= lsb);
-			size_t lsbi = lsb >> 3;
-			size_t msbi = msb >> 3;
-			if (msbi == lsbi)
-				if (msbi >= numBytes)
-					return 0;
-				else
-					return (byteArray[lsbi] >> (lsb & 7)) & ((1 << (msb - lsb + 1)) - 1);
-
-			long long result = (msbi < numBytes) ? (byteArray[msbi] & ((1 << ((msb & 7)+1)) - 1)) : 0;
-			for (size_t i=msbi-1; i>lsbi; i--) {
-				result <<= 8;
-				result |= (i < numBytes) ? byteArray[i] : 0;
-			}
-			result <<= 8 - (lsb & 7);
-			result |= (lsbi < numBytes) ? (byteArray[lsbi] >> (lsb & 7)) : 0;
-
-			return result;
-		}
-#endif
-
 		const_iterator begin() const
 		{
 			return const_iterator(*this);