00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_MATRIX_HPP
00014 #define __MMX_MATRIX_HPP
00015 #include <basix/iterator.hpp>
00016 #include <basix/vector.hpp>
00017 #include <algebramix/permutation.hpp>
00018 #include <algebramix/matrix_naive.hpp>
00019
00020 namespace mmx {
00021 #define Matrix_variant(C) matrix_variant_helper<C>::MV
00022 #define TMPL_DEF template<typename C, typename V= typename Matrix_variant(C)>
00023 #define TMPL template<typename C, typename V>
00024 #define TMPLK template<typename C, typename V, typename K>
00025 #define Format format<C>
00026 #define Matrix matrix<C,V>
00027 #define Matrix_rep matrix_rep<C,V>
00028 #define Abs_matrix matrix<Abs_type(C),V>
00029 #define Real_matrix matrix<Real_type(C),V>
00030 #define Center_matrix matrix<Center_type(C),V>
00031 #define Radius_matrix matrix<Radius_type(C),V>
00032 #define Lifted_matrix matrix<Lift_type(C)>
00033 #define Projected_matrix matrix<Project_type(C)>
00034 #define Reconstructed_matrix matrix<Reconstruct_type(C)>
00035 #define Truncated_matrix matrix<Truncate_type(C)>
00036 #define Completed_matrix matrix<Complete_type(C)>
00037 #define Evaluated_matrix matrix<Evaluate_type(C,K)>
00038 TMPL class matrix_rep;
00039 TMPL class matrix;
00040 TMPL inline nat cols (const Matrix& m);
00041 TMPL inline nat rows (const Matrix& m);
00042 TMPL inline bool is_a_scalar (const Matrix& m);
00043 TMPL inline bool is_non_scalar (const Matrix& m);
00044 TMPL inline C* tab (Matrix& m);
00045 TMPL inline const C* tab (const Matrix& m);
00046
00047
00048
00049
00050
00051 TMPL_DEF
00052 class matrix_rep REP_STRUCT_1(C) {
00053 C* a;
00054 nat nr;
00055 nat nc;
00056 bool scalar_flag;
00057 public:
00058 inline matrix_rep (C* a2, nat nr2, nat nc2, bool flag, const Format& fm):
00059 Format (fm), a (a2), nr (nr2), nc (nc2), scalar_flag (flag) {}
00060 inline ~matrix_rep () { mmx_delete<C> (a, aligned_size<C,V> (nr * nc)); }
00061 friend class Matrix;
00062 friend nat cols LESSGTR (const Matrix& m);
00063 friend nat rows LESSGTR (const Matrix& m);
00064 friend C* tab LESSGTR (Matrix& m);
00065 friend const C* tab LESSGTR (const Matrix& m);
00066 friend bool is_a_scalar LESSGTR (const Matrix& m);
00067 friend bool is_non_scalar LESSGTR (const Matrix& m);
00068 };
00069
00070 TMPL_DEF
00071 class matrix {
00072 INDIRECT_PROTO_2 (matrix, matrix_rep, C, V)
00073 public:
00074 typedef implementation<vector_linear,V> Vec;
00075 typedef implementation<matrix_linear,V> Mat;
00076 private:
00077 static C zero;
00078 public:
00079 matrix (const Format& fm= Format (no_format ())) {
00080 nat nc = Mat::def_rows;
00081 nat nr = Mat::def_cols;
00082 nat l = aligned_size<C,V> (nr * nc);
00083 C* tab= mmx_formatted_new<C> (l, fm);
00084 rep= new Matrix_rep (tab, nr, nc, false, fm);
00085 }
00086
00087 template<typename T>
00088 matrix (const T& c) {
00089 Format fm= get_format (as<C> (c));
00090 nat nr = Mat::init_rows;
00091 nat nc = Mat::init_cols;
00092 nat l = aligned_size<C,V> (nr * nc);
00093 C* tab= mmx_formatted_new<C> (l, fm);
00094 rep= new Matrix_rep (tab, nr, nc, true, fm);
00095 Mat::set (rep->a, as<C> (c), nr, nc);
00096 }
00097
00098 template<typename T>
00099 matrix (const T& c, const Format& fm) {
00100 nat nr = Mat::init_rows;
00101 nat nc = Mat::init_cols;
00102 nat l = aligned_size<C,V> (nr * nc);
00103 C* tab= mmx_formatted_new<C> (l, fm);
00104 rep= new Matrix_rep (tab, nr, nc, true, fm);
00105 Mat::set (rep->a, as<C> (c), nr, nc);
00106 }
00107
00108 template<typename T, typename W>
00109 matrix (const matrix<T,W>& m) {
00110 Format fm= as<Format > (CF(m));
00111 nat l= aligned_size<C,V> (rows(m) * cols(m));
00112 C* t= mmx_formatted_new<C> (l, fm);
00113 rep= new Matrix_rep (t, rows(m), cols(m), is_a_scalar (m), fm);
00114 Vec::cast (rep->a, tab(m), rows(m) * cols(m));
00115 }
00116
00117 template<typename T, typename W>
00118 matrix (const matrix<T,W>& m, const Format& fm) {
00119 nat l= aligned_size<C,V> (rows(m) * cols(m));
00120 C* t= mmx_formatted_new<C> (l, fm);
00121 rep= new Matrix_rep (t, rows(m), cols(m), is_a_scalar (m), fm);
00122 Vec::cast (rep->a, tab(m), rows(m) * cols(m));
00123 }
00124
00125 matrix (const C& x) {
00126 Format fm= get_format (x);
00127 nat nr = Mat::init_rows;
00128 nat nc = Mat::init_cols;
00129 nat l = aligned_size<C,V> (nr * nc);
00130 C* tab= mmx_formatted_new<C> (l, fm);
00131 rep= new Matrix_rep (tab, nr, nc, true, fm);
00132 Mat::set (rep->a, x, nr, nc);
00133 }
00134
00135 matrix (const C& x, nat nr, nat nc) {
00136 Format fm= get_format (x);
00137 nat l = aligned_size<C,V> (nr * nc);
00138 C* tab= mmx_formatted_new<C> (l, fm);
00139 rep= new Matrix_rep (tab, nr, nc, false, fm);
00140 Mat::set (rep->a, x, nr, nc);
00141 }
00142
00143 matrix (C* tab, nat nr, nat nc, const Format& fm) {
00144 rep= new Matrix_rep (tab, nr, nc, false, fm);
00145 }
00146
00147 matrix (C* tab, nat nr, nat nc, bool flag, const Format& fm) {
00148 rep= new Matrix_rep (tab, nr, nc, flag, fm);
00149 }
00150
00151 matrix (const iterator<C>& it) {
00152 Format fm= CF(it);
00153 nat i, j, nr= Mat::def_rows, nc= Mat::def_cols;
00154 if (nr==0) nr=1;
00155 nat l = aligned_size<C,V> (nr * nc);
00156 C* tab= mmx_formatted_new<C> (l, fm);
00157 for (i=0; i<nr; i++, ++it)
00158 for (j=0; j<nc; j++, ++it)
00159 tab[Mat::index (i, j, nr, nc)]= *it;
00160 rep= new Matrix_rep (tab, nc, nr, false, fm);
00161 }
00162
00163 inline const C& operator () (nat i, nat j) const {
00164 VERIFY (is_non_scalar (*this), "non-scalar matrix expected");
00165 VERIFY (i<rep->nr, "row out of range");
00166 VERIFY (j<rep->nc, "column out of range");
00167 return rep->a[Mat::index (i, j, rows (*this), cols(*this))]; }
00168 inline C& operator () (nat i, nat j) {
00169 VERIFY (is_non_scalar (*this), "non-scalar matrix expected");
00170 VERIFY (i<rep->nr, "row out of range");
00171 VERIFY (j<rep->nc, "column out of range");
00172 return rep->a[Mat::index (i, j, rows (*this), cols(*this))]; }
00173 inline C scalar () const {
00174 VERIFY (is_a_scalar (*this), "scalar matrix expected");
00175 return rep->a[0]; }
00176 inline C& scalar () {
00177 VERIFY (is_a_scalar (*this), "scalar matrix expected");
00178 return rep->a[0]; }
00179 };
00180 INDIRECT_IMPL_2 (matrix, matrix_rep, typename C, C, typename V, V)
00181 DEFINE_UNARY_FORMAT_2 (matrix)
00182
00183 TMPL C Matrix::zero= C(0);
00184 TMPL inline nat rows (const Matrix& m) { return m->nr; }
00185 TMPL inline nat cols (const Matrix& m) { return m->nc; }
00186 TMPL inline nat nbcol (const Matrix& m) {return cols(m);}
00187 TMPL inline nat nbrow (const Matrix& m) {return rows(m);}
00188 TMPL inline nat N (const Matrix& m) { return rows (m) * cols (m); }
00189 TMPL inline Format CF (const Matrix& m) { return m->tfm (); }
00190 TMPL inline bool is_a_scalar (const Matrix& m) { return m->scalar_flag; }
00191 TMPL inline bool is_non_scalar (const Matrix& m) { return !m->scalar_flag; }
00192 TMPL inline bool is_square_matrix (const Matrix& m) {
00193 return rows (m) == cols (m); }
00194 TMPL inline const C& read (const Matrix& m, nat i, nat j) { return m(i,j); }
00195 TMPL inline const C* tab (const Matrix& m) { return m->a; }
00196 TMPL inline C* tab (Matrix& m) {
00197 VERIFY (is_non_scalar (m), "non-scalar matrix expected");
00198 m.secure(); return m->a; }
00199 template<typename C, typename V, typename C2, typename V2> inline Matrix
00200 extend (const Matrix& m, const matrix<C2,V2>& n) {
00201 VERIFY (is_a_scalar (m), "scalar matrix expected");
00202 VERIFY (is_non_scalar (n), "non-scalar matrix expected");
00203 return Matrix (m.scalar(), rows (n), cols (n)); }
00204
00205 STYPE_TO_TYPE(TMPL,scalar_type,Matrix,C);
00206 UNARY_RETURN_TYPE(TMPL,abs_op,Matrix,Abs_matrix);
00207 UNARY_RETURN_TYPE(TMPL,Re_op,Matrix,Real_matrix);
00208 UNARY_RETURN_TYPE(TMPL,center_op,Matrix,Center_matrix);
00209 UNARY_RETURN_TYPE(TMPL,radius_op,Matrix,Radius_matrix);
00210 UNARY_RETURN_TYPE(TMPL,lift_op,Matrix,Lifted_matrix);
00211 UNARY_RETURN_TYPE(TMPL,project_op,Matrix,Projected_matrix);
00212 UNARY_RETURN_TYPE(TMPL,reconstruct_op,Matrix,Reconstructed_matrix);
00213 BINARY_RETURN_TYPE(TMPL,truncate_op,Matrix,nat,Truncated_matrix);
00214 UNARY_RETURN_TYPE(TMPL,complete_op,Matrix,Completed_matrix);
00215 BINARY_RETURN_TYPE(TMPLK,evaluate_op,Matrix,K,Evaluated_matrix);
00216
00217 TMPL
00218 struct fast_helper<Matrix > {
00219 typedef Fast_type(C) FC;
00220 typedef matrix<FC> fast_type;
00221 typedef typename Matrix_variant(FC) FV;
00222 static inline fast_type dd (const Matrix& m) {
00223 typedef implementation<vector_linear,FV> FVec;
00224 format<FC> fm= fast (CF(m));
00225 nat l= aligned_size<FC,FV> (rows (m) * cols (m));
00226 FC* r= mmx_formatted_new<FC> (l, fm);
00227 FVec::template vec_unary<fast_op> (r, tab (m), rows (m) * cols (m));
00228 return fast_type (r, rows (m), cols (m), is_a_scalar (m), fm); }
00229 static inline Matrix uu (const fast_type& m) {
00230 typedef implementation<vector_linear,V> Vec;
00231 Format fm= slow<C> (CF(m));
00232 nat l= aligned_size<C,V> (rows (m) * cols (m));
00233 C* r= mmx_formatted_new<C> (l, fm);
00234 Vec::template vec_unary<slow_op> (r, tab (m), rows (m) * cols (m));
00235 return Matrix (r, rows (m), cols (m), is_a_scalar (m), fm); }
00236 };
00237
00238 template<typename T, typename F, typename TV, typename FV>
00239 struct as_helper<matrix<T,TV>,matrix<F,FV> > {
00240 static inline matrix<T,TV>
00241 cv (const matrix<F,FV>& m) {
00242 format<T> fm= as<format<T> > (CF(m));
00243 nat n= rows(m) * cols(m);
00244 nat l= aligned_size<T,TV> (rows (m) * cols (m));
00245 const F* sm= tab (m);
00246 T* c= mmx_formatted_new<T> (l, fm);
00247 for (nat i=0; i<n; i++) c[i]= as<T> (sm[i]);
00248 return matrix<T,TV> (c, rows (m), cols (m), is_a_scalar (m), fm); }
00249 };
00250
00251 template<typename T, typename F, typename TV, typename FV> inline void
00252 set_as (matrix<T,TV>& r, const matrix<F,FV>& m) {
00253 r= matrix<T,TV> (m, CF(r));
00254 }
00255
00256 template<typename C, typename V, typename T> inline void
00257 set_as (Matrix& r, const T& x) {
00258 r= Matrix (x, CF(r));
00259 }
00260
00261
00262
00263
00264
00265 template<typename V, typename RS, typename CS>
00266 struct matrix_fixed: public V {
00267 typedef implementation<vector_linear,V> Vec;
00268 typedef matrix_fixed<typename V::Naive,RS,CS> Naive;
00269 typedef matrix_fixed<typename V::Positive,RS,CS> Positive;
00270 typedef matrix_fixed<typename V::No_simd,RS,CS> No_simd;
00271 typedef matrix_fixed<typename V::No_thread,RS,CS> No_thread;
00272 typedef matrix_fixed<typename V::No_scaled,RS,CS> No_scaled;
00273 };
00274
00275 template<typename F, typename V, typename RS, typename CS>
00276 struct implementation<F,matrix_fixed<V,RS,CS> >:
00277 public implementation<F,V> {};
00278
00279 template<typename V, typename RS, typename CS>
00280 struct implementation<matrix_defaults,matrix_fixed<V,RS,CS> > {
00281 static const nat def_rows = RS::val;
00282 static const nat def_cols = CS::val;
00283 static const nat init_rows= RS::val;
00284 static const nat init_cols= CS::val;
00285 };
00286
00287 #define FMatrix matrix<C,matrix_fixed<V,RS,CS> >
00288 #define FMatrix_rep matrix_rep<C,matrix_fixed<V,RS,CS> >
00289
00290 template<typename C, typename V, typename RS, typename CS>
00291 class FMatrix_rep REP_STRUCT_1(C) {
00292 C* a;
00293 public:
00294 inline matrix_rep (C* a2, nat, nat, bool, const Format& fm):
00295 Format (fm), a (a2) {}
00296 inline virtual ~matrix_rep () { mmx_delete<C> (a, RS::val * CS::val); }
00297 friend class FMatrix;
00298 friend nat cols LESSGTR (const FMatrix& m);
00299 friend nat rows LESSGTR (const FMatrix& m);
00300 friend C* tab LESSGTR (const FMatrix& m);
00301 friend bool is_a_scalar LESSGTR (const Matrix& m);
00302 friend bool is_non_scalar LESSGTR (const Matrix& m);
00303 };
00304
00305 template<typename C, typename V, typename RS, typename CS> inline nat
00306 cols (const FMatrix& m) { return CS::val; }
00307 template<typename C, typename V, typename RS, typename CS> inline nat
00308 rows (const FMatrix& m) { return RS::val; }
00309 template<typename C, typename V, typename RS, typename CS> inline bool
00310 is_a_scalar (const FMatrix& m) { (void) m; return false; }
00311 template<typename C, typename V, typename RS, typename CS> inline bool
00312 is_non_scalar (const FMatrix& m) { (void) m; return true; }
00313
00314 #undef FMatrix
00315 #undef FMatrix_rep
00316
00317
00318
00319
00320
00321 TMPL_DEF
00322 class matrix_iterator_rep: public iterator_rep<C> {
00323 public:
00324 typedef implementation<matrix_linear,V> Mat;
00325 private:
00326 const Matrix m;
00327 nat i;
00328 nat j;
00329 nat nc;
00330 protected:
00331 bool is_busy () { return i < rows(m); }
00332 void advance () { j++; if (j >= cols(m)) { j=0; i++; } }
00333 C current () {return m (i, j); }
00334 public:
00335 matrix_iterator_rep (const Matrix& m2, nat c = Mat::def_cols):
00336 iterator_rep<C> (CF(m)), m(m2), i(0), j(0), nc(c) {}
00337 };
00338
00339 TMPL iterator<C>
00340 iterate (const Matrix& m) {
00341 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
00342 return iterator<C> (new matrix_iterator_rep<C,V> (m));
00343 }
00344
00345
00346
00347
00348
00349 TMPL syntactic
00350 flatten (const Matrix& m) {
00351 if (is_a_scalar (m)) return flatten (m.scalar());
00352 int i, j, nr= rows(m), nc= cols(m);
00353 vector<syntactic> v;
00354 for (i=0; i<nr; i++) {
00355 vector<syntactic> h;
00356 for (j=0; j<nc; j++)
00357 h << flatten (m (i, j));
00358 v << apply (GEN_ROW, h);
00359 }
00360 return apply (GEN_SQTUPLE, v);
00361 }
00362
00363 TMPL
00364 struct binary_helper<Matrix >: public void_binary_helper<Matrix > {
00365 static inline string short_type_name () {
00366 return "M" * Short_type_name (C); }
00367 static inline generic full_type_name () {
00368 return gen ("Matrix", Full_type_name (C)); }
00369 static generic disassemble (const Matrix& m) {
00370 if (is_a_scalar (m)) return as<generic> (m.scalar ());
00371 vector<generic> v= vec (as<generic> (rows (m)), as<generic> (cols (m)));
00372 for (nat i=0; i<rows(m)*cols(m); i++)
00373 v << as<generic> (tab(m)[i]);
00374 return as<generic> (v); }
00375 static Matrix assemble (const generic& x) {
00376
00377 if (!is<vector<generic> > (x)) return Matrix (as<C> (x));
00378 vector<generic> v= as<vector<generic> > (x);
00379 nat rs= as<nat> (v[0]);
00380 nat cs= as<nat> (v[1]);
00381 nat l = aligned_size<C,V> (rs*cs);
00382 C* r = mmx_new<C> (l);
00383 for (nat i=0; i<rs*cs; i++)
00384 r[i]= as<C> (v[i+2]);
00385 return Matrix (r, rs, cs, Format ()); }
00386 static void write (const port& out, const Matrix& m) {
00387 binary_write<Format > (out, CF (m));
00388 binary_write<bool> (out, is_a_scalar (m));
00389 if (is_a_scalar (m)) binary_write<C> (out, m.scalar ());
00390 else {
00391 binary_write<nat> (out, rows(m));
00392 binary_write<nat> (out, cols(m));
00393 for (nat i=0; i<rows(m)*cols(m); i++)
00394 binary_write<C> (out, tab(m)[i]); } }
00395 static Matrix read (const port& in) {
00396 Format fm= binary_read<Format > (in);
00397 bool sc= binary_read<bool> (in);
00398 if (sc) return Matrix (binary_read<C> (in), fm);
00399 else {
00400 nat rs= binary_read<nat> (in);
00401 nat cs= binary_read<nat> (in);
00402 nat l = aligned_size<C,V> (rs*cs);
00403 C* r = mmx_formatted_new<C> (l, fm);
00404 for (nat i=0; i<rs*cs; i++)
00405 r[i]= binary_read<C> (in);
00406 return Matrix (r, rs, cs, fm); } }
00407 };
00408
00409
00410
00411
00412
00413 template<typename Op, typename C, typename V> nat
00414 unary_hash (const matrix<C,V>& m) {
00415 register nat i, j, h= 214365, nr= rows (m), nc= cols (m);
00416 for (i=0; i<nr; i++)
00417 for (j=0; j<nc; j++)
00418 h= (h<<1) ^ (h<<5) ^ (h>>27) ^ Op::op (m (i, j));
00419 return h;
00420 }
00421
00422 template<typename Op, typename C, typename V>
00423 matrix<Unary_return_type(Op,C),V>
00424 unary_map (const matrix<C,V>& m) {
00425 typedef implementation<vector_linear,V> Vec;
00426 typedef Unary_return_type(Op,C) T;
00427 format<T> fm= unary_map<Op> (CF(m));
00428 if (is_a_scalar (m)) return matrix<T,V> (Op::op (m.scalar()));
00429 nat nrows= rows (m);
00430 nat ncols= cols (m);
00431 nat l= aligned_size<T,V> (nrows * ncols);
00432 T* r= mmx_formatted_new<T> (l, fm);
00433 Vec::template vec_unary<Op> (r, tab (m), nrows*ncols);
00434 return matrix<T,V> (r, nrows, ncols, fm);
00435 }
00436
00437 template<typename Op, typename C1, typename C2, typename V>
00438 matrix<Binary_return_type(Op,C1,C2),V>
00439 binary_map (const matrix<C1,V>& m, const matrix<C2,V>& n) {
00440 typedef implementation<vector_linear,V> Vec;
00441 typedef Binary_return_type(Op,C1,C2) T;
00442 format<T> fm= binary_map<Op> (CF(m), CF(n));
00443 if (is_a_scalar (m) || is_a_scalar (n)) {
00444 if (is_non_scalar (m)) return binary_map<Op> (m, extend (n, m));
00445 if (is_non_scalar (n)) return binary_map<Op> (extend (m, n), n);
00446 return matrix<T,V> (Op::op (m.scalar(), n.scalar()));
00447 }
00448 nat nrows= rows (m);
00449 nat ncols= cols (m);
00450 ASSERT (rows (n) == nrows, "unequal number of rows");
00451 ASSERT (cols (n) == ncols, "unequal number of columns");
00452 nat l= aligned_size<T,V> (nrows * ncols);
00453 T* r= mmx_formatted_new<T> (l, fm);
00454 Vec::template vec_binary<Op> (r, tab (m), tab (n), nrows*ncols);
00455 return matrix<T,V> (r, nrows, ncols, fm);
00456 }
00457
00458 template<typename Op, typename C, typename V, typename X>
00459 matrix<Binary_return_type(Op,C,X),V>
00460 binary_map_scalar (const matrix<C,V>& m, const X& x) {
00461 typedef implementation<vector_linear,V> Vec;
00462 typedef Binary_return_type(Op,C,X) T;
00463 format<T> fm= binary_map_scalar<C> (CF(m), x);
00464 if (is_a_scalar (m)) return matrix<T,V> (Op::op (m.scalar(), x));
00465 nat nrows= rows (m);
00466 nat ncols= cols (m);
00467 nat l= aligned_size<T,V> (nrows * ncols);
00468 T* r= mmx_formatted_new<T> (l, fm);
00469 Vec::template vec_binary_scalar<Op> (r, tab (m), x, nrows*ncols);
00470 return matrix<T,V> (r, nrows, ncols, fm);
00471 }
00472
00473 template<typename Op, typename C, typename V> matrix<C,V>&
00474 nullary_set (matrix<C,V>& m) {
00475 typedef implementation<vector_linear,V> Vec;
00476 Vec::template vec_nullary<Op> (tab (m), rows (m) * cols (m));
00477 return m;
00478 }
00479
00480 template<typename Op, typename T, typename C, typename V> matrix<T,V>&
00481 unary_set (matrix<T,V>& m, const matrix<C,V>& n) {
00482 typedef implementation<vector_linear,V> Vec;
00483 if (is_a_scalar (m) || is_a_scalar (n)) {
00484 if (is_non_scalar (m))
00485 return unary_set<Op> (m, extend (n, m));
00486 else if (is_non_scalar (n))
00487 m= extend (m, n);
00488 else {
00489 Op::set_op (m.scalar(), n.scalar());
00490 return m;
00491 }
00492 }
00493 nat nrows= rows (m);
00494 nat ncols= cols (m);
00495 ASSERT (rows (n) == nrows, "unequal number of rows");
00496 ASSERT (cols (n) == ncols, "unequal number of columns");
00497 Vec::template vec_unary<Op> (tab (m), tab (n), nrows*ncols);
00498 return m;
00499 }
00500
00501 template<typename Op, typename T, typename V, typename X> matrix<T,V>&
00502 unary_set_scalar (matrix<T,V>& m, const X& x) {
00503 typedef implementation<vector_linear,V> Vec;
00504 if (is_a_scalar (m)) {
00505 Op::set_op (m.scalar(), x);
00506 return m;
00507 }
00508 nat nrows= rows (m);
00509 nat ncols= cols (m);
00510 Vec::template vec_unary_scalar<Op> (tab (m), x, nrows*ncols);
00511 return m;
00512 }
00513
00514 template<typename Op, typename C1, typename C2, typename V> bool
00515 binary_test (const matrix<C1,V>& m, const matrix<C2,V>& n) {
00516 typedef implementation<vector_linear,V> Vec;
00517 if (is_a_scalar (m) || is_a_scalar (n)) {
00518 if (is_non_scalar (m)) return binary_test<Op> (m, extend (n, m));
00519 if (is_non_scalar (n)) return binary_test<Op> (extend (m, n), n);
00520 return Op::op (m.scalar(), n.scalar());
00521 }
00522 nat nrows= rows (m);
00523 nat ncols= cols (m);
00524 if (rows (n) != nrows || cols (n) != ncols) return false;
00525 return Vec::template vec_binary_test<Op> (tab (m), tab (n), nrows*ncols);
00526 }
00527
00528 template<typename Op, typename C, typename V> Unary_return_type(Op,C)
00529 big (const matrix<C,V>& m) {
00530 typedef implementation<vector_linear,V> Vec;
00531 ASSERT (is_non_scalar (m), "non scalar matrix expected");
00532 return Vec::template vec_unary_big<Op> (tab (m), rows (m) * cols (m));
00533 }
00534
00535
00536
00537
00538
00539 template<typename S1, typename D, typename Fun> matrix<D>
00540 map (const Fun& fun, const matrix<S1>& m, const format<D>& fm) {
00541 nat n= rows(m) * cols(m);
00542 nat l= default_aligned_size<D> (n);
00543 D* r= mmx_formatted_new<D> (l, fm);
00544 const S1* src= tab (m);
00545 for (nat i=0; i<n; i++) r[i]= fun (src[i]);
00546 return matrix<D> (r, rows(m), cols(m), fm);
00547 }
00548
00549
00550
00551
00552
00553 TMPL void set_default (Matrix& m) { nullary_set<default_as_op> (m); }
00554 TMPL void set_pi (Matrix& m) { nullary_set<pi_as_op> (m); }
00555 TMPL void set_log2 (Matrix& m) { nullary_set<log2_as_op> (m); }
00556 TMPL void set_euler (Matrix& m) { nullary_set<euler_as_op> (m); }
00557 TMPL void set_catalan (Matrix& m) { nullary_set<catalan_as_op> (m); }
00558 TMPL void set_imaginary (Matrix& m) { nullary_set<imaginary_as_op> (m); }
00559 TMPL void set_nan (Matrix& m) { nullary_set<nan_as_op> (m); }
00560 TMPL void set_fuzz (Matrix& m) { nullary_set<fuzz_as_op> (m); }
00561 TMPL void set_smallest (Matrix& m) { nullary_set<smallest_as_op> (m); }
00562 TMPL void set_largest (Matrix& m) { nullary_set<largest_as_op> (m); }
00563 TMPL void set_accuracy (Matrix& m) { nullary_set<accuracy_as_op> (m); }
00564 TMPL void set_infinity (Matrix& m) { nullary_set<infinity_as_op> (m); }
00565 TMPL void set_maximal (Matrix& m) { nullary_set<maximal_as_op> (m); }
00566 TMPL void set_minimal (Matrix& m) { nullary_set<minimal_as_op> (m); }
00567
00568
00569
00570
00571
00572 TMPL Matrix copy (const Matrix& m) {
00573 return unary_map<id_op> (m); }
00574 TMPL Matrix duplicate (const Matrix& m) {
00575 return unary_map<duplicate_op> (m); }
00576
00577 TMPL Matrix operator - (const Matrix& m) {
00578 return unary_map<neg_op> (m); }
00579 TMPL Matrix operator + (const Matrix& m, const Matrix& n) {
00580 return binary_map<add_op> (m, n); }
00581 TMPL Matrix operator + (const Matrix& m, const C& n) {
00582 return binary_map<add_op> (m, Matrix (n)); }
00583 TMPL Matrix operator + (const C& m, const Matrix& n) {
00584 return binary_map<add_op> (Matrix (m), n); }
00585 TMPL Matrix operator - (const Matrix& m, const Matrix& n) {
00586 return binary_map<sub_op> (m, n); }
00587 TMPL Matrix operator - (const Matrix& m, const C& n) {
00588 return binary_map<sub_op> (m, Matrix (n)); }
00589 TMPL Matrix operator - (const C& m, const Matrix& n) {
00590 return binary_map<sub_op> (Matrix (m), n); }
00591 TMPL Matrix operator * (const Matrix& m, const C& c) {
00592 return binary_map_scalar<rmul_op> (m, c); }
00593 TMPL Matrix operator * (const C& c, const Matrix& m) {
00594 return binary_map_scalar<lmul_op> (m, c); }
00595 TMPL Matrix operator / (const Matrix& m, const C& c) {
00596 return binary_map_scalar<rdiv_op> (m, c); }
00597 TMPL Matrix quo (const Matrix& m, const C& c) {
00598 return binary_map_scalar<rquo_op> (m, c); }
00599 TMPL Matrix rem (const Matrix& m, const C& c) {
00600 return binary_map_scalar<rrem_op> (m, c); }
00601 TMPL Matrix& operator += (Matrix& m, const Matrix& n) {
00602 return unary_set<add_op> (m, n); }
00603 TMPL Matrix& operator -= (Matrix& m, const Matrix& n) {
00604 return unary_set<sub_op> (m, n); }
00605 TMPL Matrix& operator *= (Matrix& m, const C& x) {
00606 return unary_set_scalar<rmul_op> (m, x); }
00607 TMPL Matrix& operator /= (Matrix& m, const C& x) {
00608 return unary_set_scalar<rdiv_op> (m, x); }
00609 TMPL bool operator <= (const Matrix& m, const Matrix& n) {
00610 return binary_test<lesseq_op> (m, n); }
00611 TMPL bool operator >= (const Matrix& m, const Matrix& n) {
00612 return binary_test<gtreq_op> (m, n); }
00613 template<typename C, typename V, typename K> bool
00614 operator == (const Matrix& m, const K& c) {
00615 return m == Matrix (c, rows (m), cols (m)); }
00616 template<typename C, typename V, typename K> bool
00617 operator != (const Matrix& m, const K& c) {
00618 return m != Matrix (c, rows (m), cols (m)); }
00619 template<typename C, typename V, typename K> bool
00620 operator == (const K& c, const Matrix& m) {
00621 return m == Matrix (c, rows (m), cols (m)); }
00622 template<typename C, typename V, typename K> bool
00623 operator != (const K& c, const Matrix& m) {
00624 return m != Matrix (c, rows (m), cols (m)); }
00625 template<typename C, typename V, typename K> bool
00626 operator <= (const Matrix& m, const C& c) {
00627 return m <= Matrix (c, rows (m), cols (m)); }
00628 template<typename C, typename V, typename K> bool
00629 operator >= (const Matrix& m, const C& c) {
00630 return m >= Matrix (c, rows (m), cols (m)); }
00631
00632 TRUE_IDENTITY_OP_SUGAR(TMPL,Matrix)
00633 EXACT_IDENTITY_OP_SUGAR(TMPL,Matrix)
00634 ARITH_SCALAR_INT_SUGAR(TMPL,Matrix);
00635 STRICT_COMPARE_SUGAR(TMPL,Matrix)
00636
00637 TMPL Matrix derive (const Matrix& m) {
00638 return unary_map<derive_op> (m); }
00639 TMPL Matrix integrate (const Matrix& m) {
00640 return unary_map<integrate_op> (m); }
00641 template<typename C, typename V, typename K, typename W> Matrix
00642 integrate_init (const Matrix& m, const matrix<K,W>& n) {
00643 return binary_map<integrate_init_op> (m, as<matrix<K,V> > (n)); }
00644
00645 TMPL inline C big_add (const Matrix& m) { return big<add_op> (m); }
00646 TMPL inline C big_sup (const Matrix& m) { return big<sup_op> (m); }
00647 TMPL inline C big_max (const Matrix& m) { return big<max_op> (m); }
00648
00649
00650
00651
00652
00653 TMPL
00654 struct lift_helper<Matrix> {
00655 template<typename T,typename W> static void
00656 set_op (matrix<T,W>& n, const Matrix& m) {
00657 typedef implementation<vector_linear,V> Vec;
00658 format<T> fm= CF(n);
00659 if (is_a_scalar (m)) {
00660 T a= get_sample (fm);
00661 set_lift (a, m.scalar());
00662 n= matrix<T,W> (a, fm);
00663 } else {
00664 nat r= rows (m), c= cols (m);
00665 nat l= aligned_size<T,W> (r * c);
00666 T* t= mmx_formatted_new<T> (l, fm);
00667 Vec::template vec_unary<lift_op> (t, tab (m), r * c);
00668 n= matrix<T,W> (t, r, c, fm); } }
00669 static inline Lifted_matrix
00670 op (const Matrix& m) {
00671 Lifted_matrix n (lift (CF(m)));
00672 set_op (n, m);
00673 return n; }
00674 };
00675
00676 TMPL
00677 struct project_helper<Matrix> {
00678 template<typename T,typename W> static void
00679 set_op (matrix<T,W>& n, const Matrix& m) {
00680 typedef implementation<vector_linear,V> Vec;
00681 format<T> fm= CF(n);
00682 if (is_a_scalar (m)) {
00683 T a= get_sample (fm);
00684 set_project (a, m.scalar());
00685 n= matrix<T,W> (a, fm);
00686 } else {
00687 nat r= rows (m), c= cols (m);
00688 nat l= aligned_size<T,W> (r * c);
00689 T* t= mmx_formatted_new<T> (l, fm);
00690 Vec::template vec_unary<project_op> (t, tab (m), r * c);
00691 n= matrix<T,W> (t, r, c, fm); } }
00692 static inline Projected_matrix
00693 op (const Matrix& m) {
00694 Projected_matrix n (project (CF(m)));
00695 set_op (n, m);
00696 return n; }
00697 };
00698
00699 TMPL Reconstructed_matrix reconstruct (const Matrix& v) {
00700 return as<Reconstructed_matrix> (unary_map<reconstruct_op> (v)); }
00701
00702 template<typename C,typename V,typename W> bool
00703 is_reconstructible (const Matrix& v, matrix<Reconstruct_type(C),W>& w) {
00704 nat r= rows (v), c= cols (v);
00705 w= fill (promote (0, CF(w)), r, c);
00706 for (nat i= 0; i < r; i++)
00707 for (nat j= 0; i < c; j++)
00708 if (!is_reconstructible (v(i,j), w(i,j))) return false;
00709 return true; }
00710
00711
00712
00713
00714
00715 TMPLK Evaluated_matrix evaluate (const Matrix& v, const K& x) {
00716 return as<Evaluated_matrix> (binary_map_scalar<evaluate_op> (v, x)); }
00717
00718 template<typename C,typename V,typename K,typename W> bool
00719 is_evaluable (const Matrix& v, const K& a,
00720 matrix<Evaluate_type(C,K),W>& w) {
00721 nat r= rows (v), c= cols (v);
00722 w= matrix<Evaluate_type(C,K),W> (promote (0, CF(w)), r, c);
00723 for (nat i= 0; i < r; i++)
00724 for (nat j= 0; j < c; j++)
00725 if (!is_evaluable (v(i,j), a, w(i,j))) return false;
00726 return true; }
00727
00728
00729
00730
00731
00732 TMPL Truncated_matrix truncate (const Matrix& v, nat n) {
00733 return as<Truncated_matrix> (binary_map_scalar<truncate_op> (v, n)); }
00734
00735 TMPL Completed_matrix complete (const Matrix& v) {
00736 return as<Completed_matrix> (unary_map<complete_op> (v)); }
00737
00738
00739
00740
00741
00742 TMPL inline bool is_finite (const Matrix& m) {
00743 if (is_a_scalar (m)) return is_finite (m.scalar());
00744 return big<and_is_finite_op> (m); }
00745 TMPL inline bool is_nan (const Matrix& m) {
00746 if (is_a_scalar (m)) return is_nan (m.scalar());
00747 return big<or_is_nan_op> (m); }
00748 TMPL inline bool is_infinite (const Matrix& m) {
00749 if (is_a_scalar (m)) return is_infinite (m.scalar());
00750 return !is_nan (m) && big<or_is_infinite_op> (m); }
00751 TMPL inline bool is_fuzz (const Matrix& m) {
00752 if (is_a_scalar (m)) return is_fuzz (m.scalar());
00753 return !is_nan (m) && big<or_is_fuzz_op> (m); }
00754 TMPL inline bool is_reliable (const Matrix& m) {
00755 if (is_a_scalar (m)) return is_reliable (m.scalar());
00756 return is_reliable (C (0)); }
00757
00758 TMPL inline Matrix change_precision (const Matrix& m, xnat p) {
00759 return binary_map_scalar<change_precision_op> (m, p); }
00760 TMPL inline xnat precision (const Matrix& m) {
00761 return big<min_precision_op> (m); }
00762
00763 TMPL inline xint exponent (const Matrix& m) {
00764 return big<max_exponent_op> (m); }
00765 TMPL inline double magnitude (const Matrix& m) {
00766 return big<max_magnitude_op> (m); }
00767
00768
00769
00770
00771
00772 TMPL Abs_matrix abs (const Matrix& m) { return unary_map<abs_op> (m); }
00773 TMPL Real_matrix Re (const Matrix& m) { return unary_map<Re_op> (m); }
00774 TMPL Real_matrix Im (const Matrix& m) { return unary_map<Im_op> (m); }
00775 TMPL Matrix conj (const Matrix& m) { return unary_map<conj_op> (m); }
00776
00777 TMPL Center_matrix center (const Matrix& m) {
00778 return unary_map<center_op> (m); }
00779 TMPL Radius_matrix radius (const Matrix& m) {
00780 return unary_map<radius_op> (m); }
00781 TMPL Center_matrix lower (const Matrix& m) {
00782 return unary_map<lower_op> (m); }
00783 TMPL Center_matrix upper (const Matrix& m) {
00784 return unary_map<upper_op> (m); }
00785 TMPL inline Matrix sharpen (const Matrix& m) {
00786 return unary_map<sharpen_op> (m); }
00787 TMPLK inline Matrix blur (const Matrix& m, const K& x) {
00788 return binary_map_scalar<blur_op> (m, x); }
00789 template<typename C, typename V, typename C2, typename V2> inline Matrix
00790 blur (const Matrix& m, const matrix<C2,V2>& n) {
00791 return binary_map<blur_op> (m, n); }
00792 TMPL inline bool included (const Matrix& m, const Matrix& n) {
00793 return binary_test<included_op> (m, n); }
00794
00795
00796
00797
00798
00799 template<typename C> matrix<C>
00800 identity_matrix (const int& n, const format<C>& fm= format<C> ()) {
00801 return matrix<C> (promote (1, fm), n, n);
00802 }
00803
00804 template<typename C> matrix<C>
00805 fill_matrix (const C& x, const int& r, const int& c) {
00806 matrix<C> m (promote (0, x), (nat) r, (nat) c);
00807 for (nat i=0; i<((nat) r); i++)
00808 for (nat j=0; j<((nat) c); j++)
00809 m (i, j)= x;
00810 return m;
00811 }
00812
00813 template<typename C> matrix<C>
00814 tensor_matrix (const vector<C>& v, const vector<C>& w) {
00815 ASSERT (is_non_scalar (v), "non-scalar vector expected");
00816 ASSERT (is_non_scalar (w), "non-scalar vector expected");
00817 matrix<C> m (promote (0, CF(v)), N(v), N(w));
00818 for (nat i=0; i<N(v); i++)
00819 for (nat j=0; j<N(w); j++)
00820 m (i, j)= v[i] * w[j];
00821 return m;
00822 }
00823
00824 template<typename C> matrix<C>
00825 hilbert_matrix (const int& n, const format<C>& fm= format<C> ()) {
00826 matrix<C> m (promote (0, fm), n, n);
00827 for (nat i=0; i<(nat) n; i++)
00828 for (nat j=0; j<(nat) n; j++)
00829 m (i, j)= 1 / promote (1 + i + j, fm);
00830 return m;
00831 }
00832
00833 template<typename C> matrix<C>
00834 vandermonde (const vector<C>& v) {
00835 ASSERT (is_non_scalar (v), "non-scalar vector expected");
00836 matrix<C> m (promote (0, CF(v)), N(v), N(v));
00837 for (nat i=0; i<N(v); i++) {
00838 C p= promote (1, CF(v));
00839 for (nat j=0; j<N(v); j++, p *= v[i])
00840 m (i, j)= p;
00841 }
00842 return m;
00843 }
00844
00845 template<typename C> matrix<C>
00846 jordan_matrix (const C& c, const int& n) {
00847 matrix<C> m (c, n, n);
00848 for (nat i=0; i<((nat) (n-1)); i++)
00849 m (i, i+1)= promote (1, c);
00850 return m;
00851 }
00852
00853 template<typename C> matrix<C>
00854 toeplitz_matrix (const vector<C>& v) {
00855 ASSERT ((N(v)&1) == 1, "odd dimension expected");
00856 nat n= (N(v) >> 1) + 1;
00857 matrix<C> m (promote (0, CF(v)), n, n);
00858 for (nat i=0; i<n; i++)
00859 for (nat j=0; j<n; j++)
00860 m (i, j)= v [n-1 + i-j];
00861 return m;
00862 }
00863
00864 template<typename C> matrix<C>
00865 hankel_matrix (const vector<C>& v) {
00866 ASSERT ((N(v)&1) == 1, "odd dimension expected");
00867 nat n= (N(v) >> 1) + 1;
00868 matrix<C> m (promote (0, CF(v)), n, n);
00869 for (nat i=0; i<n; i++)
00870 for (nat j=0; j<n; j++)
00871 m (i, j)= v [i+j];
00872 return m;
00873 }
00874
00875
00876
00877
00878
00879 TMPL Matrix
00880 horizontal_join (const Matrix& m1, const Matrix& m2) {
00881 ASSERT (is_non_scalar (m1) || is_non_scalar (m2),
00882 "non-scalar matrix expected");
00883 if (!is_non_scalar (m1))
00884 return horizontal_join (Matrix (m1.scalar(), rows (m2), rows (m2)), m2);
00885 if (!is_non_scalar (m2))
00886 return horizontal_join (m1, Matrix (m2.scalar(), rows (m1), rows (m1)));
00887 ASSERT (rows (m1) == rows (m2), "unequal number of rows");
00888 Matrix r (promote (0, CF(m1)), rows (m1), cols (m1) + cols (m2));
00889 for (nat i=0; i<rows(m1); i++) {
00890 for (nat j=0; j<cols(m1); j++)
00891 r(i,j)= m1(i,j);
00892 for (nat j=0; j<cols(m2); j++)
00893 r(i,j+cols(m1))= m2(i,j);
00894 }
00895 return r;
00896 }
00897
00898 TMPL Matrix
00899 vertical_join (const Matrix& m1, const Matrix& m2) {
00900 ASSERT (is_non_scalar (m1) || is_non_scalar (m2),
00901 "non-scalar matrix expected");
00902 if (!is_non_scalar (m1))
00903 return vertical_join (Matrix (m1.scalar(), cols (m2), cols (m2)), m2);
00904 if (!is_non_scalar (m2))
00905 return vertical_join (m1, Matrix (m2.scalar(), cols (m1), cols (m1)));
00906 ASSERT (cols (m1) == cols (m2), "unequal number of columns");
00907 Matrix r (promote (0, CF(m1)), rows (m1) + rows (m2), cols (m1));
00908 for (nat j=0; j<cols(m1); j++) {
00909 for (nat i=0; i<rows(m1); i++)
00910 r(i,j)= m1(i,j);
00911 for (nat i=0; i<rows(m2); i++)
00912 r(i+rows(m1),j)= m2(i,j);
00913 }
00914 return r;
00915 }
00916
00917 TMPL Matrix
00918 transpose (const Matrix & m) {
00919 typedef implementation<matrix_linear,V> Mat;
00920 if (is_a_scalar (m)) return m;
00921 nat nrows= rows (m), ncols= cols (m);
00922 nat l= aligned_size<C,V> (nrows * ncols);
00923 C* r= mmx_formatted_new<C> (l, CF(m));
00924 Mat::transpose (r, tab (m), nrows, ncols);
00925 return Matrix (r, ncols, nrows, CF(m));
00926 }
00927
00928 TMPL inline Matrix
00929 range (const Matrix& m, nat r1, nat c1, nat r2, nat c2) {
00930 typedef implementation<matrix_linear,V> Mat;
00931 if (is_a_scalar (m)) return Matrix (m.scalar(), r2-r1, c2-c1);
00932 nat nrows= rows (m), ncols= cols (m);
00933 nat l= aligned_size<C,V> ((r2-r1) * (c2-c1));
00934 C* r= mmx_formatted_new<C> (l, CF(m));
00935 Mat::get_range (r, tab (m), r1, c1, r2, c2, nrows, ncols);
00936 return Matrix (r, r2-r1, c2-c1, CF(m));
00937 }
00938
00939 TMPL inline Matrix
00940 delete_row (const Matrix& m, nat r) {
00941 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
00942 ASSERT (r < rows (m), "index out of range");
00943 Matrix d (promote (0, CF(m)), rows (m) - 1, cols (m));
00944 for (nat i= 0; i < r; i++)
00945 for (nat j= 0; j < cols (m); j++) d(i,j)= m(i,j);
00946 for (nat i= r+1; i < rows (m); i++)
00947 for (nat j= 0; j < cols (m); j++) d(i-1,j)= m(i,j);
00948 return d;
00949 }
00950
00951 TMPL Matrix
00952 delete_col (const Matrix& m, nat c) {
00953 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
00954 ASSERT (c < cols (m), "index out of range");
00955 Matrix d (promote (0, CF(m)), rows (m), cols (m) - 1);
00956 for (nat j= 0; j < c; j++)
00957 for (nat i= 0; i < rows (m); i++) d(i,j)= m(i,j);
00958 for (nat j= c+1; j < cols (m); j++)
00959 for (nat i= 0; i < rows (m); i++) d(i,j-1)= m(i,j);
00960 return d;
00961 }
00962
00963 TMPL Matrix
00964 operator * (const Matrix& m, const Matrix& n) {
00965 typedef implementation<matrix_multiply,V> Mat;
00966 if (is_a_scalar (m) || is_a_scalar (n)) {
00967 if (is_non_scalar (m)) return m * n.scalar();
00968 if (is_non_scalar (n)) return m.scalar() * n;
00969 return Matrix (m.scalar() * n.scalar());
00970 }
00971 nat mrows= rows (m), mcols= cols (m), nrows= rows(n), ncols= cols(n);
00972 ASSERT (nrows == mcols, "numbers of rows and columns don't match");
00973 nat l= aligned_size<C,V> (mrows * ncols);
00974 C* r= mmx_formatted_new<C> (l, CF(m));
00975 Mat::mul (r, tab (m), tab (n), mrows, mcols, ncols);
00976 return Matrix (r, mrows, ncols, CF(m));
00977 }
00978
00979 template<typename C, typename V, typename W> vector<C,W>
00980 operator * (const matrix<C,V>& m, const vector<C,W>& v) {
00981 typedef implementation<vector_linear,W> Vec;
00982 typedef implementation<matrix_linear,V> Mat;
00983 if (is_a_scalar (m)) return m.scalar() * v;
00984 ASSERT (is_non_scalar (v), "non-scalar vector expected");
00985 nat rr= rows (m), cc= cols (m);
00986 ASSERT (cc == N(v), "sizes don't match");
00987 nat l= aligned_size<C,W> (rr);
00988 C* a= mmx_formatted_new<C> (l, CF(m));
00989 for (nat i=0; i<rr; i++)
00990 a[i]= Vec::template vec_binary_big_stride<mul_add_op>
00991 (tab (m) + Mat::index (i, 0, rr, cc), Mat::index (0, 1, rr, cc),
00992 seg (v), 1, cc);
00993 return vector<C,W> (a, rr, l, CF(m));
00994 }
00995
00996 template<typename C, typename V, typename W> vector<C,W>
00997 operator * (const vector<C,W>& v, const matrix<C,V>& m) {
00998 typedef implementation<vector_linear,W> Vec;
00999 typedef implementation<matrix_linear,V> Mat;
01000 if (is_a_scalar (m)) return v * m.scalar();
01001 ASSERT (is_non_scalar (v), "non-scalar vector expected");
01002 nat rr= rows (m), cc= cols (m);
01003 ASSERT (rr == N(v), "sizes don't match");
01004 nat l= aligned_size<C,W> (cc);
01005 C* a= mmx_formatted_new<C> (l, CF(m));
01006 for (nat i=0; i<cc; i++)
01007 a[i]= Vec::template vec_binary_big_stride<mul_add_op>
01008 (seg (v), 1,
01009 tab (m) + Mat::index (0, i, rr, cc), Mat::index (1, 0, rr, cc), rr);
01010 return vector<C,W> (a, cc, l, CF(m));
01011 }
01012
01013
01014
01015
01016
01017 TMPL vector<C>
01018 row (const Matrix& m, nat i) {
01019 nat n= cols (m);
01020 nat l= default_aligned_size<C> (n);
01021 C* a= mmx_formatted_new<C> (l, CF(m));
01022 for (nat j=0; j<n; j++) a[j]= m(i,j);
01023 return vector<C> (a, n, l, CF(m));
01024 }
01025
01026 TMPL vector<C>
01027 column (const Matrix& m, nat j) {
01028 nat n= rows (m);
01029 nat l= default_aligned_size<C> (n);
01030 C* a= mmx_formatted_new<C> (l, CF(m));
01031 for (nat i=0; i<n; i++) a[i]= m(i,j);
01032 return vector<C> (a, n, l, CF(m));
01033 }
01034
01035 template<typename C> vector<vector<C> >
01036 row_vectors (const matrix<C>& m) {
01037 vector<vector<C> > v=
01038 fill<vector<C> > (rows (m), format<vector<C> > (CF(m)));
01039 for (nat i=0; i<rows(m); i++)
01040 v[i]= row (m, i);
01041 return v;
01042 }
01043
01044 template<typename C> matrix<C>
01045 row_matrix (const vector<vector<C> >& v) {
01046 if (N(v) == 0) return matrix<C> (promote (0, get_format1 (CF(v))));
01047 matrix<C> r (promote (0, get_format1 (CF(v))), N(v), N(v[0]));
01048 for (nat i=0; N(v)>i; i++) {
01049 ASSERT (N(v[i]) == N(v[0]), "unequal row lengths");
01050 for (nat j=0; j<N(v[0]); j++)
01051 r (i, j)= v[i][j];
01052 }
01053 return r;
01054 }
01055
01056 TMPL void
01057 swap_row (Matrix& m, nat i, nat j) {
01058 typedef implementation<matrix_linear,V> Mat;
01059 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01060 nat mrows= rows (m), mcols= cols (m);
01061 ASSERT (i < mrows && j < mrows, "out of range");
01062 Mat::row_swap (tab (m), i, j, mrows, mcols);
01063 }
01064
01065 TMPL void
01066 swap_col (Matrix& m, nat i, nat j) {
01067 typedef implementation<matrix_linear,V> Mat;
01068 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01069 nat mrows= rows (m), mcols= cols (m);
01070 ASSERT (i < mcols && j < mcols, "out of range");
01071 Mat::col_swap (tab (m), i, j, mrows, mcols);
01072 }
01073
01074 TMPL void
01075 reverse_cols (Matrix& m) {
01076 typedef implementation<matrix_linear,V> Mat;
01077 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01078 Mat::col_reverse (tab (m), rows (m), cols (m));
01079 }
01080
01081 #define INLINE_SCAL_LINE_OP(name1,name2) \
01082 TMPL void name1 (Matrix& m, C c, nat i) { \
01083 typedef implementation<matrix_linear,V> Mat; \
01084 ASSERT (is_non_scalar (m), "non-scalar matrix expected"); \
01085 nat mrows=rows(m), mcols=cols(m); \
01086 Mat::name2 (tab(m), c, i, mrows, mcols); \
01087 }
01088
01089 #define INLINE_LINES_OP(name1, name2) \
01090 TMPL void name1 (Matrix& m, nat i, C ci, nat j, C cj) { \
01091 typedef implementation<matrix_linear,V> Mat; \
01092 ASSERT (is_non_scalar (m), "non-scalar matrix expected"); \
01093 nat mrows=rows(m), mcols=cols(m); \
01094 Mat::name2 (tab(m), i, ci, j, cj, mrows, mcols); \
01095 }
01096
01097 INLINE_SCAL_LINE_OP(row_mul, row_mul);
01098 INLINE_SCAL_LINE_OP(row_div, row_div);
01099 INLINE_SCAL_LINE_OP(col_mul, col_mul);
01100 INLINE_SCAL_LINE_OP(col_div, col_div);
01101 INLINE_LINES_OP(rows_linsub,row_combine_sub);
01102 INLINE_LINES_OP(cols_linsub,col_combine_sub);
01103
01104
01105
01106
01107
01108 TMPL matrix<C>
01109 as_matrix (const permutation& p, const Format& fm= Format ()) {
01110 nat n= N(p);
01111 matrix<C> m (promote (0, fm), n, n);
01112 for (nat i=0; i<n; i++)
01113 for (nat j=0; j<n; j++)
01114 m (i, j) = promote (i == p (j) ? 1 : 0, fm);
01115 return m;
01116 }
01117
01118 TMPL void
01119 permute_columns (Matrix& m, const permutation& p) {
01120
01121 typedef implementation<matrix_permute,V> Mat;
01122 if (is_a_scalar (m)) return;
01123 nat rs= rows (m), cs= cols (m);
01124 Mat::col_permute (tab (m), seg (*p), rs, cs);
01125 }
01126
01127 TMPL void
01128 permute_rows (Matrix& m, const permutation& p) {
01129
01130 typedef implementation<matrix_permute,V> Mat;
01131 if (is_a_scalar (m)) return;
01132 nat rs= rows (m), cs= cols (m);
01133 Mat::row_permute (tab (m), seg (*p), rs, cs);
01134 }
01135
01136 TMPL Matrix
01137 operator * (const Matrix& m, const permutation& p) {
01138
01139 typedef implementation<matrix_permute,V> Mat;
01140 if (is_a_scalar (m)) return m;
01141 nat rs= rows (m), cs= cols (m);
01142 Matrix dest (promote (0, CF(m)), rs, cs);
01143 Mat::col_permute (tab (dest), tab (m), seg (*p), rs, cs);
01144 return dest;
01145 }
01146
01147 TMPL Matrix
01148 operator * (const permutation& p, const Matrix& m) {
01149
01150 typedef implementation<matrix_permute,V> Mat;
01151 if (is_a_scalar (m)) return m;
01152 nat rs= rows (m), cs= cols (m);
01153 Matrix dest (promote (0, CF(m)), rs, cs);
01154 Mat::row_permute (tab (dest), tab (m), seg (*p), rs, cs);
01155 return dest;
01156 }
01157
01158
01159
01160
01161
01162 TMPL Matrix
01163 column_echelon (const Matrix& m, bool reduced= false) {
01164 typedef implementation<matrix_echelon,V> Mat;
01165 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01166 Matrix c= copy (m);
01167 C* k= NULL;
01168 Mat::col_echelon (tab(c), k, rows(m), cols(m), reduced);
01169 return c;
01170 }
01171
01172 TMPL Matrix
01173 column_reduced_echelon (const Matrix& m) {
01174 return column_echelon (m, true);
01175 }
01176
01177 TMPL Matrix
01178 column_reduced_echelon (const Matrix& m, permutation& permut) {
01179 typedef implementation<matrix_echelon,V> Mat;
01180 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01181 Matrix c= copy (m);
01182 C* k= NULL;
01183 nat* buf= mmx_new<nat> (default_aligned_size<nat> (cols(m)));
01184 Mat::col_echelon (tab(c), k, rows(m), cols(m), true, buf);
01185 permut= permutation (vector<nat> (buf, cols(m), format<nat> ()));
01186 return c;
01187 }
01188
01189 TMPL Matrix
01190 column_echelon (const Matrix& m, Matrix& k, bool reduced= false) {
01191 typedef implementation<matrix_echelon,V> Mat;
01192 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01193 Matrix c= copy (m);
01194 k= Matrix (promote (0, CF(m)), cols(m), cols(m));
01195 Mat::col_echelon (tab(c), tab(k), rows(m), cols(m), false);
01196 return c;
01197 }
01198
01199 TMPL Matrix
01200 column_reduced_echelon (const Matrix& m, Matrix& k) {
01201 return column_echelon (m, k, true);
01202 }
01203
01204 TMPL inline Matrix
01205 row_echelon (const Matrix& m, bool reduced= false) {
01206 return transpose (column_echelon (transpose (m), reduced));
01207 }
01208
01209 TMPL Matrix
01210 row_reduced_echelon (const Matrix& m) {
01211 return row_echelon (m, true);
01212 }
01213
01214 TMPL inline Matrix
01215 row_echelon (const Matrix& m, Matrix& k, bool reduced= false) {
01216 Matrix c= column_echelon (transpose (m), k, reduced);
01217 k= transpose (k);
01218 return transpose (c);
01219 }
01220
01221 TMPL Matrix
01222 row_reduced_echelon (const Matrix& m, Matrix& k) {
01223 return row_echelon (m, k, true);
01224 }
01225
01226
01227
01228
01229
01230 TMPL inline C
01231 det (const Matrix& m) {
01232 typedef implementation<matrix_determinant,V> Mat;
01233 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01234 ASSERT (cols(m) == rows(m), "square matrix expected");
01235 C r= promote (1, CF(m));
01236 Mat::det (r, tab(m), rows(m));
01237 return r;
01238 }
01239
01240 TMPL inline C
01241 first_minor (const Matrix& m, nat i, nat j) {
01242 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01243 ASSERT (cols(m) == rows(m), "square matrix expected");
01244 ASSERT (i < rows(m), "index out of range");
01245 ASSERT (j < cols(m), "index out of range");
01246 return det (delete_col (delete_row (m, i), j));
01247 }
01248
01249 TMPL inline C
01250 cofactor (const Matrix& m, nat i, nat j) {
01251 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01252 ASSERT (cols(m) == rows(m), "square matrix expected");
01253 ASSERT (i < rows(m), "index out of range");
01254 ASSERT (j < cols(m), "index out of range");
01255 C c= first_minor (m, i, j);
01256 return ((i + j) & 1) ? -c : c;
01257 }
01258
01259 TMPL Matrix
01260 invert (const Matrix& m) {
01261 typedef implementation<matrix_invert,V> Mat;
01262 if (is_a_scalar (m)) return Matrix (invert (m.scalar()));
01263 ASSERT (cols(m) == rows(m), "square matrix expected");
01264 Matrix a (promote (1, CF(m)), rows(m), cols(m));
01265 Mat::invert (tab(a), tab(m), rows(m));
01266 return a;
01267 }
01268
01269 TMPL Matrix operator / (const C& c, const Matrix& m) {
01270 return c * invert (m); }
01271
01272
01273
01274
01275
01276 TMPL Matrix
01277 kernel (const Matrix& m) {
01278 typedef implementation<matrix_kernel,V> Mat;
01279 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01280 Matrix k (promote (0, CF(m)), cols(m), cols(m));
01281 nat dim= Mat::kernel (tab(k), tab(m), rows(m), cols(m));
01282 return range (k, 0, 0, cols(m), dim);
01283 }
01284
01285 TMPL Matrix
01286 image (const Matrix& m) {
01287 typedef implementation<matrix_image,V> Mat;
01288 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01289 Matrix k (promote (0, CF(m)), rows(m), rows(m));
01290 nat dim= Mat::image (tab(k), tab(m), rows(m), cols(m));
01291 return range (k, 0, 0, rows(m), dim);
01292 }
01293
01294 TMPL nat
01295 rank (const Matrix& m) {
01296 typedef implementation<matrix_image,V> Mat;
01297 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01298 return Mat::rank (tab(m), rows(m), cols(m));
01299 }
01300
01301 TMPL Matrix
01302 krylov (const Matrix& m, const Matrix& v) {
01303 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01304 ASSERT (is_square_matrix (m), "square matrix expected");
01305 Matrix r= image (v);
01306 Matrix p= m;
01307 while (true) {
01308 nat rk= rows (r);
01309 r= image (vertical_join (r, r*p));
01310 if (rows (r) <= rk) return r;
01311 p= p*p;
01312 }
01313 }
01314
01315
01316
01317
01318
01319 TMPL inline Matrix
01320 row_orthogonalization (const Matrix& m) {
01321 typedef implementation<matrix_orthogonalization,V> Mat;
01322 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01323 Matrix c= copy (m);
01324 vector<C> n (promote (0, CF(m)), rows(m));
01325 Mat::row_orthogonalize (tab(c), rows(m), cols(m), seg(n));
01326 return c;
01327 }
01328
01329 TMPL inline Matrix
01330 column_orthogonalization (const Matrix& m) {
01331 typedef implementation<matrix_orthogonalization,V> Mat;
01332 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01333 Matrix c= copy (m);
01334 vector<C> n (promote (0, CF(m)), cols(m));
01335 Mat::col_orthogonalize (tab(c), rows(m), cols(m), seg(n));
01336 return c;
01337 }
01338
01339 TMPL inline Matrix
01340 row_orthogonalization (const Matrix& m, Matrix& l) {
01341 typedef implementation<matrix_orthogonalization,V> Mat;
01342 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01343 Matrix c= copy (m);
01344 vector<C> n (promote (0, CF(m)), rows(m));
01345 l= Matrix (promote (0, CF(m)), rows(m), rows(m));
01346 Mat::row_orthogonalize (tab(c), rows(m), cols(m), tab(l), seg(n));
01347 return c;
01348 }
01349
01350 TMPL inline Matrix
01351 column_orthogonalization (const Matrix& m, Matrix& l) {
01352 typedef implementation<matrix_orthogonalization,V> Mat;
01353 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01354 Matrix c= copy (m);
01355 vector<C> n (promote (0, CF(m)), cols(m));
01356 l= Matrix (promote (0, CF(m)), cols(m), cols(m));
01357 Mat::col_orthogonalize (tab(c), rows(m), cols(m), tab(l), seg(n));
01358 return c;
01359 }
01360
01361 TMPL inline Matrix
01362 row_orthonormalization (const Matrix& m) {
01363 typedef implementation<matrix_orthogonalization,V> Mat;
01364 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01365 Matrix c= copy (m);
01366 Mat::row_orthonormalize (tab(c), rows(m), cols(m));
01367 return c;
01368 }
01369
01370 TMPL inline Matrix
01371 column_orthonormalization (const Matrix& m) {
01372 typedef implementation<matrix_orthogonalization,V> Mat;
01373 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01374 Matrix c= copy (m);
01375 Mat::col_orthonormalize (tab(c), rows(m), cols(m));
01376 return c;
01377 }
01378
01379 TMPL inline Matrix
01380 row_orthonormalization (const Matrix& m, Matrix& l) {
01381 typedef implementation<matrix_orthogonalization,V> Mat;
01382 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01383 Matrix c= copy (m);
01384 l= Matrix (promote (0, CF(m)), rows(m), rows(m));
01385 Mat::row_orthonormalize (tab(c), rows(m), cols(m), tab(l));
01386 return c;
01387 }
01388
01389 TMPL inline Matrix
01390 column_orthonormalization (const Matrix& m, Matrix& l) {
01391 typedef implementation<matrix_orthogonalization,V> Mat;
01392 ASSERT (is_non_scalar (m), "non-scalar matrix expected");
01393 Matrix c= copy (m);
01394 l= Matrix (promote (0, CF(m)), cols(m), cols(m));
01395 Mat::col_orthonormalize (tab(c), rows(m), cols(m), tab(l));
01396 return c;
01397 }
01398
01399 #undef TMPL_DEF
01400 #undef TMPL
01401 #undef TMPLK
01402 #undef Format
01403 #undef Matrix
01404 #undef Matrix_rep
01405 #undef Abs_matrix
01406 #undef Real_matrix
01407 #undef Center_matrix
01408 #undef Radius_matrix
01409 #undef Lifted_matrix
01410 #undef Projected_matrix
01411 #undef Reconstructed_matrix
01412 #undef Truncated_matrix
01413 #undef Completed_matrix
01414 #undef Evaluated_matrix
01415 }
01416 #endif // __MMX_MATRIX_HPP
