diff --git a/configure.ac b/configure.ac index e97c7a699..8be1e864e 100644 --- a/configure.ac +++ b/configure.ac @@ -309,6 +309,7 @@ examples/data/Makefile doc/Makefile linbox/Makefile linbox/algorithms/Makefile +linbox/algorithms/dixon-solver/Makefile linbox/algorithms/gauss/Makefile linbox/algorithms/matrix-blas3/Makefile linbox/algorithms/opencl-kernels/Makefile diff --git a/examples/solver/t-rdisolve.C b/examples/solver/t-rdisolve.C index 34437f13d..3738322a5 100644 --- a/examples/solver/t-rdisolve.C +++ b/examples/solver/t-rdisolve.C @@ -184,7 +184,7 @@ int test() cout<<"' and Z/pZ of type '"; F.write(cout)<<"'"< QSolver; + typedef DixonSolver QSolver; typedef DiophantineSolver ZSolver; //typedef std::vector > FractionVector; diff --git a/linbox/algorithms/Makefile.am b/linbox/algorithms/Makefile.am index 542bcc874..2d0d3206f 100644 --- a/linbox/algorithms/Makefile.am +++ b/linbox/algorithms/Makefile.am @@ -21,7 +21,7 @@ pkgincludesubdir=$(pkgincludedir)/algorithms -SUBDIRS=gauss opencl-kernels matrix-blas3 polynomial-matrix +SUBDIRS=dixon-solver gauss opencl-kernels matrix-blas3 polynomial-matrix # IML noinst_LTLIBRARIES=libalgorithms.la @@ -127,7 +127,6 @@ pkgincludesub_HEADERS = \ rational-reconstruction2.h \ rational-reconstruction-base.h \ rational-reconstruction.h \ - rational-solver2.h \ rational-solver-adaptive.h \ rational-solver.h \ rational-solver.inl \ diff --git a/linbox/algorithms/det-rational.h b/linbox/algorithms/det-rational.h index 6b3f5a443..d6be12dc6 100644 --- a/linbox/algorithms/det-rational.h +++ b/linbox/algorithms/det-rational.h @@ -243,8 +243,8 @@ namespace LinBox Integer lif = 1; if ((s1 > 4*s2) && (!term)){ //cout << "lif " << std::flush; - RationalSolver < Givaro::IntegerDom , Givaro::Modular, PrimeIterator, Method::Dixon > RSolver; - LastInvariantFactor < Givaro::IntegerDom ,RationalSolver < Givaro::IntegerDom, Givaro::Modular, PrimeIterator, Method::Dixon > > LIF(RSolver); + DixonSolver < Givaro::IntegerDom , Givaro::Modular, PrimeIterator, Method::DenseElimination > RSolver; + LastInvariantFactor < Givaro::IntegerDom ,DixonSolver < Givaro::IntegerDom, Givaro::Modular, PrimeIterator, Method::DenseElimination > > LIF(RSolver); IVect r_num2 (Z,Atilde. coldim()); t1.clear(); t1.start(); diff --git a/linbox/algorithms/diophantine-solver.inl b/linbox/algorithms/diophantine-solver.inl index 96512dafd..03325deef 100644 --- a/linbox/algorithms/diophantine-solver.inl +++ b/linbox/algorithms/diophantine-solver.inl @@ -77,7 +77,13 @@ namespace LinBox SolverReturnStatus status; //this should eliminate all inconsistent systems; when level == SL_MONTECARLO maybe not. - status = _rationalSolver.monolithicSolve(x, den, A, b, (level >= SL_LASVEGAS), true, maxPrimes, level); + Method::Dixon method; + method.certifyMinimalDenominator = (level >= SL_LASVEGAS); + method.certifyInconsistency = (level >= SL_LASVEGAS); + method.singularSolutionType = SingularSolutionType::Random; + method.trialsBeforeFailure = maxPrimes; + + status = _rationalSolver.monolithicSolve(x, den, A, b, method); if (status != SS_OK) { if (status == SS_FAILED && maxPrimes > 2) std::cout << "ERROR, failed to find original solution and maxPrimes is not too small!" << std::endl; @@ -113,7 +119,14 @@ namespace LinBox int boredom = 0; //used in monte carlo, when we assume there's a diophantine solution while (! _ring.areEqual(upperDenBound, lowerDenBound)) { _rationalSolver.chooseNewPrime(); - status = _rationalSolver.monolithicSolve(x, den, A, b, (level >= SL_LASVEGAS), true, 1, level); + + Method::Dixon method; + method.certifyMinimalDenominator = (level >= SL_LASVEGAS); + method.certifyInconsistency = (level >= SL_LASVEGAS); + method.singularSolutionType = SingularSolutionType::Random; + method.trialsBeforeFailure = 1; + + status = _rationalSolver.monolithicSolve(x, den, A, b, method); numSolutionsNeeded++; #ifdef DEBUG_DIO std::cout << '.' ; diff --git a/linbox/algorithms/dixon-solver/Makefile.am b/linbox/algorithms/dixon-solver/Makefile.am new file mode 100644 index 000000000..49d67166a --- /dev/null +++ b/linbox/algorithms/dixon-solver/Makefile.am @@ -0,0 +1,26 @@ +# Copyright (c) 2010 the LinBox group +# written by Brice Boyer +# ========LICENCE======== +# This file is part of the library LinBox. +# +# LinBox is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA +# ========LICENCE======== + +pkgincludesubdir=$(pkgincludedir)/algorithms/dixon-solver + +pkgincludesub_HEADERS = \ + dixon-solver-dense.h \ + dixon-solver-dense.inl \ + dixon-solver-symbolic-numeric.h diff --git a/linbox/algorithms/dixon-solver/dixon-solver-dense.h b/linbox/algorithms/dixon-solver/dixon-solver-dense.h new file mode 100644 index 000000000..0d378e69f --- /dev/null +++ b/linbox/algorithms/dixon-solver/dixon-solver-dense.h @@ -0,0 +1,363 @@ +/* + * Copyright (C) LinBox Team + * + * ========LICENCE======== + * This file is part of the library LinBox. + * + * LinBox is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * ========LICENCE======== + */ + +#pragma once + +#include "../rational-solver.h" + +namespace LinBox { +#ifdef RSTIMING + class DixonTimer { + public: + mutable Timer ttSetup, ttRecon, ttGetDigit, ttGetDigitConvert, ttRingApply, ttRingOther; + mutable int rec_elt; + void clear() + { + ttSetup.clear(); + ttRecon.clear(); + ttGetDigit.clear(); + ttGetDigitConvert.clear(); + ttRingOther.clear(); + ttRingApply.clear(); + rec_elt = 0; + } + + template + void update(RR& rr, LC& lc) + { + ttSetup += lc.ttSetup; + ttRecon += rr.ttRecon; + rec_elt += rr._num_rec; + ttGetDigit += lc.ttGetDigit; + ttGetDigitConvert += lc.ttGetDigitConvert; + ttRingOther += lc.ttRingOther; + ttRingApply += lc.ttRingApply; + } + }; +#endif + + /** \brief partial specialization of p-adic based solver with Dixon algorithm. + * + * See the following reference for details on this algorithm: + * @bib + * - John D. Dixon Exact Solution of linear equations using p-adic + * expansions. Numerische Mathematik, volume 40, pages 137-141, + * 1982. + * + */ + template + class DixonSolver { + + public: + typedef Ring RingType; + typedef typename Ring::Element Integer; + typedef typename Field::Element Element; + typedef typename RandomPrime::Prime_Type Prime; + + // polymorphic 'certificate' generated when level >= SL_CERTIFIED + // certificate of inconsistency when any solving routine returns SS_INCONSISTENT + // certificate of minimal denominator when findRandomSolutionAndCertificate is called & + // return is SS_OK + mutable VectorFraction lastCertificate; + + // next 2 fields generated only by findRandomSolutionAndCertificate, when return is SS_OK + mutable Integer lastZBNumer; // filled in if level >= SL_CERTIFIED + mutable Integer lastCertifiedDenFactor; // filled in if level >= SL_LASVEGAS + // note: lastCertificate * b = lastZBNumer / lastCertifiedDenFactor, in lowest form + + protected: + mutable RandomPrime _genprime; + mutable Prime _prime; + Ring _ring; + Field _field; + + BlasMatrixDomain _bmdf; + +#ifdef RSTIMING + mutable Timer tSetup, ttSetup, tFastInvert, + ttFastInvert, // only done in deterministic or inconsistent + tCheckConsistency, ttCheckConsistency, // includes lifting the certificate + tMakeConditioner, ttMakeConditioner, tInvertBP, ttInvertBP, // only done in random + tCheckAnswer, ttCheckAnswer, tCertSetup, + ttCertSetup, // remaining 3 only done when makeMinDenomCert = true + tCertMaking, ttCertMaking, + + tNonsingularSetup, ttNonsingularSetup, tNonsingularInv, ttNonsingularInv, + + totalTimer; + + mutable DixonTimer ttConsistencySolve, ttSystemSolve, ttCertSolve, ttNonsingularSolve; +#endif + + public: + /** Constructor + * @param r a Ring, set by default + * @param rp a RandomPrime generator, set by default + */ + DixonSolver(const Ring& r = Ring(), const RandomPrime& rp = RandomPrime()) + : lastCertificate(r, 0) + , _genprime(rp) + , _ring(r) + { + _genprime.setBits(FieldTraits::bestBitSize()); + _prime = *_genprime; + ++_genprime; +#ifdef RSTIMING + clearTimers(); +#endif + } + + /** Constructor, trying the prime p first + * @param p a Prime + * @param r a Ring, set by default + * @param rp a RandomPrime generator, set by default + */ + DixonSolver(const Prime& p, const Ring& r = Ring(), const RandomPrime& rp = RandomPrime()) + : lastCertificate(r, 0) + , _genprime(rp) + , _prime(p) + , _ring(r) + { + _genprime.setBits(FieldTraits::bestBitSize()); +#ifdef RSTIMING + clearTimers(); +#endif + } + + /** Solve a linear system \c Ax=b over quotient field of a ring. + * + * @param num Vector of numerators of the solution + * @param den The common denominator. 1/den * num is the rational solution of \c Ax=b. + * @param A Matrix of linear system + * @param b Right-hand side of system + * @param s + * @param maxPrimes maximum number of moduli to try + * @param level level of certification to be used + * + * @return status of solution. if \c (return != SS_FAILED), and \c (level >= SL_LASVEGAS), + * solution is guaranteed correct. \c SS_FAILED - all primes used were bad \c SS_OK - + * solution found. \c SS_INCONSISTENT - system appreared inconsistent. certificate is in \p + * lastCertificate if \c (level >= SL_CERTIFIED) + */ + template + SolverReturnStatus solve(Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, const bool s = false, + const int maxPrimes = DEFAULT_MAXPRIMES, const SolverLevel level = SL_DEFAULT); + + // @fixme Can we remove that? + /** overload so that the bool 'oldMatrix' argument is not accidentally set to true */ + template + SolverReturnStatus solve(Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, const int maxPrimes, + const SolverLevel level = SL_DEFAULT) + { + return solve(num, den, A, b, false, maxPrimes, level); + } + + /** Solve a nonsingular, square linear system \c Ax=b over quotient field of a ring. + * + * @param num Vector of numerators of the solution + * @param den The common denominator. 1/den * num is the rational + * solution of Ax = b + * @param A Matrix of linear system (it must be square) + * @param b Right-hand side of system + * @param s unused + * @param maxPrimes maximum number of moduli to try + * + * @return status of solution : + * - \c SS_FAILED all primes used were bad; + * - \c SS_OK solution found, guaranteed correct; + * - \c SS_SINGULAR system appreared singular mod all primes. + * . + */ + template + SolverReturnStatus solveNonsingular(Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, bool s = false, + int maxPrimes = DEFAULT_MAXPRIMES); + + /** Solve a general rectangular linear system \c Ax=b over quotient field of a ring. + * If A is known to be square and nonsingular, calling solveNonsingular is more efficient. + * + * @param num Vector of numerators of the solution + * @param den The common denominator. 1/den * num is the rational + * solution of Ax = b + * @param A Matrix of linear system + * @param b Right-hand side of system + * @param maxPrimes maximum number of moduli to try + * @param level level of certification to be used + * + * @return status of solution. if (return != SS_FAILED), and (level >= + * SL_LASVEGAS), solution is guaranteed correct. + * - \c SS_FAILED all primes used were bad + * - \c SS_OK solution found. + * - \c SS_INCONSISTENT system appreared inconsistent. certificate is in \p + * lastCertificate if (level >= SL_CERTIFIED) + */ + template + SolverReturnStatus solveSingular(Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, + int maxPrimes = DEFAULT_MAXPRIMES, const SolverLevel level = SL_DEFAULT); + + /** Find a random solution of the general linear system \c Ax=b over quotient field of a + * ring. + * + * @param num Vector of numerators of the solution + * @param den The common denominator. 1/den * num is the rational solution of + * Ax = b. + * @param A Matrix of linear system + * @param b Right-hand side of system + * @param maxPrimes maximum number of moduli to try + * @param level level of certification to be used + * + * @return status of solution. if (return != SS_FAILED), and (level >= + * SL_LASVEGAS), solution is guaranteed correct. + * - \c SS_FAILED all primes used were bad + * - \c SS_OK solution found. + * - \c SS_INCONSISTENT system appreared inconsistent. certificate is in lastCertificate + * if (level >= SL_CERTIFIED) + */ + template + SolverReturnStatus findRandomSolution(Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, + int maxPrimes = DEFAULT_MAXPRIMES, const SolverLevel level = SL_DEFAULT); + + /** Big solving routine to perform random solving and certificate generation. + * Same arguments and return as findRandomSolution, except + * + * @param num Vector of numerators of the solution + * @param den The common denominator. 1/den * num is the rational solution of + * Ax = b + * @param A + * @param b + * @param randomSolution parameter to determine whether to randomize or not (since + * solveSingular calls this function as well) + * @param makeMinDenomCert determines whether a partial certificate for the minimal + * denominator of a rational solution is made + * @param maxPrimes + * @param level + * + * When (randomSolution == true && makeMinDenomCert == true), + * - If (level == SL_MONTECARLO) this function has the same effect as calling + * findRandomSolution. + * - If (level >= SL_LASVEGAS && return == SS_OK), \c lastCertifiedDenFactor + * contains a certified factor of the min-solution's denominator. + * - If (level >= SL_CERTIFIED && return == SS_OK), \c lastZBNumer and \c + * lastCertificate are updated as well. + * + */ + template + SolverReturnStatus monolithicSolve(Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, Method::Dixon method); + + Ring getRing() { return _ring; } + + void chooseNewPrime() + { + ++_genprime; + _prime = *_genprime; + } + +#ifdef RSTIMING + void clearTimers() + { + ttSetup.clear(); + ttFastInvert.clear(); + ttCheckConsistency.clear(); + ttMakeConditioner.clear(); + ttInvertBP.clear(); + ttCheckAnswer.clear(); + ttCertSetup.clear(); + ttCertMaking.clear(); + ttNonsingularSetup.clear(); + ttNonsingularInv.clear(); + + ttConsistencySolve.clear(); + ttSystemSolve.clear(); + ttCertSolve.clear(); + ttNonsingularSolve.clear(); + } + + public: + inline std::ostream& printTime(const Timer& timer, const char* title, std::ostream& os, const char* pref = "") + { + if (&timer != &totalTimer) totalTimer += timer; + if (timer.count() > 0) { + os << pref << title; + for (int i = strlen(title) + strlen(pref); i < 28; i++) os << ' '; + return os << timer << std::endl; + } + else + return os; + } + + inline std::ostream& printDixonTime(const DixonTimer& timer, const char* title, std::ostream& os) + { + if (timer.ttSetup.count() > 0) { + printTime(timer.ttSetup, "Setup", os, title); + printTime(timer.ttGetDigit, "Field Apply", os, title); + printTime(timer.ttGetDigitConvert, "Ring-Field-Ring Convert", os, title); + printTime(timer.ttRingApply, "Ring Apply", os, title); + printTime(timer.ttRingOther, "Ring Other", os, title); + printTime(timer.ttRecon, "Reconstruction", os, title); + os << " number of elt recontructed: " << timer.rec_elt << std::endl; + } + return os; + } + + std::ostream& reportTimes(std::ostream& os) + { + totalTimer.clear(); + printTime(ttNonsingularSetup, "NonsingularSetup", os); + printTime(ttNonsingularInv, "NonsingularInv", os); + printDixonTime(ttNonsingularSolve, "NS ", os); + printTime(ttSetup, "Setup", os); + printTime(ttFastInvert, "FastInvert", os); + printTime(ttCheckConsistency, "CheckConsistency", os); + printDixonTime(ttConsistencySolve, "INC ", os); + printTime(ttMakeConditioner, "MakeConditioner", os); + printTime(ttInvertBP, "InvertBP", os); + printDixonTime(ttSystemSolve, "SYS ", os); + printTime(ttCheckAnswer, "CheckAnswer", os); + printTime(ttCertSetup, "CertSetup", os); + printDixonTime(ttCertSolve, "CER ", os); + printTime(ttCertMaking, "CertMaking", os); + printTime(totalTimer, "TOTAL", os); + return os; + } +#endif + + private: + /// Internal usage + template + SolverReturnStatus solveApparentlyInconsistent(const BlasMatrix& A, TAS& tas, BlasMatrix* Atp_minor_inv, + size_t rank, const MethodBase& method); + + /// Internal usage + /// @note P seems to be a preconditioner: a random matrix filled with 0 or 1 + /// @note Please note that Atp_minor_inv content *might* be changed after this function + /// call. + template + void makeConditioner(BlasMatrix& A_minor, BlasMatrix*& Ap_minor_inv, BlasMatrix*& B, + BlasMatrix*& P, const BlasMatrix& A, TAS& tas, BlasMatrix* Atp_minor_inv, + size_t rank, const MethodBase& method); + + template + void certifyMinimalDenominator(const BlasMatrix& A, const Vector2& b, const TAS& tas, const BlasMatrix& B, + BlasMatrix& A_minor, BlasMatrix& Ap_minor_inv, size_t rank); + }; +} + +#include "./dixon-solver-dense.inl" diff --git a/linbox/algorithms/dixon-solver/dixon-solver-dense.inl b/linbox/algorithms/dixon-solver/dixon-solver-dense.inl new file mode 100644 index 000000000..f1ab9266e --- /dev/null +++ b/linbox/algorithms/dixon-solver/dixon-solver-dense.inl @@ -0,0 +1,803 @@ +/* + * Copyright (C) LinBox Team + * + * ========LICENCE======== + * This file is part of the library LinBox. + * + * LinBox is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * ========LICENCE======== + */ + +#include "linbox/linbox-config.h" +#include "linbox/util/debug.h" + +#include "linbox/algorithms/lifting-container.h" +#include "linbox/algorithms/matrix-inverse.h" +#include "linbox/algorithms/rational-reconstruction.h" + +namespace LinBox { + + template + template + SolverReturnStatus DixonSolver::solve(Vector1& num, Integer& den, + const IMatrix& A, const Vector2& b, + const bool old, int maxP, + const SolverLevel level) + { + SolverReturnStatus status; + int maxPrimes = maxP; + while (maxPrimes > 0) { + auto nonSingularResult = (A.rowdim() == A.coldim()) ? solveNonsingular(num, den, A, b, old, maxPrimes) : SS_SINGULAR; + switch (nonSingularResult) { + case SS_OK: return SS_OK; break; + + case SS_SINGULAR: + status = solveSingular(num, den, A, b, maxPrimes, level); + if (status != SS_FAILED) return status; + break; + + case SS_FAILED: break; + + default: throw LinboxError("Bad return value from solveNonsingular"); + } + maxPrimes--; + if (maxPrimes > 0) chooseNewPrime(); + } + return SS_FAILED; + } + + template + template + SolverReturnStatus DixonSolver::solveNonsingular( + Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, bool oldMatrix, int maxPrimes) + { + + int trials = 0, notfr; + + // history sensitive data for optimal reason + // static const IMatrix* IMP; + + BlasMatrix* FMP = NULL; + Field* F = NULL; + + do { +#ifdef RSTIMING + tNonsingularSetup.start(); +#endif + // typedef typename Field::Element Element; + // typedef typename Ring::Element Integer; + + // checking size of system + linbox_check(A.rowdim() == A.coldim()); + linbox_check(A.rowdim() == b.size()); + + LinBox::integer tmp; + + // if input matrix A is different one. + if (!oldMatrix) { + if (trials == maxPrimes) return SS_SINGULAR; + if (trials != 0) chooseNewPrime(); + ++trials; + + // Could delete a non allocated matrix -> segfault + if (FMP != NULL) delete FMP; + + // IMP = &A; + + if (F != NULL) delete F; + + F = new Field(_prime); + FMP = new BlasMatrix(*F, A.rowdim(), A.coldim()); + MatrixHom::map(*FMP, A); // use MatrixHom to reduce matrix PG 2005-06-16 + + BlasMatrix* invA = new BlasMatrix(*F, A.rowdim(), A.coldim()); + BlasMatrixDomain BMDF(*F); +#ifdef RSTIMING + tNonsingularSetup.stop(); + ttNonsingularSetup += tNonsingularSetup; + tNonsingularInv.start(); +#endif + assert(FMP != NULL); + BMDF.invin(*invA, *FMP, notfr); // notfr <- nullity + delete FMP; + FMP = invA; + +#ifdef RSTIMING + tNonsingularInv.stop(); + ttNonsingularInv += tNonsingularInv; +#endif + } + else { +#ifdef RSTIMING + tNonsingularSetup.stop(); + ttNonsingularSetup += tNonsingularSetup; +#endif + notfr = 0; + } + } while (notfr); + + typedef DixonLiftingContainer> LiftingContainer; + LiftingContainer lc(_ring, *F, A, *FMP, b, _prime); + RationalReconstruction re(lc); + if (!re.getRational(num, den, 0)) { + delete FMP; + return SS_FAILED; + } +#ifdef RSTIMING + ttNonsingularSolve.update(re, lc); +#endif + if (F != NULL) delete F; + if (FMP != NULL) delete FMP; + return SS_OK; + } + + template + template + SolverReturnStatus DixonSolver::solveSingular( + Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, int maxPrimes, const SolverLevel level) + { + Method::Dixon m; + m.certifyMinimalDenominator = false; + m.certifyInconsistency = false; + m.trialsBeforeFailure = maxPrimes; + return monolithicSolve(num, den, A, b, m); + } + + template + template + SolverReturnStatus DixonSolver::findRandomSolution( + Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, int maxPrimes, const SolverLevel level) + { + Method::Dixon m; + m.singularSolutionType = SingularSolutionType::Random; + m.certifyMinimalDenominator = false; + m.certifyInconsistency = true; + m.trialsBeforeFailure = maxPrimes; + return monolithicSolve(num, den, A, b, m); + } + + // TAS stands for Transpose Augmented System (A|b)t + // this provides a factorization (A|b) = Pt . Lt . Ut . Qt + // such that + // - Q . (A|b) . P has nonzero principal minors up to TAS.rank() + // - P permutes b to the (TAS.rank())th column of A iff the system is inconsistent mod p + + template + class TransposeAugmentedSystem { + public: + BlasMatrix* factors = nullptr; + PLUQMatrix* PLUQ = nullptr; + + BlasPermutation P; + BlasPermutation Q; + + std::vector srcRow; + std::vector srcCol; + + public: + template + TransposeAugmentedSystem(Ring& R, Field& _field, const IMatrix& A, const IVector& b) + : P(A.coldim() + 1) + , Q(A.rowdim()) + { + factors = new BlasMatrix(_field, A.coldim() + 1, A.rowdim()); + + BlasMatrix Ap(A, _field); // Getting into the field + + // Setting factors = [Ap|0]t + for (size_t i = 0; i < A.rowdim(); ++i) + for (size_t j = 0; j < A.coldim(); ++j) factors->setEntry(j, i, Ap.getEntry(i, j)); + + // And then factors = [Ap|bp]t + Integer tmpInteger; + for (size_t i = 0; i < A.rowdim(); ++i) { + typename Field::Element tmpElement; + _field.init(tmpElement, R.convert(tmpInteger, b[i])); + factors->setEntry(A.coldim(), i, tmpElement); + } + + // Getting factors -> P L U Q + PLUQ = new PLUQMatrix(*factors, P, Q); + + // Used to check consistency + srcRow.resize(A.rowdim()); + srcCol.resize(A.coldim() + 1); + FFPACK::LAPACKPerm2MathPerm(srcRow.data(), Q.getStorage().data(), srcRow.size()); + FFPACK::LAPACKPerm2MathPerm(srcCol.data(), P.getStorage().data(), srcCol.size()); + + // @note srcCol/srcRow now hold permutations (indices of (A|b)t) + } + + ~TransposeAugmentedSystem() + { + delete factors; + delete PLUQ; + } + + size_t rank() { return PLUQ->getRank(); } + }; + + // Returns: + // - FAILED if A != 0 + // - INCONSISTENT if A == 0 and b != 0 (also sets certificate) + // - OK if A == 0 and b == 0 + template + SolverReturnStatus certifyEmpty(const Matrix& A, const Vector& b, const MethodBase& method, Integer& certifiedDenFactor) + { + const Ring& R = A.field(); + + // In monte carlo, we assume A is actually empty. + bool aEmpty = true; + if (method.certifyInconsistency) { + MatrixDomain MD(R); + aEmpty = MD.isZero(A); + } + + if (!aEmpty) { + return SS_FAILED; + } + + // @fixme Use VectorDomain::isZero + for (size_t i = 0; i < b.size(); ++i) { + if (!R.areEqual(b[i], R.zero)) { + return SS_INCONSISTENT; + } + } + + // System is consistent + if (method.certifyInconsistency) { + R.assign(certifiedDenFactor, R.one); + } + + return SS_OK; + } + + // --------------------------------------------- + // SS_FAILED means that we will need to try a new prime + template + template + SolverReturnStatus DixonSolver::solveApparentlyInconsistent( + const BlasMatrix& A, TAS& tas, BlasMatrix* Atp_minor_inv, size_t rank, const MethodBase& method) + { + using LiftingContainer = DixonLiftingContainer, BlasMatrix>; + + if (!method.certifyInconsistency) return SS_INCONSISTENT; + + // @fixme Put these as class members! + BlasMatrixDomain BMDI(_ring); + BlasApply BAR(_ring); + +#ifdef RSTIMING + tCheckConsistency.start(); +#endif + + BlasVector zt(_ring, rank); + for (size_t i = 0; i < rank; ++i) _ring.assign(zt[i], A.getEntry(tas.srcRow[rank], tas.srcCol[i])); + + BlasMatrix At_minor(_ring, rank, rank); + for (size_t i = 0; i < rank; ++i) + for (size_t j = 0; j < rank; ++j) _ring.assign(At_minor.refEntry(j, i), A.getEntry(tas.srcRow[i], tas.srcCol[j])); + +#ifdef RSTIMING + tCheckConsistency.stop(); + ttCheckConsistency += tCheckConsistency; +#endif + + LiftingContainer lc(_ring, _field, At_minor, *Atp_minor_inv, zt, _prime); + RationalReconstruction re(lc); + + BlasVector shortNum(A.field(), rank); + Integer shortDen; + + // Dirty, but should not be called under normal circumstances + if (!re.getRational(shortNum, shortDen, 0)) { + return SS_FAILED; + } + +#ifdef RSTIMING + ttConsistencySolve.update(re, lc); + tCheckConsistency.start(); +#endif + + // Build up certificate + VectorFraction cert(_ring, shortNum.size()); + cert.numer = shortNum; + cert.denom = shortDen; + cert.numer.resize(A.rowdim()); + _ring.subin(cert.numer[rank], cert.denom); + _ring.assign(cert.denom, _ring.one); + + BMDI.mulin_left(cert.numer, tas.Q); + +#ifdef DEBUG_INC + cert.write(std::cout << "cert:") << std::endl; +#endif + + bool certifies = true; // check certificate + BlasVector certnumer_A(_ring, A.coldim()); + BAR.applyVTrans(certnumer_A, A, cert.numer); + typename BlasVector::iterator cai = certnumer_A.begin(); + for (size_t i = 0; certifies && i < A.coldim(); ++i, ++cai) { + certifies = certifies && _ring.isZero(*cai); + } + +#ifdef RSTIMING + tCheckConsistency.stop(); + ttCheckConsistency += tCheckConsistency; +#endif + + if (certifies) { + if (method.certifyInconsistency) lastCertificate.copy(cert); + return SS_INCONSISTENT; + } + commentator().report(Commentator::LEVEL_IMPORTANT, PARTIAL_RESULT) + << "system is suspected to be inconsistent but it was only a bad prime" << std::endl; + + // @fixme Try new prime, is SS_FAILED the right API? + // Analogous to u.A12 != A22 in Muld.+Storj. + return SS_FAILED; + } + + template + template + void DixonSolver::makeConditioner( + BlasMatrix& A_minor, BlasMatrix*& Ap_minor_inv, BlasMatrix*& B, BlasMatrix*& P, + const BlasMatrix& A, TAS& tas, BlasMatrix* Atp_minor_inv, size_t rank, const MethodBase& method) + { +#ifdef RSTIMING + tMakeConditioner.start(); +#endif + + if (method.singularSolutionType != SingularSolutionType::Random) { + // Transpose Atp_minor_inv to get Ap_minor_inv + // @note minor inv = L1\U1 + Element _rtmp; + Ap_minor_inv = Atp_minor_inv; + for (size_t i = 0; i < rank; ++i) { + for (size_t j = 0; j < i; ++j) { + Ap_minor_inv->getEntry(_rtmp, i, j); + Ap_minor_inv->setEntry(i, j, Ap_minor_inv->refEntry(j, i)); + Ap_minor_inv->setEntry(j, i, _rtmp); + } + } + + // permute original entries into A_minor + // @note A_minor = Pt A Qt + for (size_t i = 0; i < rank; ++i) + for (size_t j = 0; j < rank; ++j) _ring.assign(A_minor.refEntry(i, j), A.getEntry(tas.srcRow[i], tas.srcCol[j])); +#ifdef RSTIMING + tMakeConditioner.stop(); + ttMakeConditioner += tMakeConditioner; +#endif + + if (method.certifyMinimalDenominator) { + B = new BlasMatrix(_ring, rank, A.coldim()); + for (size_t i = 0; i < rank; ++i) + for (size_t j = 0; j < A.coldim(); ++j) _ring.assign(B->refEntry(i, j), A.getEntry(tas.srcRow[i], j)); + } + // @note B = Pt A + } + else { // if randomSolution == true + P = new BlasMatrix(_ring, A.coldim(), rank); + B = new BlasMatrix(_ring, rank, A.coldim()); + BlasMatrix Ap_minor(_field, rank, rank); + Ap_minor_inv = new BlasMatrix(_field, rank, rank); + + LinBox::integer tmp2 = 0; + size_t maxBitSize = 0; + for (size_t i = 0; i < rank; ++i) + for (size_t j = 0; j < A.coldim(); ++j) { + _ring.assign(B->refEntry(i, j), A.getEntry(tas.srcRow[i], j)); + _ring.convert(tmp2, A.getEntry(tas.srcRow[i], j)); + maxBitSize = std::max(maxBitSize, tmp2.bitsize()); + } + // @note B = Pt A +#ifdef RSTIMING + bool firstLoop = true; +#endif + // prepare B to be preconditionned through BLAS matrix mul + MatrixApplyDomain> MAD(_ring, *B); + MAD.setup(2); // @fixme Useless? + + int nullity; + do { // O(1) loops of this preconditioner expected +#ifdef RSTIMING + if (firstLoop) + firstLoop = false; + else + tMakeConditioner.start(); +#endif + // compute P a n*r random matrix of entry in [0,1] + typename BlasMatrix::Iterator iter; + for (iter = P->Begin(); iter != P->End(); ++iter) { + if (rand() % 2 == 1) + _ring.assign(*iter, _ring.one); + else + _ring.assign(*iter, _ring.zero); + } + + // compute A_minor = B.P + MAD.applyM(A_minor, *P); + + // set Ap_minor = A_minor mod p, try to compute inverse + for (size_t i = 0; i < rank; ++i) + for (size_t j = 0; j < rank; ++j) + _field.init(Ap_minor.refEntry(i, j), _ring.convert(tmp2, A_minor.getEntry(i, j))); +#ifdef RSTIMING + tMakeConditioner.stop(); + ttMakeConditioner += tMakeConditioner; + tInvertBP.start(); +#endif + + // @fixme Seems sad to be forced to specify these BlasMatrix& casts + _bmdf.inv((BlasMatrix&)*Ap_minor_inv, (BlasMatrix&)Ap_minor, nullity); + +#ifdef RSTIMING + tInvertBP.stop(); + ttInvertBP += tInvertBP; +#endif + } while (nullity > 0); + } + } + + template + template + void DixonSolver::certifyMinimalDenominator( + const BlasMatrix& A, const Vector2& b, const TAS& tas, const BlasMatrix& B, BlasMatrix& A_minor, + BlasMatrix& Ap_minor_inv, size_t rank) + { + // To make this certificate we solve with the same matrix as to get the + // solution, except transposed. +#ifdef RSTIMING + tCertSetup.start(); +#endif + + // @note We transpose Ap and A minors in-place because it won't be used anymore + Integer _rtmp; + Element _ftmp; + for (size_t i = 0; i < rank; ++i) + for (size_t j = 0; j < i; ++j) { + Ap_minor_inv.getEntry(_ftmp, i, j); + Ap_minor_inv.setEntry(i, j, Ap_minor_inv.refEntry(j, i)); + Ap_minor_inv.setEntry(j, i, _ftmp); + } + + for (size_t i = 0; i < rank; ++i) + for (size_t j = 0; j < i; ++j) { + A_minor.getEntry(_rtmp, i, j); + A_minor.setEntry(i, j, A_minor.refEntry(j, i)); + A_minor.setEntry(j, i, _rtmp); + } + + // we then try to create a partial certificate + // the correspondance with Algorithm MinimalSolution from Mulders/Storjohann: + // paper | here + // P | TAS_Q + // Q | transpose of TAS_P + // B | *B (== TAS_Q . A, but only top #rank rows) + // c | newb (== TAS_Q . b, but only top #rank rows) + // P | P + // q | q + // U | {0, 1} + // u | u + // z-hat | lastCertificate + + // we multiply the certificate by TAS_Qt at the end + // so it corresponds to b instead of newb + + // q in {0, 1}^rank + Givaro::ZRing Z; + BlasVector> q(Z, rank); + typename BlasVector>::iterator q_iter; + + bool allzero; + do { + allzero = true; + for (q_iter = q.begin(); q_iter != q.end(); ++q_iter) { + if (rand() > RAND_MAX / 2) { + _ring.assign((*q_iter), _ring.one); + allzero = false; + } + else + (*q_iter) = _ring.zero; + } + } while (allzero); + +#ifdef RSTIMING + tCertSetup.stop(); + ttCertSetup += tCertSetup; +#endif + + using LiftingContainer = DixonLiftingContainer, BlasMatrix>; + LiftingContainer lc2(_ring, _field, A_minor, Ap_minor_inv, q, _prime); + + RationalReconstruction rere(lc2); + Vector1 u_num(_ring, rank); + Integer u_den; + + // Failure + if (!rere.getRational(u_num, u_den, 0)) return; + +#ifdef RSTIMING + ttCertSolve.update(rere, lc2); + tCertMaking.start(); +#endif + + // remainder of code does z <- denom(partial_cert . Mr) * partial_cert * Qt + BlasApply BAR(_ring); + VectorFraction u_to_vf(_ring, u_num.size()); + u_to_vf.numer = u_num; + u_to_vf.denom = u_den; + BlasVector uB(_ring, A.coldim()); + BAR.applyVTrans(uB, B, u_to_vf.numer); + + Integer numergcd = _ring.zero; + vectorGcdIn(numergcd, _ring, uB); + + // denom(partial_cert . Mr) = partial_cert_to_vf.denom / numergcd + VectorFraction z(_ring, b.size()); // new constructor + u_to_vf.numer.resize(A.rowdim()); + + BlasMatrixDomain BMDI(_ring); + BMDI.mul(z.numer, u_to_vf.numer, tas.Q); + + z.denom = numergcd; + + lastCertificate.copy(z); + + // output new certified denom factor + Integer znumer_b, zbgcd; + VectorDomain VDR(_ring); + VDR.dotprod(znumer_b, z.numer, b); + _ring.gcd(zbgcd, znumer_b, z.denom); + _ring.div(lastCertifiedDenFactor, z.denom, zbgcd); + + _ring.div(lastZBNumer, znumer_b, zbgcd); +#ifdef RSTIMING + tCertMaking.stop(); + ttCertMaking += tCertMaking; +#endif + } + + // Most solving is done by the routine below. + // There used to be one for random and one for deterministic, but they have been merged to ease + // with + // repeated code (certifying inconsistency, optimization are 2 examples) + + template + template + SolverReturnStatus DixonSolver::monolithicSolve( + Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, Method::Dixon method) + { + using LiftingContainer = DixonLiftingContainer, BlasMatrix>; + + if (method.certifyMinimalDenominator && !method.certifyInconsistency) { + method.certifyInconsistency = true; + std::cerr << "WARNING: forcing certifyInconsistency due to certifyMinimalDenominator" << std::endl; + } + + size_t trials = 0; + while (trials < method.trialsBeforeFailure) { + if (trials != 0) chooseNewPrime(); + ++trials; + +#ifdef RSTIMING + tSetup.start(); +#endif + // ----- Build Transposed Augmented System (TAS) + + // checking size of system + linbox_check(A.rowdim() == b.size()); + + _field = Field(_prime); + _bmdf = BlasMatrixDomain(_field); + + LinBox::integer tmp; + BlasMatrixDomain BMDI(_ring); + BlasApply BAR(_ring); + MatrixDomain MD(_ring); + + BlasMatrix A_check(A); // used to check answer later + + // TAS stands for Transpose Augmented System (A|b)t + TransposeAugmentedSystem tas(_ring, _field, A, b); + +#ifdef RSTIMING + tSetup.stop(); + ttSetup += tSetup; +#endif + + // @note If permutation shows that b was needed, means b is not in the columns' span of + // A (=> Ax=b inconsistent) + bool appearsInconsistent = (tas.srcCol[tas.rank() - 1] == A.coldim()); + size_t rank = tas.rank() - (appearsInconsistent ? 1 : 0); + + // ----- Handle A == 0 mod p + + // Special case when A = 0, mod p. Deal with it to avoid later deadlock. + if (rank == 0) { + SolverReturnStatus status = certifyEmpty(A, b, method, lastCertifiedDenFactor); + + if (status == SS_FAILED) { + // A was empty mod p but not over Z, we try new prime. + continue; + } + else if (status == SS_INCONSISTENT) { + return SS_INCONSISTENT; + } + + // It is consistent, set 0 as a valid answer + _ring.assign(den, _ring.one); + for (typename Vector1::iterator p = num.begin(); p != num.end(); ++p) { + _ring.assign(*p, _ring.zero); + } + + return SS_OK; + } + + // ----- Certifying inconsistency + + std::unique_ptr> Atp_minor_inv = nullptr; + if ((appearsInconsistent && method.certifyInconsistency) + || method.singularSolutionType != SingularSolutionType::Random) { + // take advantage of the (PLUQ)t factorization to compute + // an inverse to the leading minor of (TAS_P . (A|b) . TAS_Q) +#ifdef RSTIMING + tFastInvert.start(); +#endif + + // @note std::make_unique is only C++14 + Atp_minor_inv = std::unique_ptr>(new BlasMatrix(_field, rank, rank)); + + FFLAS::fassign(_field, rank, rank, tas.factors->getPointer(), tas.factors->getStride(), Atp_minor_inv->getPointer(), Atp_minor_inv->getStride()); + FFPACK::ftrtri (_field, FFLAS::FflasUpper, FFLAS::FflasNonUnit, rank, Atp_minor_inv->getPointer(), Atp_minor_inv->getStride()); + FFPACK::ftrtri (_field, FFLAS::FflasLower, FFLAS::FflasUnit, rank, Atp_minor_inv->getPointer(), Atp_minor_inv->getStride()); + FFPACK::ftrtrm (_field, FFLAS::FflasLeft, FFLAS::FflasNonUnit, rank, Atp_minor_inv->getPointer(), Atp_minor_inv->getStride()); + +#ifdef RSTIMING + tFastInvert.stop(); + ttFastInvert += tFastInvert; +#endif + } + + // ----- Confirm inconsistency if it looks like it + + // If the system appears inconsistent, we either try a new prime, + // a validate the inconsistency of (A,b). + if (appearsInconsistent) { + auto status = solveApparentlyInconsistent(A_check, tas, Atp_minor_inv.get(), rank, method); + + // Failing means we should try a new prime + if (status == SS_FAILED) continue; + return status; + } + + // + // Starting from here, we know that the system is consistent mod p. + // + + // Build A_minor and Ap_minor_inv + BlasMatrix* B = nullptr; // @fixme Make that a std::unique_ptr + BlasMatrix* P = nullptr; // @fixme Make that a std::unique_ptr + BlasMatrix A_minor(_ring, rank, rank); // -- will have the full rank minor of A + BlasMatrix* Ap_minor_inv = nullptr; // -- will have inverse mod p of A_minor + makeConditioner(A_minor, Ap_minor_inv, B, P, A_check, tas, Atp_minor_inv.get(), rank, method); + + // Compute newb = (TAS_P.b)[0..(rank-1)] + BlasVector newb(b); + BMDI.mulin_right(tas.Q, newb); + newb.resize(rank); + + // ----- Do lifting on sub matrix + + BlasMatrix BBA_minor(A_minor); + LiftingContainer lc(_ring, _field, BBA_minor, *Ap_minor_inv, newb, _prime); + + // ----- Reconstruct rational + + RationalReconstruction re(lc); + VectorFraction resultVF(_ring, rank); + if (!re.getRational(resultVF.numer, resultVF.denom, 0)) { + // dirty, but should not be called + return SS_FAILED; + } + +#ifdef RSTIMING + ttSystemSolve.update(re, lc); + tCheckAnswer.start(); +#endif + + // ----- Build effective solution from sub matrix + + if (method.singularSolutionType != SingularSolutionType::Random) { + // short_answer = TAS_Q * short_answer + resultVF.numer.resize(A.coldim() + 1, _ring.zero); + BMDI.mulin_left(resultVF.numer, tas.P); + resultVF.numer.resize(A.coldim()); + } + else { + // short_answer = P * short_answer + BlasVector newNumer(_ring, A.coldim()); + BAR.applyV(newNumer, *P, resultVF.numer); + resultVF.numer = newNumer; + } + + // ----- Check consistency + + if (method.checkResult) { // check consistency + BlasVector A_times_xnumer(_ring, b.size()); + BAR.applyV(A_times_xnumer, A_check, resultVF.numer); + + Integer tmpi; + + typename Vector2::const_iterator ib = b.begin(); + typename BlasVector::iterator iAx = A_times_xnumer.begin(); + int thisrow = 0; + bool needNewPrime = false; + + for (; !needNewPrime && ib != b.end(); ++iAx, ++ib, ++thisrow) + if (!_ring.areEqual(_ring.mul(tmpi, *ib, resultVF.denom), *iAx)) { + // should attempt to certify inconsistency now + // as in "if [A31 | A32]y != b3" of step (4) + needNewPrime = true; + } + + if (needNewPrime) { + if (Ap_minor_inv != Atp_minor_inv.get()) { + delete Ap_minor_inv; + } + if (method.singularSolutionType == SingularSolutionType::Random) { + delete P; + } +#ifdef RSTIMING + tCheckAnswer.stop(); + ttCheckAnswer += tCheckAnswer; +#endif + continue; // go to start of main loop + } + } + +#ifdef RSTIMING + tCheckAnswer.stop(); + ttCheckAnswer += tCheckAnswer; +#endif + + // ----- We have the result values! + + num = resultVF.numer; + den = resultVF.denom; + + // ----- Checking minimal denominator + + if (method.certifyMinimalDenominator) { + certifyMinimalDenominator(A, b, tas, *B, A_minor, *Ap_minor_inv, rank); + } + + if (Ap_minor_inv != Atp_minor_inv.get()) { + delete Ap_minor_inv; + } + + if (method.singularSolutionType == SingularSolutionType::Random) { + delete P; + } + + // done making certificate, lets blow this popstand + return SS_OK; + } + + // All primes were bad + return SS_FAILED; + } +} diff --git a/linbox/algorithms/rational-solver2.h b/linbox/algorithms/dixon-solver/dixon-solver-symbolic-numeric.h similarity index 87% rename from linbox/algorithms/rational-solver2.h rename to linbox/algorithms/dixon-solver/dixon-solver-symbolic-numeric.h index 75d69052f..0627cc2d6 100644 --- a/linbox/algorithms/rational-solver2.h +++ b/linbox/algorithms/dixon-solver/dixon-solver-symbolic-numeric.h @@ -1,5 +1,4 @@ -/* linbox/algorithms/rational-solver2.h - * Copyright (C) 2010 LinBox +/* Copyright (C) 2010 LinBox * Author Z. Wan * * @@ -22,15 +21,15 @@ * ========LICENCE======== */ -/*! @file algorithms/rational-solver2.h +/*! @file algorithms/dixon-solver/dixon-solver-symbolic-numeric.h * @brief NO DOC * @bib * Implementation of the algorithm in manuscript, available at * http://www.cis.udel.edu/~wan/jsc_wan.ps */ -#ifndef __LINBOX_rational_solver2__H -#define __LINBOX_rational_solver2__H +#ifndef __LINBOX_dixon_solver_symbolic_numeric__H +#define __LINBOX_dixon_solver_symbolic_numeric__H #include #include @@ -56,7 +55,7 @@ namespace LinBox //template argument Field and RandomPrime are not used. //Keep it just for interface consistency. template - class RationalSolver { + class DixonSolver { protected: Ring r; @@ -64,7 +63,7 @@ namespace LinBox public: typedef typename Ring::Element Integer; - RationalSolver(const Ring& _r = Ring()) : + DixonSolver(const Ring& _r = Ring()) : r(_r) {} @@ -187,7 +186,7 @@ namespace LinBox }; #if __LINBOX_HAVE_CLAPACK template - inline int RationalSolver::cblas_dgeinv(double* M, int n) + inline int DixonSolver::cblas_dgeinv(double* M, int n) { enum CBLAS_ORDER order = CblasRowMajor; int lda = n; @@ -195,7 +194,7 @@ namespace LinBox int ierr = clapack_dgetrf (order, n, n, M, lda, P); if (ierr != 0) { commentator().report (Commentator::LEVEL_IMPORTANT, PARTIAL_RESULT) - /*std::cerr*/ << "In RationalSolver::cblas_dgeinv Matrix is not full rank" << std::endl; + /*std::cerr*/ << "In DixonSolver::cblas_dgeinv Matrix is not full rank" << std::endl; delete[] P ; return -1; } @@ -205,7 +204,7 @@ namespace LinBox } template - inline int RationalSolver::cblas_rsol (int n, const double* M, integer* numx, integer& denx, double* b) + inline int DixonSolver::cblas_rsol (int n, const double* M, integer* numx, integer& denx, double* b) { if (n < 1) return 0; double* IM = new double[n * n]; @@ -422,14 +421,14 @@ namespace LinBox template /* apply y <- Ax */ - inline int RationalSolver::cblas_dapply (int m, int n, const double* A, const double* x, double* y) + inline int DixonSolver::cblas_dapply (int m, int n, const double* A, const double* x, double* y) { cblas_dgemv (CblasRowMajor, CblasNoTrans, m, n, 1, A, n, x, 1, 0, y, 1); return 0; } template - inline int RationalSolver::cblas_mpzapply (int m, int n, const double* A, const integer* x, integer* y) + inline int DixonSolver::cblas_mpzapply (int m, int n, const double* A, const integer* x, integer* y) { const double* p_A; const integer* p_x; @@ -449,7 +448,7 @@ namespace LinBox template template - inline int RationalSolver::printvec (const Elt* v, int n) + inline int DixonSolver::printvec (const Elt* v, int n) { const Elt* p; std::cout << '['; @@ -461,7 +460,7 @@ namespace LinBox template //update num, *num <- *num * 2^shift + d - inline int RationalSolver::update_num (integer* num, int n, const double* d, int shift) + inline int DixonSolver::update_num (integer* num, int n, const double* d, int shift) { integer* p_mpz; integer tmp_mpz; @@ -476,7 +475,7 @@ namespace LinBox template //update r = r * shift - M d - inline int RationalSolver::update_r_int (double* r, int n, const double* M, const double* d, int shift) + inline int DixonSolver::update_r_int (double* r, int n, const double* M, const double* d, int shift) { double* p1; const double* p2; @@ -495,7 +494,7 @@ namespace LinBox template //update r = r * shift - M d - inline int RationalSolver::update_r_ll (double* r, int n, const double* M, const double* d, int shift) + inline int DixonSolver::update_r_ll (double* r, int n, const double* M, const double* d, int shift) { double* p1; const double* p2; @@ -513,7 +512,7 @@ namespace LinBox } template - inline double RationalSolver::cblas_dOOnorm(const double* M, int m, int n) + inline double DixonSolver::cblas_dOOnorm(const double* M, int m, int n) { double norm = 0; const double* p; @@ -528,13 +527,13 @@ namespace LinBox } template - inline double RationalSolver::cblas_dmax (const int N, const double* a, const int inc) + inline double DixonSolver::cblas_dmax (const int N, const double* a, const int inc) { return fabs(a[cblas_idamax (N, a, inc)]); } template - inline int RationalSolver::cblas_hbound (integer& b, int m, int n, const double* M) + inline int DixonSolver::cblas_hbound (integer& b, int m, int n, const double* M) { const double* p; integer tmp; @@ -551,7 +550,7 @@ namespace LinBox } }//LinBox -#endif //__LINBOX_rational_solver2__H +#endif //__LINBOX_dixon_solver_symbolic_numeric__H // Local Variables: // mode: C++ diff --git a/linbox/algorithms/hybrid-det.h b/linbox/algorithms/hybrid-det.h index 37a9a3406..d0adb83bd 100644 --- a/linbox/algorithms/hybrid-det.h +++ b/linbox/algorithms/hybrid-det.h @@ -254,13 +254,13 @@ namespace LinBox PrimeIterator genprime1(FieldTraits::bestBitSize(A.coldim())); Integers ZZ; - RationalSolver < Integers , mymodular, PrimeIterator, Method::Dixon > RSolver(A. field(), genprime); + DixonSolver < Integers , mymodular, PrimeIterator, Method::DenseElimination > RSolver(A. field(), genprime); #endif - RationalSolver < Integers , mymodular, PrimeIterator, Method::Dixon > RSolver; + DixonSolver < Integers , mymodular, PrimeIterator, Method::DenseElimination > RSolver; BlasVector r_num1 (A.field(),A. coldim()); - LastInvariantFactor < Integers ,RationalSolver < Integers, mymodular, PrimeIterator, Method::Dixon > > LIF(RSolver); + LastInvariantFactor < Integers ,DixonSolver < Integers, mymodular, PrimeIterator, Method::DenseElimination > > LIF(RSolver); #ifdef _LB_H_DET_TIMING BT.start(); #endif @@ -537,7 +537,7 @@ namespace LinBox PrimeIterator genprime1(FieldTraits::bestBitSize(A.coldim())); Integers ZZ; #endif - typedef RationalSolver < Integers , mymodular, PrimeIterator, Method::BlockHankel > Solver; + typedef DixonSolver < Integers , mymodular, PrimeIterator, Method::BlockHankel > Solver; Solver RSolver(A. field(), genprime); typename Vector:: Dense r_num1 (A. coldim()); diff --git a/linbox/algorithms/lifting-container.h b/linbox/algorithms/lifting-container.h index 6a9180092..d0527076f 100644 --- a/linbox/algorithms/lifting-container.h +++ b/linbox/algorithms/lifting-container.h @@ -37,6 +37,8 @@ #include "linbox/util/debug.h" #include "linbox/blackbox/apply.h" +#include "linbox/blackbox/diagonal.h" +#include "linbox/algorithms/matrix-hom.h" #include "linbox/algorithms/blackbox-container.h" #include "linbox/algorithms/massey-domain.h" #include "linbox/algorithms/blackbox-block-container.h" diff --git a/linbox/algorithms/matrix-hom.h b/linbox/algorithms/matrix-hom.h index cf2ce423d..5d76f6454 100644 --- a/linbox/algorithms/matrix-hom.h +++ b/linbox/algorithms/matrix-hom.h @@ -46,7 +46,7 @@ namespace LinBox { - /// \brief Limited doc so far. Used in RationalSolver. + /// \brief Limited doc so far. Used in DixonSolver. namespace MatrixHom { // function class to hanle map to BlasMatrix (needed to allow partial specialization) @@ -235,7 +235,7 @@ namespace LinBox } template - void map (SparseMatrix &Ap, + void map (SparseMatrix &Ap, const SparseMatrix & A) { typedef SparseMatrix FMatrix; @@ -244,7 +244,7 @@ namespace LinBox } template - void map (SparseMatrix &Ap, + void map (SparseMatrix &Ap, const SparseMatrix & A) { typedef SparseMatrix FMatrix; diff --git a/linbox/algorithms/matrix-inverse.h b/linbox/algorithms/matrix-inverse.h index c0572c596..3f8e00a95 100644 --- a/linbox/algorithms/matrix-inverse.h +++ b/linbox/algorithms/matrix-inverse.h @@ -28,6 +28,7 @@ #ifndef __LINBOX_matrix_inverse_H #define __LINBOX_matrix_inverse_H +#include "linbox/field/multimod-field.h" #include "linbox/util/debug.h" #include "linbox/util/error.h" #include diff --git a/linbox/algorithms/rational-solver-adaptive.h b/linbox/algorithms/rational-solver-adaptive.h index 17a21b566..7f2f3efdd 100644 --- a/linbox/algorithms/rational-solver-adaptive.h +++ b/linbox/algorithms/rational-solver-adaptive.h @@ -45,13 +45,13 @@ namespace LinBox linbox_check ((M. rowdim() == M. coldim()) && (b.size() == M.rowdim()) && (num. size() ==M.coldim())); typedef Givaro::Modular Field; // typedef Givaro::Modular Field; - RationalSolver, Method::SymbolicNumericNorm> numerical_solver; - //RationalSolver, Method::SymbolicNumericOverlap> numerical_solver; + DixonSolver, Method::SymbolicNumericNorm> numerical_solver; + //DixonSolver, Method::SymbolicNumericOverlap> numerical_solver; SolverReturnStatus ret; ret = numerical_solver. solve(num, den, M, b); if (ret != SS_OK) { - RationalSolver> solver; + DixonSolver> solver; BlasVector Ib(M.field()); Ib.reserve(b.size()); typename IRing::Element tmp; for(typename InVector::const_iterator biter = b.begin(); @@ -74,12 +74,12 @@ namespace LinBox linbox_check ((M. rowdim() == M. coldim()) && (b.size() == M.rowdim()) && (num. size() ==M.coldim())); typedef Givaro::Modular Field; // typedef Givaro::Modular Field; - RationalSolver, Method::SymbolicNumericOverlap> numerical_solver; + DixonSolver, Method::SymbolicNumericOverlap> numerical_solver; SolverReturnStatus ret; ret = numerical_solver. solve(num, den, M, b); if (ret != SS_OK) { - RationalSolver > solver; + DixonSolver > solver; ret = solver. solve(num, den, M, b); } diff --git a/linbox/algorithms/rational-solver.h b/linbox/algorithms/rational-solver.h index 8dc929f58..5ef2a37a3 100644 --- a/linbox/algorithms/rational-solver.h +++ b/linbox/algorithms/rational-solver.h @@ -125,13 +125,13 @@ namespace LinBox * - Ring: ring over which entries are defined * - Field: finite field for p-adic lifting * - RandomPrime: generator of random primes - * - MethodTraits: type of subalgorithm to use in p-adic lifting (default is Method::Dixon) + * - MethodTraits: type of subalgorithm to use in p-adic lifting (default is Method::DenseElimination) * . * * \ingroup padic */ - template - class RationalSolver { + template + class DixonSolver { public: /** Solve a linear system \c Ax=b over quotient field of a ring @@ -236,7 +236,7 @@ namespace LinBox * */ template - class RationalSolver { + class DixonSolver { public: typedef Ring RingType; @@ -265,7 +265,7 @@ namespace LinBox * @param rp a RandomPrime generator, set by default * @param traits */ - RationalSolver (const Ring& r = Ring(), + DixonSolver (const Ring& r = Ring(), const RandomPrime& rp = RandomPrime(), const Method::Wiedemann& traits=Method::Wiedemann()) : _ring(r), _genprime(rp), _traits(traits) @@ -285,7 +285,7 @@ namespace LinBox * @param rp a RandomPrime generator, set by default * @param traits */ - RationalSolver (const Prime& p, const Ring& r = Ring(), + DixonSolver (const Prime& p, const Ring& r = Ring(), const RandomPrime& rp = RandomPrime(), const Method::Wiedemann& traits=Method::Wiedemann()) : _ring(r), _genprime(rp), _prime(p), _traits(traits) @@ -444,7 +444,7 @@ namespace LinBox * */ template - class RationalSolver { + class DixonSolver { public: typedef Ring RingType; @@ -475,7 +475,7 @@ namespace LinBox * @param rp a RandomPrime generator, set by default * @param traits */ - RationalSolver (const Ring& r = Ring(), + DixonSolver (const Ring& r = Ring(), const RandomPrime& rp = RandomPrime(), const Method::BlockWiedemann& traits=Method::BlockWiedemann()) : _ring(r), _genprime(rp), _traits(traits) @@ -495,7 +495,7 @@ namespace LinBox * @param rp a RandomPrime generator, set by default * @param traits */ - RationalSolver (const Prime& p, const Ring& r = Ring(), + DixonSolver (const Prime& p, const Ring& r = Ring(), const RandomPrime& rp = RandomPrime(), const Method::BlockWiedemann& traits=Method::BlockWiedemann()) : _ring(r), _genprime(rp), _prime(p), _traits(traits) @@ -574,336 +574,15 @@ namespace LinBox } #endif }; // end of specialization for the class RationalSover with BlockWiedemann traits +} /*-------*/ - /* DIXON */ + /* Dense */ /*-------*/ -#ifdef RSTIMING - class DixonTimer { - public: - mutable Timer ttSetup, ttRecon, ttGetDigit, ttGetDigitConvert, ttRingApply, ttRingOther; - mutable int rec_elt; - void clear() const - { - ttSetup.clear(); - ttRecon.clear(); - ttGetDigit.clear(); - ttGetDigitConvert.clear(); - ttRingOther.clear(); - ttRingApply.clear(); - rec_elt=0; - } - - template - void update(RR& rr, LC& lc) const - { - ttSetup += lc.ttSetup; - ttRecon += rr.ttRecon; - rec_elt += rr._num_rec; - ttGetDigit += lc.ttGetDigit; - ttGetDigitConvert += lc.ttGetDigitConvert; - ttRingOther += lc.ttRingOther; - ttRingApply += lc.ttRingApply; - } - }; -#endif - - /** \brief partial specialization of p-adic based solver with Dixon algorithm. - * - * See the following reference for details on this algorithm: - * @bib - * - John D. Dixon Exact Solution of linear equations using p-adic - * expansions. Numerische Mathematik, volume 40, pages 137-141, - * 1982. - * - */ - - template - class RationalSolver { - - public: - - typedef Ring RingType; - typedef typename Ring::Element Integer; - typedef typename Field::Element Element; - typedef typename RandomPrime::Prime_Type Prime; - - // polymorphic 'certificate' generated when level >= SL_CERTIFIED - // certificate of inconsistency when any solving routine returns SS_INCONSISTENT - // certificate of minimal denominator when findRandomSolutionAndCertificate is called & return is SS_OK - mutable VectorFraction lastCertificate; - - //next 2 fields generated only by findRandomSolutionAndCertificate, when return is SS_OK - mutable Integer lastZBNumer; //filled in if level >= SL_CERTIFIED - mutable Integer lastCertifiedDenFactor; //filled in if level >= SL_LASVEGAS - //note: lastCertificate * b = lastZBNumer / lastCertifiedDenFactor, in lowest form - - protected: - - mutable RandomPrime _genprime; - mutable Prime _prime; - Ring _ring; -#ifdef RSTIMING - mutable Timer - tSetup, ttSetup, - tLQUP, ttLQUP, - tFastInvert, ttFastInvert, //only done in deterministic or inconsistent - tCheckConsistency,ttCheckConsistency, //includes lifting the certificate - tMakeConditioner, ttMakeConditioner, - tInvertBP, ttInvertBP, //only done in random - tCheckAnswer, ttCheckAnswer, - tCertSetup, ttCertSetup, //remaining 3 only done when makeMinDenomCert = true - tCertMaking, ttCertMaking, - - tNonsingularSetup,ttNonsingularSetup, - tNonsingularInv, ttNonsingularInv, - - totalTimer; - - mutable DixonTimer - ttConsistencySolve, ttSystemSolve, ttCertSolve, ttNonsingularSolve; -#endif - - public: - - /** Constructor - * @param r a Ring, set by default - * @param rp a RandomPrime generator, set by default - */ - RationalSolver (const Ring& r = Ring(), - const RandomPrime& rp = RandomPrime()) : - lastCertificate(r, 0), _genprime(rp), _ring(r) - { - _genprime.setBits(FieldTraits::bestBitSize()); - _prime=*_genprime; ++_genprime; -#ifdef RSTIMING - clearTimers(); -#endif - } - - - /** Constructor, trying the prime p first - * @param p a Prime - * @param r a Ring, set by default - * @param rp a RandomPrime generator, set by default - */ - RationalSolver (const Prime& p, const Ring& r = Ring(), - const RandomPrime& rp = RandomPrime()) : - lastCertificate(r, 0), _genprime(rp), _prime(p), _ring(r) - { - _genprime.setBits(FieldTraits::bestBitSize()); -#ifdef RSTIMING - clearTimers(); -#endif - } - - - /** Solve a linear system \c Ax=b over quotient field of a ring. - * - * @param num Vector of numerators of the solution - * @param den The common denominator. 1/den * num is the rational solution of \c Ax=b. - * @param A Matrix of linear system - * @param b Right-hand side of system - * @param s - * @param maxPrimes maximum number of moduli to try - * @param level level of certification to be used - * - * @return status of solution. if \c (return != SS_FAILED), and \c (level >= SL_LASVEGAS), solution is guaranteed correct. - * \c SS_FAILED - all primes used were bad - * \c SS_OK - solution found. - * \c SS_INCONSISTENT - system appreared inconsistent. certificate is in \p lastCertificate if \c (level >= SL_CERTIFIED) - */ - template - SolverReturnStatus solve(Vector1& num, Integer& den, const IMatrix& A, - const Vector2& b, const bool s = false, - const int maxPrimes = DEFAULT_MAXPRIMES, - const SolverLevel level = SL_DEFAULT) const; - - /** overload so that the bool 'oldMatrix' argument is not accidentally set to true */ - template - SolverReturnStatus solve(Vector1& num, Integer& den, - const IMatrix& A, const Vector2& b, const int maxPrimes, - const SolverLevel level = SL_DEFAULT) const - { - return solve (num, den, A, b, false, maxPrimes, level); - } - - /** Solve a nonsingular, square linear system \c Ax=b over quotient field of a ring. - * - * @param num Vector of numerators of the solution - * @param den The common denominator. 1/den * num is the rational solution of Ax = b - * @param A Matrix of linear system (it must be square) - * @param b Right-hand side of system - * @param s unused - * @param maxPrimes maximum number of moduli to try - * - * @return status of solution : - * - \c SS_FAILED all primes used were bad; - * - \c SS_OK solution found, guaranteed correct; - * - \c SS_SINGULAR system appreared singular mod all primes. - * . - */ - template - SolverReturnStatus solveNonsingular(Vector1& num, Integer& den, const IMatrix& A, - const Vector2& b, bool s = false, - int maxPrimes = DEFAULT_MAXPRIMES) const; - - /** Solve a general rectangular linear system \c Ax=b over quotient field of a ring. - * If A is known to be square and nonsingular, calling solveNonsingular is more efficient. - * - * @param num Vector of numerators of the solution - * @param den The common denominator. 1/den * num is the rational solution of Ax = b - * @param A Matrix of linear system - * @param b Right-hand side of system - * @param maxPrimes maximum number of moduli to try - * @param level level of certification to be used - * - * @return status of solution. if (return != SS_FAILED), and (level >= SL_LASVEGAS), solution is guaranteed correct. - * - \c SS_FAILED all primes used were bad - * - \c SS_OK solution found. - * - \c SS_INCONSISTENT system appreared inconsistent. certificate is in \p lastCertificate if (level >= SL_CERTIFIED) - */ - template - SolverReturnStatus solveSingular(Vector1& num, Integer& den, const IMatrix& A, - const Vector2& b, - int maxPrimes = DEFAULT_MAXPRIMES, const SolverLevel level = SL_DEFAULT) const; - - /** Find a random solution of the general linear system \c Ax=b over quotient field of a ring. - * - * @param num Vector of numerators of the solution - * @param den The common denominator. 1/den * num is the rational solution of Ax = b. - * @param A Matrix of linear system - * @param b Right-hand side of system - * @param maxPrimes maximum number of moduli to try - * @param level level of certification to be used - * - * @return status of solution. if (return != SS_FAILED), and (level >= SL_LASVEGAS), solution is guaranteed correct. - * - \c SS_FAILED all primes used were bad - * - \c SS_OK solution found. - * - \c SS_INCONSISTENT system appreared inconsistent. certificate is in lastCertificate if (level >= SL_CERTIFIED) - */ - template - SolverReturnStatus findRandomSolution(Vector1& num, Integer& den, const IMatrix& A, - const Vector2& b, - int maxPrimes = DEFAULT_MAXPRIMES, - const SolverLevel level = SL_DEFAULT) const; - - /** Big solving routine to perform random solving and certificate generation. - * Same arguments and return as findRandomSolution, except - * - * @param num Vector of numerators of the solution - * @param den The common denominator. 1/den * num is the rational solution of Ax = b - * @param A - * @param b - * @param randomSolution parameter to determine whether to randomize or not (since solveSingular calls this function as well) - * @param makeMinDenomCert determines whether a partial certificate for the minimal denominator of a rational solution is made - * @param maxPrimes - * @param level - * - * When (randomSolution == true && makeMinDenomCert == true), - * - If (level == SL_MONTECARLO) this function has the same effect as calling findRandomSolution. - * - If (level >= SL_LASVEGAS && return == SS_OK), \c lastCertifiedDenFactor contains a certified factor of the min-solution's denominator. - * - If (level >= SL_CERTIFIED && return == SS_OK), \c lastZBNumer and \c lastCertificate are updated as well. - * - */ - template - SolverReturnStatus monolithicSolve (Vector1& num, Integer& den, const IMatrix& A, - const Vector2& b, - bool makeMinDenomCert, bool randomSolution, - int maxPrimes = DEFAULT_MAXPRIMES, - const SolverLevel level = SL_DEFAULT) const; - - Ring getRing() const - { - return _ring; - } - - void chooseNewPrime() const - { - ++_genprime; - _prime = *_genprime; - } - -#ifdef RSTIMING - void clearTimers() const - { - ttSetup.clear(); - ttLQUP.clear(); - ttFastInvert.clear(); - ttCheckConsistency.clear(); - ttMakeConditioner.clear(); - ttInvertBP.clear(); - ttCheckAnswer.clear(); - ttCertSetup.clear(); - ttCertMaking.clear(); - ttNonsingularSetup.clear(); - ttNonsingularInv.clear(); - - ttConsistencySolve.clear(); - ttSystemSolve.clear(); - ttCertSolve.clear(); - ttNonsingularSolve.clear(); - } - - public: - - inline std::ostream& printTime(const Timer& timer, const char* title, - std::ostream& os, const char* pref = "") const - { - if (&timer != &totalTimer) - totalTimer += timer; - if (timer.count() > 0) { - os << pref << title; - for (int i=strlen(title)+strlen(pref); i<28; i++) - os << ' '; - return os << timer << std::endl; - } - else - return os; - } - - inline std::ostream& printDixonTime(const DixonTimer& timer, const char* title, - std::ostream& os) const - { - if (timer.ttSetup.count() > 0) { - printTime(timer.ttSetup, "Setup", os, title); - printTime(timer.ttGetDigit, "Field Apply", os, title); - printTime(timer.ttGetDigitConvert, "Ring-Field-Ring Convert", os, title); - printTime(timer.ttRingApply, "Ring Apply", os, title); - printTime(timer.ttRingOther, "Ring Other", os, title); - printTime(timer.ttRecon, "Reconstruction", os, title); - os<<" number of elt recontructed: "< - class RationalSolver; + class DixonSolver; /** \brief solver using a hybrid Numeric/Symbolic computation. * @@ -929,7 +608,7 @@ namespace LinBox //template argument Field and RandomPrime are not used. //Keep it just for interface consistency. template - class RationalSolver ; + class DixonSolver ; /*--------------*/ /* BLOCK HANKEL */ @@ -939,7 +618,7 @@ namespace LinBox * NO DOC */ template - class RationalSolver { + class DixonSolver { public: typedef Ring RingType; typedef typename Ring::Element Integer; @@ -958,7 +637,7 @@ namespace LinBox * @param r a Ring, set by default * @param rp a RandomPrime generator, set by default */ - RationalSolver (const Ring& r = Ring(), + DixonSolver (const Ring& r = Ring(), const RandomPrime& rp = RandomPrime()) : _genprime(rp), _ring(r) { @@ -973,7 +652,7 @@ namespace LinBox * @param r a Ring, set by default * @param rp a RandomPrime generator, set by default */ - RationalSolver (const Prime& p, const Ring& r = Ring(), + DixonSolver (const Prime& p, const Ring& r = Ring(), const RandomPrime& rp = RandomPrime()) : _genprime(rp), _prime(p), _ring(r) { @@ -997,7 +676,7 @@ namespace LinBox * NO DOC */ template - class RationalSolver { + class DixonSolver { public: typedef Ring RingType; typedef typename Ring::Element Integer; @@ -1016,7 +695,7 @@ namespace LinBox * @param r a Ring, set by default * @param rp a RandomPrime generator, set by default */ - RationalSolver (const Ring& r = Ring(), + DixonSolver (const Ring& r = Ring(), const RandomPrime& rp = RandomPrime()) : _genprime(rp), _ring(r) { @@ -1031,7 +710,7 @@ namespace LinBox * @param r a Ring, set by default * @param rp a RandomPrime generator, set by default */ - RationalSolver (const Prime& p, const Ring& r = Ring(), + DixonSolver (const Prime& p, const Ring& r = Ring(), const RandomPrime& rp = RandomPrime()) : _genprime(rp), _prime(p), _ring(r) {} diff --git a/linbox/algorithms/rational-solver.inl b/linbox/algorithms/rational-solver.inl index 8950f5d86..2df6b331a 100644 --- a/linbox/algorithms/rational-solver.inl +++ b/linbox/algorithms/rational-solver.inl @@ -65,30 +65,6 @@ namespace LinBox { - - /*! @brief NO DOC ! - * @bug why is this hard coded ? - */ - template - inline bool checkBlasPrime(const Prime p) - { - return p < Prime(67108863); - } - - template<> - inline bool checkBlasPrime(const BlasVector > p) - { - bool tmp=true; - for (size_t i=0;i= integer(67108863)) { - tmp=false; - break; - } - - return tmp; - } - - // SPECIALIZATION FOR WIEDEMANN // note: if Vector1 != Vector2 compilation of solve or solveSingluar will fail (via an invalid call to sparseprecondition)! @@ -97,7 +73,7 @@ namespace LinBox template template SolverReturnStatus - RationalSolver::solve (Vector1& num, + DixonSolver::solve (Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, @@ -132,7 +108,7 @@ namespace LinBox template template SolverReturnStatus - RationalSolver::solveNonsingular( Vector1& num, + DixonSolver::solveNonsingular( Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, @@ -206,7 +182,7 @@ namespace LinBox template template SolverReturnStatus - RationalSolver:: solveSingular (Vector1& num, + DixonSolver:: solveSingular (Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, @@ -341,7 +317,7 @@ namespace LinBox template template void - RationalSolver::sparseprecondition (const Field& F, + DixonSolver::sparseprecondition (const Field& F, const IMatrix *A, Compose, Compose > > *&PAQ, const FMatrix *Ap, @@ -403,7 +379,6 @@ namespace LinBox } - // SPECIALIZATION FOR BLOCK WIEDEMANN // note: if Vector1 != Vector2 compilation of solve or solveSingluar will fail (via an invalid call to sparseprecondition)! @@ -412,7 +387,7 @@ namespace LinBox template template SolverReturnStatus - RationalSolver::solve (Vector1& num, Integer& den, + DixonSolver::solve (Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, const bool old, @@ -446,7 +421,7 @@ namespace LinBox template template SolverReturnStatus - RationalSolver::solveNonsingular (Vector1& num, + DixonSolver::solveNonsingular (Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, @@ -496,844 +471,6 @@ namespace LinBox // END OF SPECIALIZATION FOR BLOCK WIEDEMANN - - // SPECIALIZATION FOR DIXON - - template - template - SolverReturnStatus - RationalSolver::solve (Vector1& num, - Integer& den, - const IMatrix& A, - const Vector2& b, - const bool old, - int maxP, - const SolverLevel level ) const - { - SolverReturnStatus status; - int maxPrimes=maxP; - while (maxPrimes > 0) - { -#ifdef SKIP_NONSINGULAR - switch (SS_SINGULAR) -#else - switch (A.rowdim() == A.coldim() ? solveNonsingular(num, den,A,b,old,maxPrimes) : SS_SINGULAR) -#endif - { - - case SS_OK: -#ifdef DEBUG_DIXON - std::cout <<"nonsingular worked\n"; -#endif - return SS_OK; - break; - - case SS_SINGULAR: -#ifdef DEBUG_DIXON - std::cout<<"switching to singular\n"; -#endif - status = solveSingular(num, den,A,b,maxPrimes,level); - if (status != SS_FAILED) - return status; - break; - - case SS_FAILED: - //std::cout <<"nonsingular failed\n"; // BDS: in what sense is this part of the spec of solve? - break; - - default: - throw LinboxError ("Bad return value from solveNonsingular"); - - } - maxPrimes--; - if (maxPrimes > 0) chooseNewPrime(); - } - return SS_FAILED; - } - - - template - template - SolverReturnStatus - RationalSolver::solveNonsingular(Vector1& num, - Integer& den, - const IMatrix& A, - const Vector2& b, - bool oldMatrix, - int maxPrimes) const - { - - //std::cout<<"DIXON\n\n\n\n"; -#ifdef DEBUG_DIXON - std::cout << "entering nonsingular solver\n"; -#endif - int trials = 0, notfr; - - // history sensitive data for optimal reason - // static const IMatrix* IMP; - - BlasMatrix* FMP = NULL; - Field *F=NULL; - - do - { -#if 0 - if (trials == maxPrimes) return SS_SINGULAR; - if (trials != 0) chooseNewPrime(); - ++trials; -#endif -#ifdef DEBUG_DIXON - //std::cout << "_prime: "<<_prime<<"\n"; - std::cout<<"A:=\n"; - A.write(std::cout, Tag::FileFormat::Maple); - std::cout<<"b:=[\n"; - for (size_t i=0;i segfault - if (FMP != NULL) delete FMP; - - // IMP = &A; - - if (F != NULL) delete F; - - F= new Field (_prime); - - FMP = new BlasMatrix(*F, A.rowdim(),A.coldim()); - - MatrixHom::map (*FMP, A ); // use MatrixHom to reduce matrix PG 2005-06-16 -#if 0 - typename BlasMatrix::Iterator iter_p = FMP->Begin(); - typename IMatrix::ConstIterator iter = A.Begin(); - for (;iter != A.End();++iter,++iter_p) - F->init(*iter_p, _ring.convert(tmp,*iter)); -#endif - -#ifdef DEBUG_DIXON - std::cout<< "p = "; - F->write(std::cout) << std::endl; - std::cout<<" Ap:="; - FMP->write(std::cout, Tag::FileFormat::Maple) << std::endl; -#endif - - if (!checkBlasPrime(_prime)){ - if (FMP != NULL) delete FMP; - FMP = new BlasMatrix(*F, A.rowdim(),A.coldim()); - notfr = (int)MatrixInverse::matrixInverseIn(*F,*FMP); - } - else { - BlasMatrix *invA = new BlasMatrix(*F, A.rowdim(),A.coldim()); - BlasMatrixDomain BMDF(*F); -#ifdef RSTIMING - tNonsingularSetup.stop(); - ttNonsingularSetup += tNonsingularSetup; - tNonsingularInv.start(); -#endif - assert(FMP != NULL); - BMDF.invin(*invA, *FMP, notfr); //notfr <- nullity - // if (FMP != NULL) // useless - delete FMP; - FMP = invA; -#if 0 - std::cout << "notfr = " << notfr << std::endl; - std::cout << "inverse mod p: " << std::endl; - FMP->write(std::cout, *F); -#endif -#ifdef RSTIMING - tNonsingularInv.stop(); - ttNonsingularInv += tNonsingularInv; -#endif - } - } - else { -#ifdef RSTIMING - tNonsingularSetup.stop(); - ttNonsingularSetup += tNonsingularSetup; -#endif - notfr = 0; - } - } while (notfr); - -#ifdef DEBUG_DIXON - std::cout<<"Ainvmodp:="; - FMP->write(std::cout, Tag::FileFormat::Maple)< > LiftingContainer; - LiftingContainer lc(_ring, *F, A, *FMP, b, _prime); - RationalReconstruction re(lc); - if (!re.getRational(num, den,0)){ - delete FMP; - return SS_FAILED; - } -#ifdef RSTIMING - ttNonsingularSolve.update(re, lc); -#endif - if (F!=NULL) - delete F; - if (FMP != NULL) - delete FMP; - return SS_OK; - } - - template - template - SolverReturnStatus - RationalSolver::solveSingular (Vector1& num, - Integer& den, - const IMatrix& A, - const Vector2& b, - int maxPrimes, - const SolverLevel level) const - { - return monolithicSolve (num, den, A, b, false, false, maxPrimes, level); - } - - template - template - SolverReturnStatus - RationalSolver::findRandomSolution (Vector1& num, - Integer& den, - const IMatrix& A, - const Vector2& b, - int maxPrimes, - const SolverLevel level ) const - { - - return monolithicSolve (num, den, A, b, false, true, maxPrimes, level); - } - - - // Most solving is done by the routine below. - // There used to be one for random and one for deterministic, but they have been merged to ease with - // repeated code (certifying inconsistency, optimization are 2 examples) - - template - template - SolverReturnStatus - RationalSolver::monolithicSolve (Vector1& num, - Integer& den, - const IMatrix& A, - const Vector2& b, - bool makeMinDenomCert, - bool randomSolution, - int maxPrimes, - const SolverLevel level) const - { - - // if (level == SL_MONTECARLO && maxPrimes > 1) - // std::cout << "WARNING: Even if maxPrimes > 1, SL_MONTECARLO uses just one prime." << std::endl; -#if 0 - if (makeMinDenomCert && !randomSolution) - std::cout << "WARNING: Will not compute a certificate of minimal denominator deterministically." << std::endl; -#endif - if (makeMinDenomCert && level == SL_MONTECARLO) - std::cout << "WARNING: No certificate of min-denominality generated due to level=SL_MONTECARLO" << std::endl; - int trials = 0; - while (trials < maxPrimes){ - if (trials != 0) chooseNewPrime(); - ++trials; -#ifdef DEBUG_DIXON - std::cout << "_prime:= "<<_prime<<";\n"; -#endif -#ifdef RSTIMING - tSetup.start(); -#endif - - // typedef typename Field::Element Element; - // typedef typename Ring::Element Integer_t; - typedef DixonLiftingContainer, BlasMatrix > LiftingContainer; - - // checking size of system - linbox_check(A.rowdim() == b.size()); - - LinBox::integer tmp; - Field F (_prime); - BlasMatrixDomain BMDI(_ring); - BlasMatrixDomain BMDF(F); - BlasApply BAR(_ring); - MatrixDomain MD(_ring); - VectorDomain VDR(_ring); - - BlasMatrix A_check(A); // used to check answer later - - // TAS_xxx stands for Transpose Augmented System (A|b)t - // this provides a factorization (A|b) = TAS_Pt . TAS_Ut . TAS_Qt . TAS_Lt - // such that - // - TAS_Q . (A|b) . TAS_P has nonzero principal minors up to TAS_rank - // - TAS_P permutes b to the (TAS_rank)th column of A iff the system is inconsistent mod p - BlasMatrix* TAS_factors = new BlasMatrix(F, A.coldim()+1, A.rowdim()); - Hom Hmap(_ring, F); - - BlasMatrix Ap(F, A.rowdim(), A.coldim()); // @fixme Shouldn't Ap(F, A) work without map below? - MatrixHom::map(Ap, A); - for (size_t i=0;isetEntry(j,i, Ap.getEntry(i,j)); - - for (size_t i=0;isetEntry(A.coldim(),i, tmpe); - } - - // @note TAS_factors now holds (A|b)t - -#ifdef RSTIMING - tSetup.stop(); - ttSetup += tSetup; - tPLUQ.start(); -#endif - BlasPermutation TAS_P(TAS_factors->rowdim()) ; - BlasPermutation TAS_Q(TAS_factors->coldim()) ; - - PLUQMatrix* TAS_PLUQ = new PLUQMatrix(*TAS_factors,TAS_P,TAS_Q); - size_t TAS_rank = TAS_PLUQ->getRank(); - - // check consistency. note, getQ returns Qt. - // BlasPermutation TAS_P = TAS_PLUQ->getP(); - // BlasPermutation TAS_Qt = TAS_PLUQ->getQ(); - std::vector srcRow(A.rowdim()), srcCol(A.coldim()+1); - std::vector::iterator sri = srcRow.begin(), sci = srcCol.begin(); - for (size_t i=0; i BMDs(iDom); - BMDs.mulin_right(TAS_P, srcCol); - BMDs.mulin_right(TAS_Q, srcRow); - - // @note srcCol/srcRow hold permutations (indices of (A|b)t) - // @fixme FFLAS has something for this? - -#ifdef DEBUG_INC - std::cout << "Q takes (0 1 ...) to ("; - for (size_t i=0; i= SL_LASVEGAS) { // in monte carlo, we assume A is actually empty - typename BlasMatrix::Iterator iter = A_check.Begin(); - for (; aEmpty && iter != A_check.End(); ++iter) - aEmpty &= _ring.isZero(*iter); - } - if (aEmpty) { - for (size_t i=0; i= SL_CERTIFIED) { - lastCertificate.clearAndResize(b.size()); - _ring.assign(lastCertificate.numer[i], _ring.one); - } - return SS_INCONSISTENT; - } -#if 0 - // both A and b are all zero. - for (size_t i=0; i= SL_LASVEGAS) - _ring.assign(lastCertifiedDenFactor, _ring.one); - if (level == SL_CERTIFIED) { - _ring.assign(lastZBNumer, _ring.zero); - lastCertificate.clearAndResize(b.size()); - } - return SS_OK; - } - // so a was empty mod p but not over Z. - - continue; //try new prime - } - - BlasMatrix* Atp_minor_inv = NULL; - - // @fixme What is the goal of this? - if ((appearsInconsistent && level > SL_MONTECARLO) || randomSolution == false) { - // take advantage of the (PLUQ)t factorization to compute - // an inverse to the leading minor of (TAS_Q . (A|b) . TAS_P) -#ifdef RSTIMING - tFastInvert.start(); -#endif - Atp_minor_inv = new BlasMatrix(F, rank, rank); - - FFLAS::fassign(F, rank, rank, TAS_factors->getPointer(), TAS_factors->getStride(), Atp_minor_inv->getPointer(), Atp_minor_inv->getStride()); - FFPACK::ftrtri (F, FFLAS::FflasUpper, FFLAS::FflasNonUnit, rank, Atp_minor_inv->getPointer(), Atp_minor_inv->getStride()); - FFPACK::ftrtri (F, FFLAS::FflasLower, FFLAS::FflasUnit, rank, Atp_minor_inv->getPointer(), Atp_minor_inv->getStride()); - FFPACK::ftrtrm (F, FFLAS::FflasLeft, FFLAS::FflasNonUnit, rank, Atp_minor_inv->getPointer(), Atp_minor_inv->getStride()); - // FFPACK::PLUQtoInverseOfFullRankMinor(F, rank, TAS_factors->getPointer(), A.rowdim(), - // TAS_Qt.getPointer(), - // Atp_minor_inv->getPointer(), rank); -#ifdef RSTIMING - tFastInvert.stop(); - ttFastInvert += tFastInvert; -#endif - } - - delete TAS_PLUQ; - delete TAS_factors; - - if (appearsInconsistent && level <= SL_MONTECARLO) - return SS_INCONSISTENT; - - if (appearsInconsistent) { -#ifdef RSTIMING - tCheckConsistency.start(); -#endif - Givaro::ZRing Z; - BlasVector > zt(Z,rank); - for (size_t i=0; i At_minor(_ring, rank, rank); - for (size_t i=0; iwrite(std::cout << "Atp_minor_inv:" << std::endl);//, F); - std::cout << "zt: "; for (size_t i=0; i re(lc); - - Vector1 short_num(A.field(),rank); - Integer short_den; - - if (!re.getRational(short_num, short_den,0)) - return SS_FAILED; // dirty, but should not be called - // under normal circumstances -#ifdef RSTIMING - ttConsistencySolve.update(re, lc); - tCheckConsistency.start(); -#endif - VectorFraction cert(_ring, short_num. size()); - cert. numer = short_num; - cert. denom = short_den; - cert.numer.resize(b.size()); - _ring.subin(cert.numer[rank], cert.denom); - _ring.assign(cert.denom, _ring.one); - BMDI.mulin_left(cert.numer, TAS_P); -#ifdef DEBUG_INC - cert.write(std::cout << "cert:") << std::endl; -#endif - - bool certifies = true; //check certificate - BlasVector certnumer_A(_ring,A.coldim()); - BAR.applyVTrans(certnumer_A, A_check, cert.numer); - typename BlasVector::iterator cai = certnumer_A.begin(); - for (size_t i=0; certifies && i A_minor(_ring, rank, rank); // -- will have the full rank minor of A - BlasMatrix *Ap_minor_inv; // -- will have inverse mod p of A_minor - BlasMatrix *P = NULL, *B = NULL; // -- only used in random case - - if (!randomSolution) { - // use shortcut - transpose Atp_minor_inv to get Ap_minor_inv - Element _rtmp; - Ap_minor_inv = Atp_minor_inv; - for (size_t i=0; igetEntry(_rtmp, i, j); - Ap_minor_inv->setEntry(i, j, Ap_minor_inv->refEntry(j, i)); - Ap_minor_inv->setEntry(j, i, _rtmp); - } - // @note minor inv = L1\U1 - - // permute original entries into A_minor - // @note A_minor = Pt A Qt - for (size_t i=0; i= SL_LASVEGAS){ // @fixme makeMinDenomCert always false? - B = new BlasMatrix(_ring, rank, A.coldim()); - for (size_t i=0; irefEntry(i, j), A_check.getEntry(srcRow[i],j)); - } - // @note B = Pt A - } - else { // if randomSolution == true - P = new BlasMatrix(_ring, A.coldim(), rank); - B = new BlasMatrix(_ring, rank,A.coldim()); - BlasMatrix Ap_minor(F, rank, rank); - Ap_minor_inv = new BlasMatrix(F, rank, rank); - - LinBox::integer tmp2=0; - size_t maxBitSize = 0; - for (size_t i=0; irefEntry(i, j), A_check.getEntry(srcRow[i], j)); - _ring.convert(tmp2, A_check.getEntry(srcRow[i], j)); - maxBitSize = std::max(maxBitSize, tmp2.bitsize()); - } - // @note B = Pt A -#ifdef RSTIMING - bool firstLoop = true; -#endif - // prepare B to be preconditionned through BLAS matrix mul - MatrixApplyDomain > MAD(_ring,*B); - MAD.setup(2); // @fixme Useless? - - int nullity; - do { // O(1) loops of this preconditioner expected -#ifdef RSTIMING - if (firstLoop) - firstLoop = false; - else - tMakeConditioner.start(); -#endif - // compute P a n*r random matrix of entry in [0,1] - typename BlasMatrix::Iterator iter; - for (iter = P->Begin(); iter != P->End(); ++iter) { - if (rand() > RAND_MAX/2) - _ring.assign(*iter, _ring.one); - else - _ring.assign(*iter, _ring.zero); - } - - // compute A_minor = B.P -#if 0 - if (maxBitSize * log((double)A.coldim()) > 53) - MD.mul(A_minor, *B, *P); - else { - double *B_dbl= new double[rank*A.coldim()]; - double *P_dbl= new double[A.coldim()*rank]; - double *A_minor_dbl = new double[rank*rank]; - for (size_t i=0;igetEntry(i,j)); - _ring.convert(P_dbl[i+j*rank], P->getEntry(j,i)); - } - cblas_dgemm(CblasRowMajor, CblasNoTrans, - CblasNoTrans, - rank, rank, A.coldim(), 1, - B_dbl, A.coldim(), P_dbl, rank, 0,A_minor_dbl, rank); - - for (size_t i=0;i&)*Ap_minor_inv, (BlasMatrix&)Ap_minor, nullity); -#ifdef RSTIMING - tInvertBP.stop(); - ttInvertBP += tInvertBP; -#endif - } while (nullity > 0); - } - // Compute newb = (TAS_Q.b)[0..(rank-1)] - BlasVector newb(b); - BMDI.mulin_right(TAS_Q, newb); - newb.resize(rank); - - BlasMatrix BBA_minor(A_minor); -#if 0 - BlasMatrix BBA_inv(F,*Ap_minor_inv); - BlasMatrix BBA_minor(A_minor); - BlasMatrix BBA_inv(*Ap_minor_inv); - LiftingContainer lc(_ring, F, BBA_minor, BBA_inv, newb, _prime); -#endif - LiftingContainer lc(_ring, F, BBA_minor, *Ap_minor_inv, newb, _prime); - -#ifdef DEBUG_DIXON - std::cout<<"length of lifting: "< re(lc); - - Vector1 short_num(_ring,rank); Integer short_den; - - - if (!re.getRational(short_num, short_den,0)) - return SS_FAILED; // dirty, but should not be called - // under normal circumstances - -#ifdef RSTIMING - ttSystemSolve.update(re, lc); - tCheckAnswer.start(); -#endif - - VectorFraction answer_to_vf(_ring, short_num. size()); - answer_to_vf. numer = short_num; - answer_to_vf. denom = short_den; - - if (!randomSolution) { - // short_answer = TAS_P * short_answer - answer_to_vf.numer.resize(A.coldim()+1,_ring.zero); - BMDI.mulin_left(answer_to_vf.numer, TAS_P); - answer_to_vf.numer.resize(A.coldim()); - } - else { - // short_answer = P * short_answer - BlasVector newNumer(_ring,A.coldim()); - BAR.applyV(newNumer, *P, answer_to_vf.numer); - //BAR.applyVspecial(newNumer, *P, answer_to_vf.numer); - - answer_to_vf.numer = newNumer; - } - - // @fixme Debug only? - if (level >= SL_LASVEGAS) { //check consistency - - BlasVector A_times_xnumer(_ring,b.size()); - - BAR.applyV(A_times_xnumer, A_check, answer_to_vf.numer); - - Integer tmpi; - - typename Vector2::const_iterator ib = b.begin(); - typename BlasVector::iterator iAx = A_times_xnumer.begin(); - int thisrow = 0; - bool needNewPrime = false; - - for (; !needNewPrime && ib != b.end(); ++iAx, ++ib, ++thisrow) - if (!_ring.areEqual(_ring.mul(tmpi, *ib, answer_to_vf.denom), *iAx)) { - // should attempt to certify inconsistency now - // as in "if [A31 | A32]y != b3" of step (4) - needNewPrime = true; - } - - if (needNewPrime) { - delete Ap_minor_inv; - if (randomSolution) {delete P; delete B;} -#ifdef RSTIMING - tCheckAnswer.stop(); - ttCheckAnswer += tCheckAnswer; -#endif - continue; //go to start of main loop - } - } - - //answer_to_vf.toFVector(answer); - num = answer_to_vf. numer; - den = answer_to_vf. denom; -#ifdef RSTIMING - tCheckAnswer.stop(); - ttCheckAnswer += tCheckAnswer; -#endif - // @fixme makeMinDenomCert always used as false? - if (makeMinDenomCert && level >= SL_LASVEGAS) // && randomSolution - { - // To make this certificate we solve with the same matrix as to get the - // solution, except transposed. -#ifdef RSTIMING - tCertSetup.start(); -#endif - Integer _rtmp; - Element _ftmp; - for (size_t i=0; igetEntry(_ftmp, i, j); - Ap_minor_inv->setEntry(i, j, Ap_minor_inv->refEntry(j, i)); - Ap_minor_inv->setEntry(j, i, _ftmp); - } - - for (size_t i=0; i Z; - BlasVector > q(Z,rank); - typename BlasVector >::iterator q_iter; - - bool allzero; - do { - allzero = true; - for (q_iter = q.begin(); q_iter != q.end(); ++q_iter) { - if (rand() > RAND_MAX/2) { - _ring.assign((*q_iter), _ring.one); - allzero = false; - } - else - (*q_iter) = _ring.zero; - } - } while (allzero); -#ifdef RSTIMING - tCertSetup.stop(); - ttCertSetup += tCertSetup; -#endif - //LiftingContainer lc2(_ring, F, BBA_minor, BBA_inv, q, _prime); - LiftingContainer lc2(_ring, F, A_minor, *Ap_minor_inv, q, _prime); - - RationalReconstruction rere(lc2); - Vector1 u_num(_ring,rank); Integer u_den; - if (!rere.getRational(u_num, u_den,0)) return SS_FAILED; - -#ifdef RSTIMING - ttCertSolve.update(rere, lc2); - tCertMaking.start(); -#endif - // remainder of code does z <- denom(partial_cert . Mr) * partial_cert * Qt - VectorFraction u_to_vf(_ring, u_num.size()); - u_to_vf. numer = u_num; - u_to_vf. denom = u_den; - BlasVector uB(_ring,A.coldim()); - BAR.applyVTrans(uB, *B, u_to_vf.numer); - -#if 0 - std::cout << "BP: "; - A_minor.write(std::cout, _ring) << std::endl; - std::cout << "q: "; - for (size_t i=0; i z(_ring, b.size()); //new constructor - u_to_vf.numer.resize(A.rowdim()); - - BMDI.mul(z.numer, u_to_vf.numer, TAS_Q); - - z.denom = numergcd; - - // z.write(std::cout << "z: ") << std::endl; - - if (level >= SL_CERTIFIED) - lastCertificate.copy(z); - - // output new certified denom factor - Integer znumer_b, zbgcd; - VDR.dotprod(znumer_b, z.numer, b); - _ring.gcd(zbgcd, znumer_b, z.denom); - _ring.div(lastCertifiedDenFactor, z.denom, zbgcd); - - if (level >= SL_CERTIFIED) - _ring.div(lastZBNumer, znumer_b, zbgcd); -#ifdef RSTIMING - tCertMaking.stop(); - ttCertMaking += tCertMaking; -#endif - } - - delete Ap_minor_inv; - delete B; - - if (randomSolution) {delete P;} - - // done making certificate, lets blow this popstand - return SS_OK; - } - - std::cout << "ouch" << std::endl; - return SS_FAILED; //all primes were bad - } - - - /* * Specialization for Block Hankel method */ @@ -1341,7 +478,7 @@ namespace LinBox template template SolverReturnStatus - RationalSolver::solveNonsingular(Vector1& num, + DixonSolver::solveNonsingular(Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, @@ -1435,7 +572,7 @@ namespace LinBox template template SolverReturnStatus - RationalSolver::solve(Vector1& num, + DixonSolver::solve(Vector1& num, Integer& den, const IMatrix& A, const Vector2& b, @@ -1502,7 +639,7 @@ namespace LinBox } //end of namespace LinBox //BB : moved the following "guarded" code in a new file, verbatim : -#include "rational-solver2.h" +#include "./dixon-solver/dixon-solver-symbolic-numeric.h" #endif //__LINBOX_rational_solver_INL diff --git a/linbox/blackbox/apply.h b/linbox/blackbox/apply.h index 672bf9c5b..786174215 100644 --- a/linbox/blackbox/apply.h +++ b/linbox/blackbox/apply.h @@ -119,7 +119,7 @@ namespace LinBox } inline Vector& applyVTrans(Vector &y, - BlasMatrix &A, + const BlasMatrix &A, const Vector &x) const { diff --git a/linbox/solutions/methods.h b/linbox/solutions/methods.h index bc06e0b78..881de1959 100644 --- a/linbox/solutions/methods.h +++ b/linbox/solutions/methods.h @@ -219,6 +219,8 @@ namespace LinBox { // ----- For Dixon method. // @fixme SingularSolutionType::Deterministic fails with Dense Dixon SingularSolutionType singularSolutionType = SingularSolutionType::Random; + bool certifyMinimalDenominator = false; //!< Whether the solver should try to find a certificate + //! that the provided denominator is minimal. // ----- For random-based systems. size_t trialsBeforeFailure = LINBOX_DEFAULT_TRIALS_BEFORE_FAILURE; //!< Maximum number of trials before giving up. diff --git a/linbox/solutions/solve.h b/linbox/solutions/solve.h index 224215e38..e827167c6 100644 --- a/linbox/solutions/solve.h +++ b/linbox/solutions/solve.h @@ -80,8 +80,8 @@ namespace LinBox { * - Otherwise > Error * - Method::Dixon * - IntegerTag - * | - DenseMatrix > `RationalSolver<..., Method::Dixon>` - * | - SparseMatrix > `RationalSolver<..., Method::SparseElimination>` + * | - DenseMatrix > `DixonSolver<..., Method::DenseElimination>` + * | - SparseMatrix > `DixonSolver<..., Method::SparseElimination>` * | - Otherwise > Error * - Otherwise > Error * - Method::Blackbox > Method::Wiedemann @@ -107,7 +107,7 @@ namespace LinBox { * | - Otherwise > Error * - Otherwise > Error * - Method::SymbolicNumericNorm - * - IntegerTag > `RationalSolver<..., Method::SymbolicNumericNorm>` + * - IntegerTag > `DixonSolver<..., Method::SymbolicNumericNorm>` * - Otherwise > Error * * @param [out] x solution, can be a rational solution (vector of numerators and one denominator) diff --git a/linbox/solutions/solve/solve-dixon.h b/linbox/solutions/solve/solve-dixon.h index fb9ab4c73..3afddc0ca 100644 --- a/linbox/solutions/solve/solve-dixon.h +++ b/linbox/solutions/solve/solve-dixon.h @@ -23,16 +23,12 @@ #pragma once #include -#include // @todo Rename dixon-rational-solver? +#include #include #include namespace LinBox { namespace { - // @fixme Remove when RationalSolver is renamed - template - using DixonRationalSolver = RationalSolver; - template struct MethodForMatrix { using type = Method::Wiedemann; // For blackboxes @@ -40,7 +36,7 @@ namespace LinBox { template struct MethodForMatrix> { - using type = Method::Dixon; + using type = Method::DenseElimination; }; template @@ -64,13 +60,12 @@ namespace LinBox { using PrimeGenerator = PrimeIterator; PrimeGenerator primeGenerator(FieldTraits::bestBitSize(A.coldim())); - using Solver = - DixonRationalSolver::type>; + using Solver = DixonSolver::type>; Solver dixonSolve(A.field(), primeGenerator); // @fixme I'm still bit sad that we cannot use generically the function below, // just because RationalSolve<..., SparseElimination> has not the same - // API (i.e. no solveNonSingular) than RationalSolver<..., Dixon> + // API (i.e. no solveNonSingular) than DixonSolver<..., DenseElimination> int maxTrials = m.trialsBeforeFailure; SolverReturnStatus status = dixonSolve.solve(xNum, xDen, A, b, maxTrials); @@ -78,7 +73,8 @@ namespace LinBox { if (status == SS_INCONSISTENT) { throw LinboxMathInconsistentSystem("From Dixon method."); - } else if (status == SS_FAILED || status == SS_BAD_PRECONDITIONER) { + } + else if (status == SS_FAILED || status == SS_BAD_PRECONDITIONER) { throw LinboxError("From Dixon method."); } } @@ -99,7 +95,7 @@ namespace LinBox { using PrimeGenerator = PrimeIterator; PrimeGenerator primeGenerator(FieldTraits::bestBitSize(A.coldim())); - using Solver = DixonRationalSolver::type>; + using Solver = DixonSolver::type>; Solver dixonSolve(A.field(), primeGenerator); // Either A is known to be non-singular, or we just don't know yet. @@ -113,15 +109,13 @@ namespace LinBox { // Either A is known to be singular, or we just failed trying to solve it as non-singular. if (singular) { - SolverLevel level = (m.certifyInconsistency ? SL_LASVEGAS : SL_MONTECARLO); - if (m.singularSolutionType == SingularSolutionType::Diophantine) { + SolverLevel level = (m.certifyInconsistency ? SL_LASVEGAS : SL_MONTECARLO); DiophantineSolver diophantineSolve(dixonSolve); status = diophantineSolve.diophantineSolve(xNum, xDen, A, b, maxTrials, level); } else { - bool randomSolutionType = (m.singularSolutionType == SingularSolutionType::Random); - status = dixonSolve.monolithicSolve(xNum, xDen, A, b, false, randomSolutionType, maxTrials, level); + status = dixonSolve.monolithicSolve(xNum, xDen, A, b, m); } } @@ -129,7 +123,8 @@ namespace LinBox { if (status == SS_INCONSISTENT) { throw LinboxMathInconsistentSystem("From Dixon method."); - } else if (status == SS_FAILED || status == SS_BAD_PRECONDITIONER) { + } + else if (status == SS_FAILED || status == SS_BAD_PRECONDITIONER) { throw LinboxError("From Dixon method."); } } @@ -150,13 +145,12 @@ namespace LinBox { using PrimeGenerator = PrimeIterator; PrimeGenerator primeGenerator(FieldTraits::bestBitSize(A.coldim())); - using Solver = - DixonRationalSolver::type>; + using Solver = DixonSolver::type>; Solver dixonSolve(A.field(), primeGenerator); // @fixme I'm a bit sad that we cannot use generically the function above, // just because RationalSolve<..., SparseElimination> has not the same - // API (i.e. no solveNonSingular) than RationalSolver<..., Dixon> + // API (i.e. no solveNonSingular) than DixonSolver<..., DenseElimination> int maxTrials = m.trialsBeforeFailure; SolverReturnStatus status = dixonSolve.solve(xNum, xDen, A, b, maxTrials); @@ -164,7 +158,8 @@ namespace LinBox { if (status == SS_INCONSISTENT) { throw LinboxMathInconsistentSystem("From Dixon method."); - } else if (status == SS_FAILED || status == SS_BAD_PRECONDITIONER) { + } + else if (status == SS_FAILED || status == SS_BAD_PRECONDITIONER) { throw LinboxError("From Dixon method."); } } diff --git a/linbox/solutions/solve/solve-numeric-symbolic.h b/linbox/solutions/solve/solve-numeric-symbolic.h index b9c6e97b5..b2ca51333 100644 --- a/linbox/solutions/solve/solve-numeric-symbolic.h +++ b/linbox/solutions/solve/solve-numeric-symbolic.h @@ -100,7 +100,7 @@ namespace LinBox { using Field = Givaro::Modular; // @fixme Why not double? using PrimeGenerator = PrimeIterator; - RationalSolver rsolver(b.field()); + DixonSolver rsolver(b.field()); SolverReturnStatus status = rsolver.solve(xNum, xDen, A, b); diff --git a/tests/test-last-invariant-factor.C b/tests/test-last-invariant-factor.C index a4a069caf..f3eed1833 100644 --- a/tests/test-last-invariant-factor.C +++ b/tests/test-last-invariant-factor.C @@ -207,8 +207,8 @@ int main(int argc, char** argv) RandomDenseStream s1 (R, gen, n, iterations); - typedef RationalSolver, PrimeIterator > Solver; - // typedef RationalSolver, LinBox::PrimeIterator > Solver; + typedef DixonSolver, PrimeIterator > Solver; + // typedef DixonSolver, LinBox::PrimeIterator > Solver; typedef LastInvariantFactor LIF; diff --git a/tests/test-one-invariant-factor.C b/tests/test-one-invariant-factor.C index e6d654afe..998400ffc 100644 --- a/tests/test-one-invariant-factor.C +++ b/tests/test-one-invariant-factor.C @@ -218,8 +218,8 @@ int main(int argc, char** argv) commentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (5); - typedef RationalSolver, PrimeIterator > Solver; - // typedef RationalSolver, LinBox::PrimeIterator > Solver; + typedef DixonSolver, PrimeIterator > Solver; + // typedef DixonSolver, LinBox::PrimeIterator > Solver; typedef LastInvariantFactor LIF; typedef OneInvariantFactor OIF; diff --git a/tests/test-rational-solver.C b/tests/test-rational-solver.C index a2e3c2c62..c5ad46e37 100644 --- a/tests/test-rational-solver.C +++ b/tests/test-rational-solver.C @@ -93,7 +93,7 @@ bool testRandomSolve (const Ring& R, for(int i = 0; i < n; ++i) R.init (D[(size_t)i][(size_t)i], d[(size_t)i]); - typedef RationalSolver > RSolver; + typedef DixonSolver > RSolver; RSolver rsolver; BlasVector num(R,(size_t)n); @@ -144,7 +144,7 @@ int main(int argc, char** argv) using Ring = Givaro::IntegerDom; Ring R; Ring::RandIter gen(R); - Field F(101); + Field F(101); RandomDenseStream s1 (R, gen, n, (unsigned int)iterations), s2 (R, gen, n, (unsigned int)iterations); if (!testRandomSolve(R, F, s1, s2)) pass = false; diff --git a/tests/test-smith-form-binary.C b/tests/test-smith-form-binary.C index 938a764bc..e3349a333 100644 --- a/tests/test-smith-form-binary.C +++ b/tests/test-smith-form-binary.C @@ -74,7 +74,7 @@ int main(int argc, char** argv) PIR R; typedef Givaro::Modular Field; - typedef RationalSolver > Solver; + typedef DixonSolver > Solver; typedef LastInvariantFactor LIF; typedef OneInvariantFactor OIF; typedef SmithFormBinary > SF; @@ -126,7 +126,7 @@ int main(int argc, char** argv) RandomDenseStream s1 (R, gen, n, 1); typedef Givaro::Modular Field; - typedef RationalSolver > Solver; + typedef DixonSolver > Solver; typedef LastInvariantFactor LIF; typedef OneInvariantFactor OIF; typedef SmithFormBinary > SF; diff --git a/tests/test-solve-nonsingular.C b/tests/test-solve-nonsingular.C index c73a56b43..e4468c042 100644 --- a/tests/test-solve-nonsingular.C +++ b/tests/test-solve-nonsingular.C @@ -387,7 +387,7 @@ int main(int argc, char** argv) { report << "zw: not done. Requires 64 bit architecture (maybe, needs checking -bds)." << std::endl << std::endl; } else { - RationalSolver, Method::SymbolicNumericNorm> rsolver(R); + DixonSolver, Method::SymbolicNumericNorm> rsolver(R); part_pass = testRandomSolve(R, rsolver, A, b); report << "zw: " << (part_pass ? "pass" : "fail") << std::endl << std::endl; } @@ -395,7 +395,7 @@ int main(int argc, char** argv) { pass = pass && part_pass; if(run & 4){ PrimeIterator genprime(FieldTraits::bestBitSize(A.coldim())); - RationalSolver, Method::Dixon> rsolver(R, genprime); + DixonSolver, Method::DenseElimination> rsolver(R, genprime); part_pass = testRandomSolve(R, rsolver, A, b); report << "dixon: " << (part_pass ? "pass" : "fail") << std::endl << std::endl; }