|
| 1 | +/* Copyright (C) The LinBox group |
| 2 | + * ========LICENCE======== |
| 3 | + * This file is part of the library LinBox. |
| 4 | + * |
| 5 | + * LinBox is free software: you can redistribute it and/or modify |
| 6 | + * it under the terms of the GNU Lesser General Public |
| 7 | + * License as published by the Free Software Foundation; either |
| 8 | + * version 2.1 of the License, or (at your option) any later version. |
| 9 | + * |
| 10 | + * This library is distributed in the hope that it will be useful, |
| 11 | + * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 12 | + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 13 | + * Lesser General Public License for more details. |
| 14 | + * |
| 15 | + * You should have received a copy of the GNU Lesser General Public |
| 16 | + * License along with this library; if not, write to the Free Software |
| 17 | + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
| 18 | + * ========LICENCE======== |
| 19 | + */ |
| 20 | + |
| 21 | +/**\file examples/dixonsolve.C |
| 22 | + @example examples/dixonsolve.C |
| 23 | + |
| 24 | + |
| 25 | + * \brief Dixon System Solving via Lifting using dense LU or sparse LU |
| 26 | + * \ingroup examples |
| 27 | + */ |
| 28 | +#include <iostream> |
| 29 | +#include <omp.h> |
| 30 | + |
| 31 | +#include "linbox/algorithms/rational-solver.h" |
| 32 | +#include "linbox/solutions/methods.h" |
| 33 | +#include "linbox/solutions/solve.h" |
| 34 | +#include "linbox/util/args-parser.h" |
| 35 | +#include "linbox/util/error.h" |
| 36 | +#include "linbox/util/matrix-stream.h" |
| 37 | +#include "linbox/randiter/random-prime.h" |
| 38 | +#include <givaro/givrandom.h> |
| 39 | + |
| 40 | +// DenseLU |
| 41 | +#include "linbox/matrix/dense-matrix.h" |
| 42 | + |
| 43 | +// SparseElim |
| 44 | +#include "linbox/matrix/sparse-matrix.h" |
| 45 | + |
| 46 | +using namespace LinBox; |
| 47 | +typedef Givaro::ZRing<Givaro::Integer> Ints; |
| 48 | +typedef DenseVector<Ints> ZVector; |
| 49 | + |
| 50 | + |
| 51 | +template<typename _Matrix, typename _EliminationMethod> |
| 52 | +int test(_Matrix A, std::string vector_file, bool inv, bool pp, bool sparse_elim) { |
| 53 | + // TODO : git rm dixondenseLU dixonsparseelim |
| 54 | + |
| 55 | + std::cout << "A is " << A.rowdim() << " by " << A.coldim() << std::endl; |
| 56 | + |
| 57 | + if (pp) |
| 58 | + { |
| 59 | + // Print Matrix |
| 60 | + |
| 61 | + // Matrix Market |
| 62 | + // std::cout << "A is " << A << std::endl; |
| 63 | + |
| 64 | + // Maple |
| 65 | + A.write(std::cout << "A:=", Tag::FileFormat::Maple) << ';' << std::endl; |
| 66 | + } |
| 67 | + |
| 68 | + // Vector File |
| 69 | + Ints ZZ; |
| 70 | + std::ifstream invect; |
| 71 | + |
| 72 | + ZVector B(ZZ, A.rowdim()); |
| 73 | + |
| 74 | + bool createB = vector_file.empty(); |
| 75 | + if (!createB) { |
| 76 | + invect.open (vector_file, std::ifstream::in); |
| 77 | + if (!invect) { |
| 78 | + createB = true; |
| 79 | + } else { |
| 80 | + for(ZVector::iterator it=B.begin(); it != B.end(); ++it) |
| 81 | + invect >> *it; |
| 82 | + } |
| 83 | + } |
| 84 | + |
| 85 | + // Vectors |
| 86 | + ZVector X(ZZ, A.coldim()); |
| 87 | + |
| 88 | + if (createB) { |
| 89 | + ZVector U(ZZ, A.coldim()); |
| 90 | + Givaro::GivRandom bgen( BaseTimer::seed() ); |
| 91 | + if (inv) { |
| 92 | + std::cerr << "Creating a random {-1,1} vector " << std::endl; |
| 93 | + for(auto& it:B) it = (bgen.brand()?1:-1); |
| 94 | + } else { |
| 95 | + std::cerr << "Creating a random consistant {-1,1} vector " << std::endl; |
| 96 | + for(FFPACK::rns_double::integer& it:U) it = (bgen.brand()?1:-1); |
| 97 | + |
| 98 | + // B = A U |
| 99 | + A.apply(B,U); |
| 100 | + } |
| 101 | + } |
| 102 | + |
| 103 | + if(pp) |
| 104 | + { |
| 105 | + // Print RHS |
| 106 | + B.write(std::cout << "B:=", Tag::FileFormat::Maple) << ';' << std::endl; |
| 107 | + } |
| 108 | + |
| 109 | + std::cout << "B is " << B.size() << "x1" << std::endl; |
| 110 | + |
| 111 | + Timer chrono; |
| 112 | + |
| 113 | + // BlasElimination |
| 114 | + // TODO : à vérifier si cela marche avec Method::DenseElimination |
| 115 | + _EliminationMethod M; |
| 116 | + if (inv){ |
| 117 | + M.singularity = Singularity::NonSingular; |
| 118 | + } |
| 119 | + |
| 120 | + |
| 121 | + //==================================================== |
| 122 | + // BEGIN Replacement solve with fixed prime |
| 123 | + |
| 124 | + Ints::Element d; |
| 125 | + Method::Dixon m(M); |
| 126 | + typedef Givaro::Modular<double> Field; |
| 127 | + |
| 128 | + const size_t bitsize((size_t) FieldTraits<Field>::bestBitSize(A.rowdim())); |
| 129 | + Givaro::Integer randomPrime( *(PrimeIterator<>(bitsize)) ); |
| 130 | + |
| 131 | + FixedPrimeIterator fixedprime( randomPrime ); |
| 132 | + DixonSolver<Ints, Field, FixedPrimeIterator, _EliminationMethod> rsolve(A.field(), fixedprime); |
| 133 | + std::cout << "Using: " << *fixedprime << " as the fixed p-adic." << std::endl; |
| 134 | + |
| 135 | + chrono.start(); |
| 136 | + if (!sparse_elim){ |
| 137 | + // Dense Elimination |
| 138 | + if (inv) |
| 139 | + { |
| 140 | + std::cout << "Solving using Dense Elimination for non singular system" << std::endl; |
| 141 | + SolverReturnStatus ss = rsolve.solveNonsingular(X, d, A, B); |
| 142 | + if (ss != SS_OK) { |
| 143 | + std::cerr << "Error during solveNonsingular (possibly singular matrix or p-adic precision too small)" << std::endl; |
| 144 | + exit(-1); |
| 145 | + } |
| 146 | + } |
| 147 | + else |
| 148 | + { |
| 149 | + std::cout << "Solving using Dense Elimination for any system" << std::endl; |
| 150 | + SolverReturnStatus ss = rsolve.solve(X, d, A, B); |
| 151 | + if (ss == SS_FAILED){ |
| 152 | + std::cerr << "Error during solve (all primes used were bad)" << std::endl; |
| 153 | + exit(-1); |
| 154 | + } |
| 155 | + if (ss == SS_INCONSISTENT){ |
| 156 | + std::cerr << "Error: system appeared inconsistent" << std::endl; |
| 157 | + exit(-1); |
| 158 | + } |
| 159 | + } |
| 160 | + } else { |
| 161 | + // Sparse Elimination |
| 162 | + try |
| 163 | + { |
| 164 | + std::cout << "Solving using Sparse Elimination for any system" << std::endl; |
| 165 | + rsolve.solve(X, d, A, B); |
| 166 | + } |
| 167 | + catch(LinboxError& e) |
| 168 | + { |
| 169 | + std::cerr << e << '\n'; |
| 170 | + exit(-1); |
| 171 | + } |
| 172 | + } |
| 173 | + chrono.stop(); |
| 174 | + |
| 175 | + std::cout << "CPU time (seconds): " << chrono.usertime() << std::endl; |
| 176 | + |
| 177 | + { |
| 178 | + // Solution size |
| 179 | + |
| 180 | + std::cout<<"Reduced solution: \n"; |
| 181 | + size_t maxbits=0; |
| 182 | + for (size_t i=0;i<A.coldim();++i){ |
| 183 | + maxbits=(maxbits > X[i].bitsize() ? maxbits: X[i].bitsize()); |
| 184 | + } |
| 185 | + std::cout<<" numerators of size "<<maxbits<<" bits" << std::endl |
| 186 | + <<" denominators hold over "<<d.bitsize()<<" bits\n"; |
| 187 | + } |
| 188 | + |
| 189 | + |
| 190 | + { |
| 191 | + // Check Solution |
| 192 | + |
| 193 | + VectorDomain<Ints> VD(ZZ); |
| 194 | + MatrixDomain<Ints> MD(ZZ); |
| 195 | + ZVector LHS(ZZ, A.rowdim()), RHS(ZZ, B); |
| 196 | + // check that Ax = d.b |
| 197 | + MD.vectorMul(LHS, A, X); |
| 198 | + VD.mulin(RHS, d); |
| 199 | + if (VD.areEqual(LHS, RHS)) |
| 200 | + std::cout << "Ax=d.b : Yes" << std::endl; |
| 201 | + else |
| 202 | + std::cout << "Ax=d.b : No" << std::endl; |
| 203 | + } |
| 204 | + |
| 205 | + { |
| 206 | + // Print Solution |
| 207 | + |
| 208 | + std::cout << "Solution is ["; |
| 209 | + for(auto it:X) ZZ.write(std::cout, it) << " "; |
| 210 | + std::cout << "] / "; |
| 211 | + ZZ.write(std::cout, d)<< std::endl; |
| 212 | + } |
| 213 | + return 0; |
| 214 | +} |
| 215 | + |
| 216 | +int main (int argc, char **argv) { |
| 217 | + |
| 218 | + // TODO : seed ? |
| 219 | + std::string matrix_file = ""; |
| 220 | + std::string vector_file = ""; |
| 221 | + bool inv = false; |
| 222 | + bool pp = false; |
| 223 | + bool sparse_elim = false; |
| 224 | + |
| 225 | + |
| 226 | + Argument as[] = { |
| 227 | + { 'm', "-m FILE", "Set the input file for the matrix.", TYPE_STR , &matrix_file }, |
| 228 | + { 'v', "-v FILE", "Set the input file for the vector.", TYPE_STR , &vector_file }, |
| 229 | + { 'i', "-i" , "whether the matrix is known to be invertible.", TYPE_BOOL , &inv }, |
| 230 | + { 'p', "-p" , "whether you want to pretty print the matrix.", TYPE_BOOL , &pp }, |
| 231 | + { 's', "-s" , "whether to use sparse elimination.", TYPE_BOOL , &sparse_elim }, |
| 232 | + END_OF_ARGUMENTS |
| 233 | + }; |
| 234 | + |
| 235 | + FFLAS::parseArguments(argc,argv,as); |
| 236 | + |
| 237 | + // Matrix File |
| 238 | + if (matrix_file.empty()) { |
| 239 | + std::cerr << "You must specify an input file for the matrix with -m" << std::endl; |
| 240 | + exit(-1); |
| 241 | + } |
| 242 | + std::ifstream input (matrix_file); |
| 243 | + if (!input) { std::cerr << "Error opening matrix file " << argv[1] << std::endl; exit(-1); } |
| 244 | + |
| 245 | + |
| 246 | + // Read Integral matrix from File |
| 247 | + Ints ZZ; |
| 248 | + MatrixStream< Ints > ms( ZZ, input ); |
| 249 | + |
| 250 | + if (sparse_elim){ |
| 251 | + SparseMatrix<Ints> A(ms); |
| 252 | + return test<SparseMatrix<Ints>,Method::SparseElimination> |
| 253 | + (A, vector_file, inv , pp, sparse_elim); |
| 254 | + } else { |
| 255 | + DenseMatrix<Ints> A(ms); |
| 256 | + return test<DenseMatrix<Ints>,Method::DenseElimination> |
| 257 | + (A, vector_file, inv , pp, sparse_elim); |
| 258 | + } |
| 259 | +} |
| 260 | + |
| 261 | + |
| 262 | +// Local Variables: |
| 263 | +// mode: C++ |
| 264 | +// tab-width: 4 |
| 265 | +// indent-tabs-mode: nil |
| 266 | +// c-basic-offset: 4 |
| 267 | +// End: |
| 268 | +// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s |
0 commit comments