1 /* 2 Copyright (C) 2013 Tom Bachmann 3 4 This file is part of FLINT. 5 6 FLINT is free software: you can redistribute it and/or modify it under 7 the terms of the GNU Lesser General Public License (LGPL) as published 8 by the Free Software Foundation; either version 2.1 of the License, or 9 (at your option) any later version. See <http://www.gnu.org/licenses/>. 10 */ 11 12 // Common code shared among matrix classes 13 14 #ifndef FLINTXX_MATRIX_H 15 #define FLINTXX_MATRIX_H 16 17 #include "flint_classes.h" 18 #include "mp.h" 19 #include "rules.h" 20 #include "stdmath.h" 21 #include "ltuple.h" 22 #include "traits.h" 23 #include "tuple.h" 24 #include "../permxx.h" 25 26 namespace flint { 27 FLINT_DEFINE_BINOP(solve) 28 FLINT_DEFINE_BINOP(solve_fflu) 29 FLINT_DEFINE_THREEARY(mat_at) 30 FLINT_DEFINE_THREEARY(solve_fflu_precomp) 31 FLINT_DEFINE_UNOP(charpoly) 32 FLINT_DEFINE_UNOP(det) 33 FLINT_DEFINE_UNOP(det_fflu) 34 FLINT_DEFINE_UNOP(det_interpolate) 35 FLINT_DEFINE_UNOP(nullspace) 36 FLINT_DEFINE_UNOP(rref) 37 FLINT_DEFINE_UNOP(trace) 38 FLINT_DEFINE_UNOP(transpose) 39 40 FLINT_DEFINE_THREEARY(fflu) 41 FLINT_DEFINE_THREEARY_HERE_2DEFAULT(fflu, permxx*, 0, bool, false) 42 43 template<class T> struct matrix_traits : rules::UNIMPLEMENTED { }; 44 // override this for the non-reference type of your choice 45 // { 46 // template<class M> static slong rows(const M&); 47 // template<class M> static slong cols(const M&); 48 // template<class M> static ??? at(const M&, slong, slong); 49 // template<class M> static ??? at(M&, slong, slong); 50 // }; 51 namespace traits { 52 template <class T, class Enable = void> struct is_mat : is_implemented< 53 matrix_traits<typename flint_classes::to_nonref<T>::type> > { }; 54 template<class T> struct is_mat<T, 55 typename mp::disable_if<flint_classes::is_flint_class<T> >::type> 56 : mp::false_ { }; 57 } // traits 58 namespace matrices { 59 namespace mdetail { 60 // Some helper traits used below. 61 template<class Data> struct second_is_mat_data : mp::false_ { }; 62 template<class Data1, class Data2> 63 struct second_is_mat_data<tuple<Data1, tuple<Data2, empty_tuple> > > 64 : traits::is_mat<typename traits::basetype<Data2>::type> { }; 65 template<class Expr> struct second_is_mat 66 : second_is_mat_data<typename Expr::data_t> { }; 67 68 template<class Mat> 69 struct both_mat : mp::and_< 70 traits::is_mat< 71 typename traits::basetype<typename Mat::data_t::head_t>::type>, 72 traits::is_mat< 73 typename traits::basetype<typename Mat::data_t::tail_t::head_t>::type> 74 > { }; 75 76 // A more convenient way to obtain the traits associated to a non-immediate 77 // or non-nonref expression. 78 template<class Expr> struct immediate_traits 79 : matrix_traits<typename flint_classes::to_nonref<Expr>::type> { }; 80 } // mdetail 81 82 // For matrix expressions to create temporaries, it is necessary to know the 83 // dimensions of the result of a computation. This is a generic implementation, 84 // which assumes that the output dimensions are the same as the dimensions of 85 // the first argument, which is assumed to be a matrix, except if there are 86 // precisely two arguments only the second of which is a matrix, in which case 87 // we assume its the dimension of that. 88 // This implementation works correctly in many cases, e.g. matrix-addition or 89 // matrix-scalar multiplication. 90 template<class Operation> 91 struct outsize_generic 92 { 93 template<class Mat> 94 static slong rows(const Mat& m, 95 typename mp::disable_if<mdetail::second_is_mat<Mat> >::type* = 0) 96 { 97 return m._data().head.rows(); 98 } 99 template<class Mat> 100 static slong cols(const Mat& m, 101 typename mp::disable_if<mdetail::second_is_mat<Mat> >::type* = 0) 102 { 103 return m._data().head.cols(); 104 } 105 106 template<class Mat> 107 static slong rows(const Mat& m, 108 typename mp::enable_if<mdetail::second_is_mat<Mat> >::type* = 0) 109 { 110 return m._data().tail.head.rows(); 111 } 112 template<class Mat> 113 static slong cols(const Mat& m, 114 typename mp::enable_if<mdetail::second_is_mat<Mat> >::type* = 0) 115 { 116 return m._data().tail.head.cols(); 117 } 118 }; 119 120 // This is the expression template used for computing the dimensions of an 121 // operation. Without further specialisation, it is just the generic 122 // implementation described above. 123 // If you introduce a new operation where the generic implementation is 124 // incorrect, you must specialise this template. 125 template<class Operation> 126 struct outsize : outsize_generic<Operation> { }; 127 128 // Specialise immediates, where the dimensions are stored with the object. 129 template<> 130 struct outsize<operations::immediate> 131 { 132 template<class Mat> 133 static slong rows(const Mat& m) 134 { 135 return mdetail::immediate_traits<Mat>::rows(m); 136 } 137 template<class Mat> 138 static slong cols(const Mat& m) 139 { 140 return mdetail::immediate_traits<Mat>::cols(m); 141 } 142 }; 143 144 // Specialise multiplication. For matrix-matrix multiplication, use 145 // the usual formula. For matrix-scalar multiplication, use the generic 146 // implementation. 147 template<> 148 struct outsize<operations::times> 149 { 150 template<class Mat> 151 static slong rows(const Mat& m, 152 typename mp::enable_if<mdetail::both_mat<Mat> >::type* = 0) 153 { 154 return m._data().head.rows(); 155 } 156 template<class Mat> 157 static slong cols(const Mat& m, 158 typename mp::enable_if<mdetail::both_mat<Mat> >::type* = 0) 159 { 160 return m._data().tail.head.cols(); 161 } 162 163 template<class Mat> 164 static slong rows(const Mat& m, 165 typename mp::disable_if<mdetail::both_mat<Mat> >::type* = 0) 166 { 167 return outsize_generic<operations::times>::rows(m); 168 } 169 template<class Mat> 170 static slong cols(const Mat& m, 171 typename mp::disable_if<mdetail::both_mat<Mat> >::type* = 0) 172 { 173 return outsize_generic<operations::times>::cols(m); 174 } 175 }; 176 // Any particular multipication algorithm also has to be specialised. 177 template<> 178 struct outsize<operations::mul_classical_op> 179 : outsize<operations::times> { }; 180 181 // Specialise transpose. 182 template<> 183 struct outsize<operations::transpose_op> 184 { 185 template<class Mat> 186 static slong rows(const Mat& m) 187 { 188 return m._data().head.cols(); 189 } 190 template<class Mat> 191 static slong cols(const Mat& m) 192 { 193 return m._data().head.rows(); 194 } 195 }; 196 197 // Specialise nullspace. Note that the nullspace computation functions in 198 // flint return a matrix the columns of which span the nullspace. Since the 199 // nullity is not known in advance in general, we have to allocate a square 200 // matrix. 201 template<> 202 struct outsize<operations::nullspace_op> 203 { 204 template<class Mat> 205 static slong rows(const Mat& m) 206 { 207 return m._data().head.cols(); 208 } 209 template<class Mat> 210 static slong cols(const Mat& m) 211 { 212 return m._data().head.cols(); 213 } 214 }; 215 216 // This is a bit of a hack. Matrix operations returning a tuple typically 217 // only return one matrix. We key outsize on the inner operation to find out 218 // the dimensions. So e.g. solve(A, X).get<1>() (say) will invoke outsize 219 // with ltuple_get_op<1> as argument, which then invokes outsize with solve_op 220 // and (A, X) as argument. 221 template<unsigned n> 222 struct outsize<operations::ltuple_get_op<n> > 223 { 224 template<class Mat> 225 static slong rows(const Mat& m) 226 { 227 return outsize< 228 typename Mat::data_t::head_t::operation_t>::rows(m._data().head); 229 } 230 template<class Mat> 231 static slong cols(const Mat& m) 232 { 233 return outsize< 234 typename Mat::data_t::head_t::operation_t>::cols(m._data().head); 235 } 236 }; 237 238 // This is not actually a matrix expression, but called by the above ... 239 template<> 240 struct outsize<operations::solve_op> 241 { 242 template<class Mat> 243 static slong rows(const Mat& m) 244 { 245 return m._data().second().rows(); 246 } 247 template<class Mat> 248 static slong cols(const Mat& m) 249 { 250 return m._data().second().cols(); 251 } 252 }; 253 template<> struct outsize<operations::solve_fflu_op> 254 : outsize<operations::solve_op> { }; 255 template<> 256 struct outsize<operations::solve_fflu_precomp_op> 257 { 258 template<class Mat> 259 static slong rows(const Mat& m) 260 { 261 return m._data().tail.second().rows(); 262 } 263 template<class Mat> 264 static slong cols(const Mat& m) 265 { 266 return m._data().tail.second().cols(); 267 } 268 }; 269 270 namespace mdetail { 271 struct base_traits 272 { 273 template<class M> 274 static slong rows(const M& m) 275 { 276 return matrices::outsize<typename M::operation_t>::rows(m); 277 } 278 template<class M> 279 static slong cols(const M& m) 280 { 281 return matrices::outsize<typename M::operation_t>::cols(m); 282 } 283 }; 284 } // mdetail 285 286 // These traits classes are useful for implementing unified coefficient access. 287 // See fmpz_matxx etc for example usage. 288 template<class Mat> 289 struct generic_traits : mdetail::base_traits 290 { 291 template<class T, class U> 292 struct at 293 { 294 typedef FLINT_THREEARY_ENABLE_RETTYPE(mat_at, Mat, T, U) 295 entry_ref_t; 296 typedef entry_ref_t entry_srcref_t; 297 298 static entry_srcref_t get(const Mat& m, T i, U j) 299 { 300 return mat_at(m, i, j); 301 } 302 }; 303 }; 304 305 template<class Srcref> 306 struct generic_traits_srcref : mdetail::base_traits 307 { 308 template<class T, class U> 309 struct at 310 { 311 typedef Srcref entry_ref_t; 312 typedef Srcref entry_srcref_t; 313 314 template<class M> 315 static Srcref get(M m, T i, U j) 316 { 317 return mdetail::immediate_traits<M>::at(m, i, j); 318 } 319 }; 320 }; 321 322 template<class Ref> 323 struct generic_traits_ref : mdetail::base_traits 324 { 325 template<class T, class U> 326 struct at 327 { 328 typedef Ref entry_ref_t; 329 typedef Ref entry_srcref_t; 330 331 template<class M> 332 static Ref get(M m, T i, U j) 333 { 334 return mdetail::immediate_traits<M>::at(m, i, j); 335 } 336 }; 337 }; 338 339 template<class Ref, class Srcref> 340 struct generic_traits_nonref : mdetail::base_traits 341 { 342 template<class T, class U> 343 struct at 344 { 345 typedef Ref entry_ref_t; 346 typedef Srcref entry_srcref_t; 347 348 template<class M> 349 static Ref get(M& m, T i, U j) 350 { 351 return mdetail::immediate_traits<M>::at(m, i, j); 352 } 353 354 template<class M> 355 static Srcref get(const M& m, T i, U j) 356 { 357 return mdetail::immediate_traits<M>::at(m, i, j); 358 } 359 }; 360 }; 361 } // matrices 362 363 // immediate functions 364 template<class Mat> 365 inline typename mp::enable_if<traits::is_mat<Mat>, slong>::type 366 rank(const Mat& m) 367 { 368 return m.rank(); 369 } 370 371 template<class Mat> 372 inline typename mp::enable_if<traits::is_mat<Mat>, slong>::type 373 find_pivot_any(const Mat& m, slong start, slong end, slong c) 374 { 375 return m.find_pivot_any(start, end, c); 376 } 377 378 template<class Mat> 379 inline typename mp::enable_if<traits::is_mat<Mat>, slong>::type 380 find_pivot_partial(const Mat& m, slong start, slong end, slong c) 381 { 382 return m.find_pivot_partial(start, end, c); 383 } 384 } // flint 385 386 // Define rows(), cols(), create_temporary() and at() methods. 387 // For this to work, Traits must behave like the above generic traits. 388 // Also your matrix class must have a static create_temporary_rowscols function. 389 // See fmpz_mat for an example. 390 #define FLINTXX_DEFINE_MATRIX_METHODS(Traits) \ 391 template<class T, class U> \ 392 typename Traits::template at<T, U>::entry_ref_t at(T i, U j) \ 393 {return Traits::template at<T, U>::get(*this, i, j);} \ 394 template<class T, class U> \ 395 typename Traits::template at<T, U>::entry_srcref_t at(T i, U j) const \ 396 {return Traits::template at<T, U>::get(*this, i, j);} \ 397 \ 398 slong rows() const {return Traits::rows(*this);} \ 399 slong cols() const {return Traits::cols(*this);} \ 400 evaluated_t create_temporary() const \ 401 { \ 402 return create_temporary_rowscols(*this, rows(), cols()); \ 403 } 404 405 // Disable temporary merging. Requires create_temporary_rowscols. 406 // TODO do we really need the ltuple code everywhere? 407 #define FLINTXX_DEFINE_TEMPORARY_RULES(Matrix) \ 408 namespace traits { \ 409 template<> struct use_temporary_merging<Matrix> : mp::false_ { }; \ 410 } /* traits */ \ 411 namespace rules { \ 412 template<class Expr> \ 413 struct use_default_temporary_instantiation<Expr, Matrix> : mp::false_ { }; \ 414 template<class Expr> \ 415 struct instantiate_temporaries<Expr, Matrix, \ 416 typename mp::disable_if<flint_classes::is_Base<Matrix, Expr> >::type> \ 417 { \ 418 /* The only case where this should ever happen is if Expr is an ltuple */ \ 419 /* TODO static assert this */ \ 420 static Matrix get(const Expr& e) \ 421 { \ 422 typedef typename Expr::operation_t op_t; \ 423 return Matrix::create_temporary_rowscols(e, \ 424 matrices::outsize<op_t>::rows(e), \ 425 matrices::outsize<op_t>::cols(e)); \ 426 } \ 427 }; \ 428 template<class Expr> \ 429 struct instantiate_temporaries<Expr, Matrix, \ 430 typename mp::enable_if<flint_classes::is_Base<Matrix, Expr> >::type> \ 431 { \ 432 static Matrix get(const Expr& e) \ 433 { \ 434 return e.create_temporary(); \ 435 } \ 436 }; \ 437 } /* rules */ 438 439 // Add a fflu() member function to the matrix class. 440 #define FLINTXX_DEFINE_MEMBER_FFLU \ 441 template<class T> typename detail::nary_op_helper2<operations::fflu_op, \ 442 typename base_t::derived_t, permxx*, T>::enable::type \ 443 fflu(permxx* p, const T& t) const \ 444 { \ 445 return flint::fflu(*this, p, t); \ 446 } 447 448 #endif 449