1 // -*- c++ -*- (enables emacs c++ mode) 2 //=========================================================================== 3 // 4 // Copyright (C) 2002-2008 Yves Renard 5 // 6 // This file is a part of GETFEM++ 7 // 8 // Getfem++ is free software; you can redistribute it and/or modify it 9 // under the terms of the GNU Lesser General Public License as published 10 // by the Free Software Foundation; either version 2.1 of the License, or 11 // (at your option) any later version. 12 // This program is distributed in the hope that it will be useful, but 13 // WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 14 // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 15 // License for more details. 16 // You should have received a copy of the GNU Lesser General Public License 17 // along with this program; if not, write to the Free Software Foundation, 18 // Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. 19 // 20 // As a special exception, you may use this file as part of a free software 21 // library without restriction. Specifically, if other files instantiate 22 // templates or use macros or inline functions from this file, or you compile 23 // this file and link it with other files to produce an executable, this 24 // file does not by itself cause the resulting executable to be covered by 25 // the GNU General Public License. This exception does not however 26 // invalidate any other reasons why the executable file might be covered by 27 // the GNU General Public License. 28 // 29 //=========================================================================== 30 31 /**@file gmm_sub_matrix.h 32 @author Yves Renard <Yves.Renard@insa-lyon.fr> 33 @date October 13, 2002. 34 @brief Generic sub-matrices. 35 */ 36 37 #ifndef GMM_SUB_MATRIX_H__ 38 #define GMM_SUB_MATRIX_H__ 39 40 #include "gmm_sub_vector.h" 41 42 namespace gmm { 43 44 /* ********************************************************************* */ 45 /* sub row matrices type */ 46 /* ********************************************************************* */ 47 48 template <typename PT, typename SUBI1, typename SUBI2> 49 struct gen_sub_row_matrix { 50 typedef gen_sub_row_matrix<PT, SUBI1, SUBI2> this_type; 51 typedef typename std::iterator_traits<PT>::value_type M; 52 typedef M * CPT; 53 typedef typename std::iterator_traits<PT>::reference ref_M; 54 typedef typename select_ref<typename linalg_traits<M> 55 ::const_row_iterator, typename linalg_traits<M>::row_iterator, 56 PT>::ref_type iterator; 57 typedef typename linalg_traits<this_type>::reference reference; 58 typedef typename linalg_traits<this_type>::porigin_type porigin_type; 59 60 SUBI1 si1; 61 SUBI2 si2; 62 iterator begin_; 63 porigin_type origin; 64 operatorgen_sub_row_matrix65 reference operator()(size_type i, size_type j) const 66 { return linalg_traits<M>::access(begin_ + si1.index(i), si2.index(j)); } 67 nrowsgen_sub_row_matrix68 size_type nrows(void) const { return si1.size(); } ncolsgen_sub_row_matrix69 size_type ncols(void) const { return si2.size(); } 70 gen_sub_row_matrixgen_sub_row_matrix71 gen_sub_row_matrix(ref_M m, const SUBI1 &s1, const SUBI2 &s2) 72 : si1(s1), si2(s2), begin_(mat_row_begin(m)), 73 origin(linalg_origin(m)) {} gen_sub_row_matrixgen_sub_row_matrix74 gen_sub_row_matrix() {} gen_sub_row_matrixgen_sub_row_matrix75 gen_sub_row_matrix(const gen_sub_row_matrix<CPT, SUBI1, SUBI2> &cr) : 76 si1(cr.si1), si2(cr.si2), begin_(cr.begin_),origin(cr.origin) {} 77 }; 78 79 template <typename PT, typename SUBI1, typename SUBI2> 80 struct gen_sub_row_matrix_iterator { 81 typedef gen_sub_row_matrix<PT, SUBI1, SUBI2> this_type; 82 typedef typename modifiable_pointer<PT>::pointer MPT; 83 typedef typename std::iterator_traits<PT>::value_type M; 84 typedef typename select_ref<typename linalg_traits<M> 85 ::const_row_iterator, typename linalg_traits<M>::row_iterator, 86 PT>::ref_type ITER; 87 typedef ITER value_type; 88 typedef ITER *pointer; 89 typedef ITER &reference; 90 typedef ptrdiff_t difference_type; 91 typedef size_t size_type; 92 typedef std::random_access_iterator_tag iterator_category; 93 typedef gen_sub_row_matrix_iterator<PT, SUBI1, SUBI2> iterator; 94 95 ITER it; 96 SUBI1 si1; 97 SUBI2 si2; 98 size_type ii; 99 100 iterator operator ++(int) { iterator tmp = *this; ii++; return tmp; } 101 iterator operator --(int) { iterator tmp = *this; ii--; return tmp; } 102 iterator &operator ++() { ii++; return *this; } 103 iterator &operator --() { ii--; return *this; } 104 iterator &operator +=(difference_type i) { ii += i; return *this; } 105 iterator &operator -=(difference_type i) { ii -= i; return *this; } 106 iterator operator +(difference_type i) const 107 { iterator itt = *this; return (itt += i); } 108 iterator operator -(difference_type i) const 109 { iterator itt = *this; return (itt -= i); } 110 difference_type operator -(const iterator &i) const { return ii - i.ii; } 111 112 ITER operator *() const { return it + si1.index(ii); } 113 ITER operator [](int i) { return it + si1.index(ii+i); } 114 115 bool operator ==(const iterator &i) const { return (ii == i.ii); } 116 bool operator !=(const iterator &i) const { return !(i == *this); } 117 bool operator < (const iterator &i) const { return (ii < i.ii); } 118 gen_sub_row_matrix_iteratorgen_sub_row_matrix_iterator119 gen_sub_row_matrix_iterator(void) {} gen_sub_row_matrix_iteratorgen_sub_row_matrix_iterator120 gen_sub_row_matrix_iterator(const 121 gen_sub_row_matrix_iterator<MPT, SUBI1, SUBI2> &itm) 122 : it(itm.it), si1(itm.si1), si2(itm.si2), ii(itm.ii) {} gen_sub_row_matrix_iteratorgen_sub_row_matrix_iterator123 gen_sub_row_matrix_iterator(const ITER &iter, const SUBI1 &s1, 124 const SUBI2 &s2, size_type i) 125 : it(iter), si1(s1), si2(s2), ii(i) { } 126 127 }; 128 129 template <typename PT, typename SUBI1, typename SUBI2> 130 struct linalg_traits<gen_sub_row_matrix<PT, SUBI1, SUBI2> > { 131 typedef gen_sub_row_matrix<PT, SUBI1, SUBI2> this_type; 132 typedef typename std::iterator_traits<PT>::value_type M; 133 typedef typename which_reference<PT>::is_reference is_reference; 134 typedef abstract_matrix linalg_type; 135 typedef typename linalg_traits<M>::origin_type origin_type; 136 typedef typename select_ref<const origin_type *, origin_type *, 137 PT>::ref_type porigin_type; 138 typedef typename linalg_traits<M>::value_type value_type; 139 typedef typename select_ref<value_type, 140 typename linalg_traits<M>::reference, PT>::ref_type reference; 141 typedef abstract_null_type sub_col_type; 142 typedef abstract_null_type col_iterator; 143 typedef abstract_null_type const_sub_col_type; 144 typedef abstract_null_type const_col_iterator; 145 typedef typename sub_vector_type<const typename 146 linalg_traits<M>::const_sub_row_type *, SUBI2>::vector_type 147 const_sub_row_type; 148 typedef typename select_ref<abstract_null_type, 149 typename sub_vector_type<typename linalg_traits<M>::sub_row_type *, 150 SUBI2>::vector_type, PT>::ref_type sub_row_type; 151 typedef gen_sub_row_matrix_iterator<typename const_pointer<PT>::pointer, 152 SUBI1, SUBI2> const_row_iterator; 153 typedef typename select_ref<abstract_null_type, 154 gen_sub_row_matrix_iterator<PT, SUBI1, SUBI2>, PT>::ref_type 155 row_iterator; 156 typedef typename linalg_traits<const_sub_row_type>::storage_type 157 storage_type; 158 typedef row_major sub_orientation; 159 typedef linalg_true index_sorted; 160 static size_type nrows(const this_type &m) { return m.nrows(); } 161 static size_type ncols(const this_type &m) { return m.ncols(); } 162 static const_sub_row_type row(const const_row_iterator &it) 163 { return const_sub_row_type(linalg_traits<M>::row(*it), it.si2); } 164 static sub_row_type row(const row_iterator &it) 165 { return sub_row_type(linalg_traits<M>::row(*it), it.si2); } 166 static const_row_iterator row_begin(const this_type &m) 167 { return const_row_iterator(m.begin_, m.si1, m.si2, 0); } 168 static row_iterator row_begin(this_type &m) 169 { return row_iterator(m.begin_, m.si1, m.si2, 0); } 170 static const_row_iterator row_end(const this_type &m) 171 { return const_row_iterator(m.begin_, m.si1, m.si2, m.nrows()); } 172 static row_iterator row_end(this_type &m) 173 { return row_iterator(m.begin_, m.si1, m.si2, m.nrows()); } 174 static origin_type* origin(this_type &v) { return v.origin; } 175 static const origin_type* origin(const this_type &v) { return v.origin; } 176 static void do_clear(this_type &m) { 177 row_iterator it = mat_row_begin(m), ite = mat_row_end(m); 178 for (; it != ite; ++it) clear(row(it)); 179 } 180 static value_type access(const const_row_iterator &itrow, size_type i) 181 { return linalg_traits<M>::access(*itrow, itrow.si2.index(i)); } 182 static reference access(const row_iterator &itrow, size_type i) 183 { return linalg_traits<M>::access(*itrow, itrow.si2.index(i)); } 184 }; 185 186 template <typename PT, typename SUBI1, typename SUBI2> 187 std::ostream &operator <<(std::ostream &o, 188 const gen_sub_row_matrix<PT, SUBI1, SUBI2>& m) 189 { gmm::write(o,m); return o; } 190 191 192 /* ********************************************************************* */ 193 /* sub column matrices type */ 194 /* ********************************************************************* */ 195 196 template <typename PT, typename SUBI1, typename SUBI2> 197 struct gen_sub_col_matrix { 198 typedef gen_sub_col_matrix<PT, SUBI1, SUBI2> this_type; 199 typedef typename std::iterator_traits<PT>::value_type M; 200 typedef M * CPT; 201 typedef typename std::iterator_traits<PT>::reference ref_M; 202 typedef typename select_ref<typename linalg_traits<M> 203 ::const_col_iterator, typename linalg_traits<M>::col_iterator, 204 PT>::ref_type iterator; 205 typedef typename linalg_traits<this_type>::reference reference; 206 typedef typename linalg_traits<this_type>::porigin_type porigin_type; 207 208 SUBI1 si1; 209 SUBI2 si2; 210 iterator begin_; 211 porigin_type origin; 212 213 reference operator()(size_type i, size_type j) const 214 { return linalg_traits<M>::access(begin_ + si2.index(j), si1.index(i)); } 215 216 size_type nrows(void) const { return si1.size(); } 217 size_type ncols(void) const { return si2.size(); } 218 219 gen_sub_col_matrix(ref_M m, const SUBI1 &s1, const SUBI2 &s2) 220 : si1(s1), si2(s2), begin_(mat_col_begin(m)), 221 origin(linalg_origin(m)) {} 222 gen_sub_col_matrix() {} 223 gen_sub_col_matrix(const gen_sub_col_matrix<CPT, SUBI1, SUBI2> &cr) : 224 si1(cr.si1), si2(cr.si2), begin_(cr.begin_),origin(cr.origin) {} 225 }; 226 227 template <typename PT, typename SUBI1, typename SUBI2> 228 struct gen_sub_col_matrix_iterator { 229 typedef gen_sub_col_matrix<PT, SUBI1, SUBI2> this_type; 230 typedef typename modifiable_pointer<PT>::pointer MPT; 231 typedef typename std::iterator_traits<PT>::value_type M; 232 typedef typename select_ref<typename linalg_traits<M>::const_col_iterator, 233 typename linalg_traits<M>::col_iterator, 234 PT>::ref_type ITER; 235 typedef ITER value_type; 236 typedef ITER *pointer; 237 typedef ITER &reference; 238 typedef ptrdiff_t difference_type; 239 typedef size_t size_type; 240 typedef std::random_access_iterator_tag iterator_category; 241 typedef gen_sub_col_matrix_iterator<PT, SUBI1, SUBI2> iterator; 242 243 ITER it; 244 SUBI1 si1; 245 SUBI2 si2; 246 size_type ii; 247 248 iterator operator ++(int) { iterator tmp = *this; ii++; return tmp; } 249 iterator operator --(int) { iterator tmp = *this; ii--; return tmp; } 250 iterator &operator ++() { ii++; return *this; } 251 iterator &operator --() { ii--; return *this; } 252 iterator &operator +=(difference_type i) { ii += i; return *this; } 253 iterator &operator -=(difference_type i) { ii -= i; return *this; } 254 iterator operator +(difference_type i) const 255 { iterator itt = *this; return (itt += i); } 256 iterator operator -(difference_type i) const 257 { iterator itt = *this; return (itt -= i); } 258 difference_type operator -(const iterator &i) const { return ii - i.ii; } 259 260 ITER operator *() const { return it + si2.index(ii); } 261 ITER operator [](int i) { return it + si2.index(ii+i); } 262 263 bool operator ==(const iterator &i) const { return (ii == i.ii); } 264 bool operator !=(const iterator &i) const { return !(i == *this); } 265 bool operator < (const iterator &i) const { return (ii < i.ii); } 266 267 gen_sub_col_matrix_iterator(void) {} 268 gen_sub_col_matrix_iterator(const 269 gen_sub_col_matrix_iterator<MPT, SUBI1, SUBI2> &itm) 270 : it(itm.it), si1(itm.si1), si2(itm.si2), ii(itm.ii) {} 271 gen_sub_col_matrix_iterator(const ITER &iter, const SUBI1 &s1, 272 const SUBI2 &s2, size_type i) 273 : it(iter), si1(s1), si2(s2), ii(i) { } 274 }; 275 276 template <typename PT, typename SUBI1, typename SUBI2> 277 struct linalg_traits<gen_sub_col_matrix<PT, SUBI1, SUBI2> > { 278 typedef gen_sub_col_matrix<PT, SUBI1, SUBI2> this_type; 279 typedef typename std::iterator_traits<PT>::value_type M; 280 typedef typename linalg_traits<M>::origin_type origin_type; 281 typedef typename select_ref<const origin_type *, origin_type *, 282 PT>::ref_type porigin_type; 283 typedef typename which_reference<PT>::is_reference is_reference; 284 typedef abstract_matrix linalg_type; 285 typedef typename linalg_traits<M>::value_type value_type; 286 typedef typename select_ref<value_type, 287 typename linalg_traits<M>::reference, PT>::ref_type reference; 288 typedef abstract_null_type sub_row_type; 289 typedef abstract_null_type row_iterator; 290 typedef abstract_null_type const_sub_row_type; 291 typedef abstract_null_type const_row_iterator; 292 typedef typename sub_vector_type<const typename 293 linalg_traits<M>::const_sub_col_type *, SUBI1>::vector_type 294 const_sub_col_type; 295 typedef typename select_ref<abstract_null_type, 296 typename sub_vector_type<typename linalg_traits<M>::sub_col_type *, 297 SUBI1>::vector_type, PT>::ref_type sub_col_type; 298 typedef gen_sub_col_matrix_iterator<typename const_pointer<PT>::pointer, 299 SUBI1, SUBI2> const_col_iterator; 300 typedef typename select_ref<abstract_null_type, 301 gen_sub_col_matrix_iterator<PT, SUBI1, SUBI2>, PT>::ref_type 302 col_iterator; 303 typedef col_major sub_orientation; 304 typedef linalg_true index_sorted; 305 typedef typename linalg_traits<const_sub_col_type>::storage_type 306 storage_type; 307 static size_type nrows(const this_type &m) { return m.nrows(); } 308 static size_type ncols(const this_type &m) { return m.ncols(); } 309 static const_sub_col_type col(const const_col_iterator &it) 310 { return const_sub_col_type(linalg_traits<M>::col(*it), it.si1); } 311 static sub_col_type col(const col_iterator &it) 312 { return sub_col_type(linalg_traits<M>::col(*it), it.si1); } 313 static const_col_iterator col_begin(const this_type &m) 314 { return const_col_iterator(m.begin_, m.si1, m.si2, 0); } 315 static col_iterator col_begin(this_type &m) 316 { return col_iterator(m.begin_, m.si1, m.si2, 0); } 317 static const_col_iterator col_end(const this_type &m) 318 { return const_col_iterator(m.begin_, m.si1, m.si2, m.ncols()); } 319 static col_iterator col_end(this_type &m) 320 { return col_iterator(m.begin_, m.si1, m.si2, m.ncols()); } 321 static origin_type* origin(this_type &v) { return v.origin; } 322 static const origin_type* origin(const this_type &v) { return v.origin; } 323 static void do_clear(this_type &m) { 324 col_iterator it = mat_col_begin(m), ite = mat_col_end(m); 325 for (; it != ite; ++it) clear(col(it)); 326 } 327 static value_type access(const const_col_iterator &itcol, size_type i) 328 { return linalg_traits<M>::access(*itcol, itcol.si1.index(i)); } 329 static reference access(const col_iterator &itcol, size_type i) 330 { return linalg_traits<M>::access(*itcol, itcol.si1.index(i)); } 331 }; 332 333 template <typename PT, typename SUBI1, typename SUBI2> std::ostream &operator << 334 (std::ostream &o, const gen_sub_col_matrix<PT, SUBI1, SUBI2>& m) 335 { gmm::write(o,m); return o; } 336 337 /* ******************************************************************** */ 338 /* sub matrices */ 339 /* ******************************************************************** */ 340 341 template <typename PT, typename SUBI1, typename SUBI2, typename ST> 342 struct sub_matrix_type_ { 343 typedef abstract_null_type return_type; 344 }; 345 template <typename PT, typename SUBI1, typename SUBI2> 346 struct sub_matrix_type_<PT, SUBI1, SUBI2, col_major> 347 { typedef gen_sub_col_matrix<PT, SUBI1, SUBI2> matrix_type; }; 348 template <typename PT, typename SUBI1, typename SUBI2> 349 struct sub_matrix_type_<PT, SUBI1, SUBI2, row_major> 350 { typedef gen_sub_row_matrix<PT, SUBI1, SUBI2> matrix_type; }; 351 template <typename PT, typename SUBI1, typename SUBI2> 352 struct sub_matrix_type { 353 typedef typename std::iterator_traits<PT>::value_type M; 354 typedef typename sub_matrix_type_<PT, SUBI1, SUBI2, 355 typename principal_orientation_type<typename 356 linalg_traits<M>::sub_orientation>::potype>::matrix_type matrix_type; 357 }; 358 359 template <typename M, typename SUBI1, typename SUBI2> inline 360 typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI2> 361 ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI2>::matrix_type, 362 M *>::return_type 363 sub_matrix(M &m, const SUBI1 &si1, const SUBI2 &si2) { 364 GMM_ASSERT2(si1.last() <= mat_nrows(m) && si2.last() <= mat_ncols(m), 365 "sub matrix too large"); 366 return typename select_return<typename sub_matrix_type<const M *, SUBI1, 367 SUBI2>::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI2> 368 ::matrix_type, M *>::return_type(linalg_cast(m), si1, si2); 369 } 370 371 template <typename M, typename SUBI1, typename SUBI2> inline 372 typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI2> 373 ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI2>::matrix_type, 374 const M *>::return_type 375 sub_matrix(const M &m, const SUBI1 &si1, const SUBI2 &si2) { 376 GMM_ASSERT2(si1.last() <= mat_nrows(m) && si2.last() <= mat_ncols(m), 377 "sub matrix too large"); 378 return typename select_return<typename sub_matrix_type<const M *, SUBI1, 379 SUBI2>::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI2> 380 ::matrix_type, const M *>::return_type(linalg_cast(m), si1, si2); 381 } 382 383 template <typename M, typename SUBI1> inline 384 typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI1> 385 ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type, 386 M *>::return_type 387 sub_matrix(M &m, const SUBI1 &si1) { 388 GMM_ASSERT2(si1.last() <= mat_nrows(m) && si1.last() <= mat_ncols(m), 389 "sub matrix too large"); 390 return typename select_return<typename sub_matrix_type<const M *, SUBI1, 391 SUBI1>::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI1> 392 ::matrix_type, M *>::return_type(linalg_cast(m), si1, si1); 393 } 394 395 template <typename M, typename SUBI1> inline 396 typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI1> 397 ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type, 398 const M *>::return_type 399 sub_matrix(const M &m, const SUBI1 &si1) { 400 GMM_ASSERT2(si1.last() <= mat_nrows(m) && si1.last() <= mat_ncols(m), 401 "sub matrix too large"); 402 return typename select_return<typename sub_matrix_type<const M *, SUBI1, 403 SUBI1>::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI1> 404 ::matrix_type, const M *>::return_type(linalg_cast(m), si1, si1); 405 } 406 407 } 408 409 #endif // GMM_SUB_MATRIX_H__ 410