1 /* Copyright (C) 2010 LinBox 2 * Written by JG Dumas 3 * 4 * 5 * 6 * ========LICENCE======== 7 * This file is part of the library LinBox. 8 * 9 * LinBox is free software: you can redistribute it and/or modify 10 * it under the terms of the GNU Lesser General Public 11 * License as published by the Free Software Foundation; either 12 * version 2.1 of the License, or (at your option) any later version. 13 * 14 * This library is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 * Lesser General Public License for more details. 18 * 19 * You should have received a copy of the GNU Lesser General Public 20 * License along with this library; if not, write to the Free Software 21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22 * ========LICENCE======== 23 */ 24 25 /*! @file algorithms/matrix-hom.h 26 * @ingroup algorithms 27 * @brief Matrix Homomorphism 28 * A map function converts a matrix on a field/ring 29 * to its natural image in another field/ring. 30 * @sa 31 * \c rebind operator. 32 */ 33 34 #ifndef __LINBOX_matrix_hom_H 35 #define __LINBOX_matrix_hom_H 36 37 //! @bug it is dangerous to include matrices defs that include hom for their rebind... 38 #include "linbox/integer.h" 39 #include "linbox/field/hom.h" 40 #include "linbox/matrix/matrix-category.h" 41 #include "linbox/matrix/dense-matrix.h" 42 #include "linbox/matrix/sparse-matrix.h" 43 #include "linbox/blackbox/blackbox.h" 44 #include "linbox/vector/blas-vector.h" 45 46 namespace LinBox 47 { 48 49 /// \brief Limited doc so far. Used in RationalSolver. 50 namespace MatrixHom 51 { 52 // function class to hanle map to BlasMatrix (needed to allow partial specialization) 53 template< class Field, class IMatrix, class Type> 54 class BlasMatrixMAP { 55 public: 56 template<class _Rep> 57 void operator() (BlasMatrix<Field,_Rep> &Ap, const IMatrix& A, Type type); 58 }; 59 60 template<class Field, class IMatrix> 61 class BlasMatrixMAP<Field, IMatrix, MatrixContainerCategory::Blackbox> { 62 public: 63 template<class _Rep> operator()64 void operator() (BlasMatrix<Field,_Rep> &Ap, const IMatrix &A, MatrixContainerCategory::Blackbox type) 65 { 66 67 68 typedef typename IMatrix::Field Ring; 69 Ring r = A.field(); 70 71 72 std::vector<typename Ring::Element> e(A.coldim(), r.zero), tmp(A.rowdim()); 73 74 typename BlasMatrix<Field,_Rep>::ColIterator col_p; 75 76 typename BlasMatrix<Field,_Rep>::Col::iterator elt_p; 77 78 typename std::vector<typename Ring::Element>::iterator e_p, tmp_p; 79 80 Hom<Ring, Field> hom(A. field(), Ap.field()); 81 82 for (col_p = Ap.colBegin(), e_p = e.begin(); 83 e_p != e.end(); ++ col_p, ++ e_p) { 84 85 r.assign(*e_p, r.one); 86 A.apply (tmp, e); 87 88 for (tmp_p = tmp.begin(), elt_p = col_p -> begin(); 89 tmp_p != tmp.end(); ++ tmp_p, ++ elt_p) 90 hom.image (*elt_p, *tmp_p); 91 r.assign(*e_p, r.zero); 92 } 93 } 94 }; 95 96 template<class Field, class IMatrix> 97 class BlasMatrixMAP<Field, IMatrix, MatrixContainerCategory::Container> { 98 public: 99 template<class _Rep> operator()100 void operator() (BlasMatrix<Field,_Rep> &Ap, const IMatrix &A, MatrixContainerCategory::Container type) 101 { 102 Hom<typename IMatrix::Field , Field> hom(A.field(), Ap.field()); 103 typename Field::Element e; 104 for( typename IMatrix::ConstIndexedIterator indices = A.IndexedBegin(); 105 (indices != A.IndexedEnd()) ; 106 ++indices ) { 107 108 hom. image (e, A.getEntry(indices.rowIndex(),indices.colIndex()) ); 109 110 if (!Ap.field().isZero(e)) 111 Ap.setEntry (indices.rowIndex(), 112 indices.colIndex(), e); 113 else 114 Ap.setEntry (indices.rowIndex(), 115 indices.colIndex(), Ap.field().zero); 116 } 117 118 } 119 }; 120 121 template<class Field, class IMatrix> 122 class BlasMatrixMAP<Field, IMatrix, MatrixContainerCategory::BlasContainer> { 123 public: 124 template<class _Rep> operator()125 void operator() (BlasMatrix<Field,_Rep> &Ap, const IMatrix &A, MatrixContainerCategory::BlasContainer type) 126 { 127 Hom<typename IMatrix::Field , Field> hom(A.field(), Ap.field()); 128 129 typename IMatrix::ConstIterator iterA = A.Begin(); 130 typename BlasMatrix<Field,_Rep>::Iterator iterAp = Ap.Begin(); 131 132 for(; iterA != A.End(); ++iterA, ++iterAp) 133 hom. image (*iterAp, *iterA); 134 } 135 }; 136 137 template<class Field, class _Rep> 138 class BlasMatrixMAP<Field, BlasMatrix<Givaro::ZRing<Integer>,_Rep>, MatrixContainerCategory::BlasContainer> { 139 public: 140 template<class _Rep2> operator()141 void operator() (BlasMatrix<Field,_Rep2> &Ap, const BlasMatrix<Givaro::ZRing<Integer>,_Rep> &A, MatrixContainerCategory::BlasContainer type) 142 { 143 Givaro::ZRing<Integer> ZZ ; 144 Hom<Givaro::ZRing<Integer> , Field> hom(ZZ, Ap.field()); 145 146 typename BlasMatrix<Givaro::ZRing<Integer>,_Rep>::ConstIterator iterA = A.Begin(); 147 typename BlasMatrix<Field,_Rep2>::Iterator iterAp = Ap.Begin(); 148 149 for(; iterA != A.End(); ++iterA, ++iterAp) 150 hom. image (*iterAp, *iterA); 151 } 152 }; 153 154 #ifdef __LINBOX_blas_matrix_multimod_H 155 template< class IMatrix> 156 class BlasMatrixMAP<MultiModDouble, IMatrix, MatrixContainerCategory::BlasContainer > { 157 public: 158 template<class _Rep> operator()159 void operator() (BlasMatrix<MultiModDouble,_Rep> &Ap, const IMatrix &A, MatrixContainerCategory::BlasContainer type) 160 { 161 for (size_t i=0; i<Ap.field().size();++i) 162 MatrixHom::map(Ap.getMatrix(i), A); 163 } 164 }; 165 166 template< class IMatrix> 167 class BlasMatrixMAP<MultiModDouble, IMatrix, MatrixContainerCategory::Container > { 168 public: 169 template<class _Rep> operator()170 void operator() (BlasMatrix<MultiModDouble,_Rep> &Ap, const IMatrix &A, MatrixContainerCategory::Container type) 171 { 172 for (size_t i=0; i<Ap.field().size();++i) 173 MatrixHom::map(Ap.getMatrix(i), A); 174 } 175 }; 176 177 template< class IMatrix> 178 class BlasMatrixMAP<MultiModDouble, IMatrix, MatrixContainerCategory::Blackbox > { 179 public: 180 template<class _Rep> operator()181 void operator() (BlasMatrix<MultiModDouble,_Rep> &Ap, const IMatrix &A, MatrixContainerCategory::Blackbox type) 182 { 183 for (size_t i=0; i<Ap.field().size();++i) 184 MatrixHom::map(Ap.getMatrix(i), A); 185 } 186 }; 187 #endif 188 } // MatrixHom 189 190 namespace MatrixHom 191 { 192 193 194 template<class FMatrix, class IMatrix> map(FMatrix & Ap,const IMatrix & A)195 void map (FMatrix & Ap, const IMatrix& A) 196 { 197 typename IMatrix::template rebind<typename FMatrix::Field>()( Ap, A); 198 } 199 200 // construct a sparse matrix over finite field, such that Ap = A mod p, where F = Ring / <p> 201 template<class Ring, class Vect1, class Field, class Vect2> map(SparseMatrix<Field,Vect2> & Ap,const SparseMatrix<Ring,Vect1> & A)202 void map (SparseMatrix<Field, Vect2>& Ap, const SparseMatrix<Ring, Vect1>& A) 203 { 204 // typename SparseMatrix<Ring,Vect1>::template rebind<Field, typename SparseVectorTranslate<Field,Vect2>::other_t >()( Ap, A); 205 typename SparseMatrix<Ring,Vect1>::template rebind<Field,Vect2>()( Ap, A); 206 } 207 208 209 // construct a BlasMatrix over finite field, such that Ap - A mod p, where F = Ring / <p> 210 template<class Ring, class Field, class _Rep> map(BlasMatrix<Field,_Rep> & Ap,const BlasMatrix<Ring,_Rep> & A)211 void map (BlasMatrix<Field,_Rep> &Ap, const BlasMatrix<Ring,_Rep>& A ) 212 { 213 typename BlasMatrix<Ring,_Rep>::template rebind<Field>()( Ap, A); 214 } 215 216 template <class Field, class IMatrix, class _Rep> map(BlasMatrix<Field,_Rep> & Ap,const IMatrix & A)217 void map (BlasMatrix<Field,_Rep> &Ap, const IMatrix &A) 218 { 219 BlasMatrixMAP<Field, IMatrix, typename MatrixContainerTrait<IMatrix>::Type> ()(Ap, A, typename MatrixContainerTrait<IMatrix>::Type()); 220 } 221 222 template <class Field, class IPoly, class IMatrix> map(PolynomialBB<typename IMatrix::template rebind<Field>::other,typename IPoly::template rebind<Field>::other> & Ap,const PolynomialBB<IMatrix,IPoly> & A)223 void map (PolynomialBB< typename IMatrix::template rebind<Field>::other, 224 typename IPoly::template rebind<Field>::other> &Ap, 225 const PolynomialBB<IMatrix, IPoly> &A) 226 { 227 typename PolynomialBB<IMatrix,IPoly>::template rebind<Field>() (Ap, A); 228 } 229 230 template <class Field, class Ring> map(ScalarMatrix<Field> & Ap,const ScalarMatrix<Ring> & A)231 void map (ScalarMatrix<Field> &Ap, 232 const ScalarMatrix<Ring> &A) 233 { 234 typename ScalarMatrix<Ring>::template rebind<Field>() (Ap, A); 235 } 236 237 template <class Field, class Field2, class Vect> map(SparseMatrix<Field,Vect> & Ap,const SparseMatrix<Field2,Vect> & A)238 void map (SparseMatrix<Field, Vect> &Ap, 239 const SparseMatrix<Field2, Vect> & A) 240 { 241 typedef SparseMatrix<Field, Vect> FMatrix; 242 typedef SparseMatrix<Field2,Vect> IMatrix; 243 typename IMatrix::template rebind<typename FMatrix::Field>()( Ap, A); 244 } 245 246 template <class Field, class Field2, class Vect> map(SparseMatrix<Field,Vect> & Ap,const SparseMatrix<Field,Vect> & A)247 void map (SparseMatrix<Field, Vect> &Ap, 248 const SparseMatrix<Field, Vect> & A) 249 { 250 typedef SparseMatrix<Field, Vect> FMatrix; 251 typedef SparseMatrix<Field2,Vect> IMatrix; 252 typename IMatrix::template rebind<typename FMatrix::Field>()( Ap, A); 253 } 254 255 256 // construct a sparse matrix over finite field, such that Ap = A mod p, where F = Ring / <p> 257 template <class Field, class Vect, class IMatrix> map(SparseMatrix<Field,Vect> & Ap,const IMatrix & A)258 void map (SparseMatrix<Field, Vect> &Ap, const IMatrix& A) 259 { 260 261 typedef typename IMatrix::Field Ring; 262 263 Ring r = A.field(); 264 265 BlasVector<Ring> e(r,A.coldim()), tmp(r,A.rowdim()); 266 267 typename BlasVector<Ring>::iterator iter, e_p; 268 269 typename Field::Element val; 270 271 int i = 0; 272 273 Hom<Ring, Field> hom(A. field(), Ap.field()); 274 275 for (e_p=e.begin();e_p != e.end(); ++e_p,++i){ 276 r.assign(*e_p, r.one); 277 A.apply(tmp,e); 278 int j; 279 for (iter=tmp.begin(),j=0; iter != tmp.end(); ++iter,++j) { 280 hom. image (val, *iter); 281 if (!Ap.field().isZero(val)) 282 Ap.setEntry ((size_t)j,(size_t)i, val); 283 284 } 285 r.assign(*e_p, r.zero); 286 } 287 288 } 289 290 291 } // MatrixHom 292 293 } // LinBox 294 295 #endif //__LINBOX_matrix_hom_H 296 297 // Local Variables: 298 // mode: C++ 299 // tab-width: 4 300 // indent-tabs-mode: nil 301 // c-basic-offset: 4 302 // End: 303 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 304