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