GetFEM  5.5
gmm_sub_matrix.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-2026 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 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program. If not, see https://www.gnu.org/licenses/.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 
31 /**@file gmm_sub_matrix.h
32  @author Yves Renard <[email protected]>
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>::const_row_iterator,
55  typename linalg_traits<M>::row_iterator,
56  PT>::ref_type
57  iterator;
58  typedef typename linalg_traits<this_type>::reference reference;
59  typedef typename linalg_traits<this_type>::porigin_type porigin_type;
60 
61  SUBI1 si1;
62  SUBI2 si2;
63  iterator begin_;
64  porigin_type origin;
65 
66  reference operator()(size_type i, size_type j) const
67  { return linalg_traits<M>::access(begin_ + si1.index(i), si2.index(j)); }
68 
69  size_type nrows() const { return si1.size(); }
70  size_type ncols() const { return si2.size(); }
71 
72  gen_sub_row_matrix(ref_M m, const SUBI1 &s1, const SUBI2 &s2)
73  : si1(s1), si2(s2), begin_(mat_row_begin(m)),
74  origin(linalg_origin(m)) {}
75  gen_sub_row_matrix() {}
76  gen_sub_row_matrix(const gen_sub_row_matrix<CPT, SUBI1, SUBI2> &cr) :
77  si1(cr.si1), si2(cr.si2), begin_(cr.begin_),origin(cr.origin) {}
78  };
79 
80  template <typename PT, typename SUBI1, typename SUBI2>
81  struct gen_sub_row_matrix_iterator {
82  typedef gen_sub_row_matrix<PT, SUBI1, SUBI2> this_type;
83  typedef typename modifiable_pointer<PT>::pointer MPT;
84  typedef typename const_pointer<PT>::pointer CPT;
85  typedef typename std::iterator_traits<PT>::value_type M;
86  typedef typename select_ref<typename linalg_traits<M>::const_row_iterator,
87  typename linalg_traits<M>::row_iterator,
88  PT>::ref_type
89  value_type;
90  typedef value_type *pointer;
91  typedef ptrdiff_t difference_type;
92  typedef size_t size_type;
93  typedef std::random_access_iterator_tag iterator_category;
94  typedef gen_sub_row_matrix_iterator<PT, SUBI1, SUBI2> iterator;
95 
96  value_type it;
97  SUBI1 si1;
98  SUBI2 si2;
99  size_type ii;
100 
101  iterator operator ++(int) { iterator tmp = *this; ii++; return tmp; }
102  iterator operator --(int) { iterator tmp = *this; ii--; return tmp; }
103  iterator &operator ++() { ii++; return *this; }
104  iterator &operator --() { ii--; return *this; }
105  iterator &operator +=(difference_type i) { ii += i; return *this; }
106  iterator &operator -=(difference_type i) { ii -= i; return *this; }
107  iterator operator +(difference_type i) const
108  { iterator itt = *this; return (itt += i); }
109  iterator operator -(difference_type i) const
110  { iterator itt = *this; return (itt -= i); }
111  difference_type operator -(const iterator &i) const { return ii - i.ii; }
112 
113  value_type operator *() const { return it + si1.index(ii); }
114  value_type operator [](int i) { return it + si1.index(ii+i); }
115 
116  bool operator ==(const iterator &i) const { return (ii == i.ii); }
117  bool operator !=(const iterator &i) const { return !(i == *this); }
118  bool operator < (const iterator &i) const { return (ii < i.ii); }
119  bool operator > (const iterator &i) const { return (ii > i.ii); }
120  bool operator >=(const iterator &i) const { return (ii >= i.ii); }
121 
122  gen_sub_row_matrix_iterator() {}
123  // one of the following two constructors is the copy constructor
124  gen_sub_row_matrix_iterator
125  (const gen_sub_row_matrix_iterator<MPT, SUBI1, SUBI2> &itm)
126  : it(itm.it), si1(itm.si1), si2(itm.si2), ii(itm.ii) {}
127  gen_sub_row_matrix_iterator
128  (const gen_sub_row_matrix_iterator<CPT, SUBI1, SUBI2> &itm)
129  : it(itm.it), si1(itm.si1), si2(itm.si2), ii(itm.ii) {}
130  gen_sub_row_matrix_iterator(const value_type &iter, const SUBI1 &s1,
131  const SUBI2 &s2, size_type i)
132  : it(iter), si1(s1), si2(s2), ii(i) { }
133  };
134 
135  template <typename PT, typename SUBI1, typename SUBI2>
136  struct linalg_traits<gen_sub_row_matrix<PT, SUBI1, SUBI2> > {
137  typedef gen_sub_row_matrix<PT, SUBI1, SUBI2> this_type;
138  typedef typename std::iterator_traits<PT>::value_type M;
139  typedef typename which_reference<PT>::is_reference is_reference;
140  typedef abstract_matrix linalg_type;
141  typedef typename linalg_traits<M>::origin_type origin_type;
142  typedef typename select_ref<const origin_type *,
143  origin_type *,
144  PT>::ref_type
145  porigin_type;
146  typedef typename linalg_traits<M>::value_type value_type;
147  typedef typename select_ref<value_type,
148  typename linalg_traits<M>::reference,
149  PT>::ref_type
150  reference;
151  typedef abstract_null_type sub_col_type;
152  typedef abstract_null_type col_iterator;
153  typedef abstract_null_type const_sub_col_type;
154  typedef abstract_null_type const_col_iterator;
155  typedef typename sub_vector_type
156  <const typename org_type<typename linalg_traits<M>
157  ::const_sub_row_type>::t *,
158  SUBI2>::vector_type
159  const_sub_row_type;
160  typedef typename select_ref
161  <abstract_null_type,
162  typename sub_vector_type<typename org_type
163  <typename linalg_traits<M>
164  ::sub_row_type>::t *,
165  SUBI2>::vector_type,
166  PT>::ref_type
167  sub_row_type;
168  typedef gen_sub_row_matrix_iterator<typename const_pointer<PT>::pointer,
169  SUBI1, SUBI2>
170  const_row_iterator;
171  typedef typename select_ref<abstract_null_type,
172  gen_sub_row_matrix_iterator<PT, SUBI1, SUBI2>,
173  PT>::ref_type
174  row_iterator;
175  typedef typename linalg_traits<const_sub_row_type>::storage_type
176  storage_type;
177  typedef row_major sub_orientation;
178  typedef linalg_true index_sorted;
179 
180  static size_type nrows(const this_type &m) { return m.nrows(); }
181 
182  static size_type ncols(const this_type &m) { return m.ncols(); }
183 
184  static const_sub_row_type row(const const_row_iterator &it)
185  { return const_sub_row_type(linalg_traits<M>::row(*it), it.si2); }
186 
187  static sub_row_type row(const row_iterator &it)
188  { return sub_row_type(linalg_traits<M>::row(*it), it.si2); }
189 
190  static const_row_iterator row_begin(const this_type &m)
191  { return const_row_iterator(m.begin_, m.si1, m.si2, 0); }
192 
193  static row_iterator row_begin(this_type &m)
194  { return row_iterator(m.begin_, m.si1, m.si2, 0); }
195 
196  static const_row_iterator row_end(const this_type &m)
197  { return const_row_iterator(m.begin_, m.si1, m.si2, m.nrows()); }
198 
199  static row_iterator row_end(this_type &m)
200  { return row_iterator(m.begin_, m.si1, m.si2, m.nrows()); }
201 
202  static origin_type* origin(this_type &v) { return v.origin; }
203 
204  static const origin_type* origin(const this_type &v) { return v.origin; }
205 
206  static void do_clear(this_type &m) {
207  row_iterator it = mat_row_begin(m), ite = mat_row_end(m);
208  for (; it != ite; ++it) clear(row(it));
209  }
210 
211  static value_type access(const const_row_iterator &itrow, size_type i)
212  { return linalg_traits<M>::access(*itrow, itrow.si2.index(i)); }
213 
214  static reference access(const row_iterator &itrow, size_type i)
215  { return linalg_traits<M>::access(*itrow, itrow.si2.index(i)); }
216  };
217 
218  template <typename PT, typename SUBI1, typename SUBI2>
219  std::ostream &operator <<(std::ostream &o,
220  const gen_sub_row_matrix<PT, SUBI1, SUBI2>& m)
221  { gmm::write(o,m); return o; }
222 
223 
224  /* ********************************************************************* */
225  /* sub column matrices type */
226  /* ********************************************************************* */
227 
228  template <typename PT, typename SUBI1, typename SUBI2>
229  struct gen_sub_col_matrix {
230  typedef gen_sub_col_matrix<PT, SUBI1, SUBI2> this_type;
231  typedef typename std::iterator_traits<PT>::value_type M;
232  typedef M * CPT;
233  typedef typename std::iterator_traits<PT>::reference ref_M;
234  typedef typename select_ref<typename linalg_traits<M>::const_col_iterator,
235  typename linalg_traits<M>::col_iterator,
236  PT>::ref_type
237  iterator;
238  typedef typename linalg_traits<this_type>::reference reference;
239  typedef typename linalg_traits<this_type>::porigin_type porigin_type;
240 
241  SUBI1 si1;
242  SUBI2 si2;
243  iterator begin_;
244  porigin_type origin;
245 
246  reference operator()(size_type i, size_type j) const
247  { return linalg_traits<M>::access(begin_ + si2.index(j), si1.index(i)); }
248 
249  size_type nrows() const { return si1.size(); }
250  size_type ncols() const { return si2.size(); }
251 
252  gen_sub_col_matrix(ref_M m, const SUBI1 &s1, const SUBI2 &s2)
253  : si1(s1), si2(s2), begin_(mat_col_begin(m)),
254  origin(linalg_origin(m)) {}
255  gen_sub_col_matrix() {}
256  gen_sub_col_matrix(const gen_sub_col_matrix<CPT, SUBI1, SUBI2> &cr) :
257  si1(cr.si1), si2(cr.si2), begin_(cr.begin_),origin(cr.origin) {}
258  };
259 
260  template <typename PT, typename SUBI1, typename SUBI2>
261  struct gen_sub_col_matrix_iterator {
262  typedef gen_sub_col_matrix<PT, SUBI1, SUBI2> this_type;
263  typedef typename modifiable_pointer<PT>::pointer MPT;
264  typedef typename const_pointer<PT>::pointer CPT;
265  typedef typename std::iterator_traits<PT>::value_type M;
266  typedef typename select_ref<typename linalg_traits<M>::const_col_iterator,
267  typename linalg_traits<M>::col_iterator,
268  PT>::ref_type
269  value_type;
270  typedef value_type *pointer;
271  typedef ptrdiff_t difference_type;
272  typedef size_t size_type;
273  typedef std::random_access_iterator_tag iterator_category;
274  typedef gen_sub_col_matrix_iterator<PT, SUBI1, SUBI2> iterator;
275 
276  value_type it;
277  SUBI1 si1;
278  SUBI2 si2;
279  size_type ii;
280 
281  iterator operator ++(int) { iterator tmp = *this; ii++; return tmp; }
282  iterator operator --(int) { iterator tmp = *this; ii--; return tmp; }
283  iterator &operator ++() { ii++; return *this; }
284  iterator &operator --() { ii--; return *this; }
285  iterator &operator +=(difference_type i) { ii += i; return *this; }
286  iterator &operator -=(difference_type i) { ii -= i; return *this; }
287  iterator operator +(difference_type i) const
288  { iterator itt = *this; return (itt += i); }
289  iterator operator -(difference_type i) const
290  { iterator itt = *this; return (itt -= i); }
291  difference_type operator -(const iterator &i) const { return ii - i.ii; }
292 
293  value_type operator *() const { return it + si2.index(ii); }
294  value_type operator [](int i) { return it + si2.index(ii+i); }
295 
296  bool operator ==(const iterator &i) const { return (ii == i.ii); }
297  bool operator !=(const iterator &i) const { return !(i == *this); }
298  bool operator < (const iterator &i) const { return (ii < i.ii); }
299  bool operator > (const iterator &i) const { return (ii > i.ii); }
300  bool operator >=(const iterator &i) const { return (ii >= i.ii); }
301 
302  gen_sub_col_matrix_iterator() {}
303  // one of the following two constructors is the copy constructor
304  gen_sub_col_matrix_iterator
305  (const gen_sub_col_matrix_iterator<MPT, SUBI1, SUBI2> &itm)
306  : it(itm.it), si1(itm.si1), si2(itm.si2), ii(itm.ii) {}
307  gen_sub_col_matrix_iterator
308  (const gen_sub_col_matrix_iterator<CPT, SUBI1, SUBI2> &itm)
309  : it(itm.it), si1(itm.si1), si2(itm.si2), ii(itm.ii) {}
310  gen_sub_col_matrix_iterator(const value_type &iter, const SUBI1 &s1,
311  const SUBI2 &s2, size_type i)
312  : it(iter), si1(s1), si2(s2), ii(i) { }
313  };
314 
315  template <typename PT, typename SUBI1, typename SUBI2>
316  struct linalg_traits<gen_sub_col_matrix<PT, SUBI1, SUBI2> > {
317  typedef gen_sub_col_matrix<PT, SUBI1, SUBI2> this_type;
318  typedef typename std::iterator_traits<PT>::value_type M;
319  typedef typename which_reference<PT>::is_reference is_reference;
320  typedef abstract_matrix linalg_type;
321  typedef typename linalg_traits<M>::origin_type origin_type;
322  typedef typename select_ref<const origin_type *,
323  origin_type *,
324  PT>::ref_type
325  porigin_type;
326  typedef typename linalg_traits<M>::value_type value_type;
327  typedef typename select_ref<value_type,
328  typename linalg_traits<M>::reference,
329  PT>::ref_type
330  reference;
331  typedef abstract_null_type sub_row_type;
332  typedef abstract_null_type row_iterator;
333  typedef abstract_null_type const_sub_row_type;
334  typedef abstract_null_type const_row_iterator;
335  typedef typename sub_vector_type
336  <const typename org_type<typename linalg_traits<M>
337  ::const_sub_col_type>::t *,
338  SUBI1>::vector_type
339  const_sub_col_type;
340  typedef typename select_ref
341  <abstract_null_type,
342  typename sub_vector_type<typename org_type
343  <typename linalg_traits<M>
344  ::sub_col_type>::t *,
345  SUBI1>::vector_type,
346  PT>::ref_type
347  sub_col_type;
348  typedef gen_sub_col_matrix_iterator<typename const_pointer<PT>::pointer,
349  SUBI1, SUBI2>
350  const_col_iterator;
351  typedef typename select_ref<abstract_null_type,
352  gen_sub_col_matrix_iterator<PT, SUBI1, SUBI2>,
353  PT>::ref_type
354  col_iterator;
355  typedef typename linalg_traits<const_sub_col_type>::storage_type
356  storage_type;
357  typedef col_major sub_orientation;
358  typedef linalg_true index_sorted;
359 
360  static size_type nrows(const this_type &m) { return m.nrows(); }
361 
362  static size_type ncols(const this_type &m) { return m.ncols(); }
363 
364  static const_sub_col_type col(const const_col_iterator &it)
365  { return const_sub_col_type(linalg_traits<M>::col(*it), it.si1); }
366 
367  static sub_col_type col(const col_iterator &it)
368  { return sub_col_type(linalg_traits<M>::col(*it), it.si1); }
369 
370  static const_col_iterator col_begin(const this_type &m)
371  { return const_col_iterator(m.begin_, m.si1, m.si2, 0); }
372 
373  static col_iterator col_begin(this_type &m)
374  { return col_iterator(m.begin_, m.si1, m.si2, 0); }
375 
376  static const_col_iterator col_end(const this_type &m)
377  { return const_col_iterator(m.begin_, m.si1, m.si2, m.ncols()); }
378 
379  static col_iterator col_end(this_type &m)
380  { return col_iterator(m.begin_, m.si1, m.si2, m.ncols()); }
381 
382  static origin_type* origin(this_type &v) { return v.origin; }
383 
384  static const origin_type* origin(const this_type &v) { return v.origin; }
385 
386  static void do_clear(this_type &m) {
387  col_iterator it = mat_col_begin(m), ite = mat_col_end(m);
388  for (; it != ite; ++it) clear(col(it));
389  }
390 
391  static value_type access(const const_col_iterator &itcol, size_type i)
392  { return linalg_traits<M>::access(*itcol, itcol.si1.index(i)); }
393 
394  static reference access(const col_iterator &itcol, size_type i)
395  { return linalg_traits<M>::access(*itcol, itcol.si1.index(i)); }
396  };
397 
398  template <typename PT, typename SUBI1, typename SUBI2>
399  std::ostream &operator <<(std::ostream &o,
400  const gen_sub_col_matrix<PT, SUBI1, SUBI2>& m)
401  { gmm::write(o,m); return o; }
402 
403 
404  /* ******************************************************************** */
405  /* sub matrices */
406  /* ******************************************************************** */
407 
408  template <typename PT, typename SUBI1, typename SUBI2, typename ST>
409  struct sub_matrix_type_ {
410  typedef abstract_null_type return_type;
411  };
412 
413  template <typename PT, typename SUBI1, typename SUBI2>
414  struct sub_matrix_type_<PT, SUBI1, SUBI2, col_major>
415  { typedef gen_sub_col_matrix<PT, SUBI1, SUBI2> matrix_type; };
416 
417  template <typename PT, typename SUBI1, typename SUBI2>
418  struct sub_matrix_type_<PT, SUBI1, SUBI2, row_major>
419  { typedef gen_sub_row_matrix<PT, SUBI1, SUBI2> matrix_type; };
420 
421  template <typename PT, typename SUBI1, typename SUBI2>
422  struct sub_matrix_type {
423  typedef typename std::iterator_traits<PT>::value_type M;
424  typedef typename sub_matrix_type_<PT, SUBI1, SUBI2,
425  typename principal_orientation_type<typename
426  linalg_traits<M>::sub_orientation>::potype>::matrix_type matrix_type;
427  };
428 
429  template <typename M, typename SUBI1, typename SUBI2> inline
430  typename select_return
431  <typename sub_matrix_type<const M *, SUBI1, SUBI2>::matrix_type,
432  typename sub_matrix_type<M *, SUBI1, SUBI2>::matrix_type,
433  M *>::return_type
434  sub_matrix(M &m, const SUBI1 &si1, const SUBI2 &si2) {
435  GMM_ASSERT2(si1.last() <= mat_nrows(m) && si2.last() <= mat_ncols(m),
436  "sub matrix too large");
437  return
438  typename select_return
439  <typename sub_matrix_type<const M *, SUBI1,SUBI2>::matrix_type,
440  typename sub_matrix_type<M *, SUBI1, SUBI2>::matrix_type,
441  M *>::return_type(linalg_cast(m), si1, si2);
442  }
443 
444  template <typename M, typename SUBI1, typename SUBI2> inline
445  typename select_return
446  <typename sub_matrix_type<const M *, SUBI1, SUBI2>::matrix_type,
447  typename sub_matrix_type<M *, SUBI1, SUBI2>::matrix_type,
448  const M *>::return_type
449  sub_matrix(const M &m, const SUBI1 &si1, const SUBI2 &si2) {
450  GMM_ASSERT2(si1.last() <= mat_nrows(m) && si2.last() <= mat_ncols(m),
451  "sub matrix too large");
452  return
453  typename select_return
454  <typename sub_matrix_type<const M *, SUBI1, SUBI2>::matrix_type,
455  typename sub_matrix_type<M *, SUBI1, SUBI2>::matrix_type,
456  const M *>::return_type(linalg_cast(m), si1, si2);
457  }
458 
459  template <typename M, typename SUBI1> inline
460  typename select_return
461  <typename sub_matrix_type<const M *, SUBI1, SUBI1>::matrix_type,
462  typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type,
463  M *>::return_type
464  sub_matrix(M &m, const SUBI1 &si1) {
465  GMM_ASSERT2(si1.last() <= mat_nrows(m) && si1.last() <= mat_ncols(m),
466  "sub matrix too large");
467  return
468  typename select_return
469  <typename sub_matrix_type<const M *, SUBI1, SUBI1>::matrix_type,
470  typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type,
471  M *>::return_type(linalg_cast(m), si1, si1);
472  }
473 
474  template <typename M, typename SUBI1> inline
475  typename select_return
476  <typename sub_matrix_type<const M *, SUBI1, SUBI1>::matrix_type,
477  typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type,
478  const M *>::return_type
479  sub_matrix(const M &m, const SUBI1 &si1) {
480  GMM_ASSERT2(si1.last() <= mat_nrows(m) && si1.last() <= mat_ncols(m),
481  "sub matrix too large");
482  return
483  typename select_return
484  <typename sub_matrix_type<const M *, SUBI1, SUBI1>::matrix_type,
485  typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type,
486  const M *>::return_type(linalg_cast(m), si1, si1);
487  }
488 
489 }
490 
491 #endif // GMM_SUB_MATRIX_H__
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
gmm_sub_vector.h
Generic sub-vectors.
bgeot::operator-
rational_fraction< T > operator-(const polynomial< T > &P, const rational_fraction< T > &Q)
Subtract Q from P.
Definition: bgeot_poly.h:755
bgeot::operator+
rational_fraction< T > operator+(const polynomial< T > &P, const rational_fraction< T > &Q)
Add Q to P.
Definition: bgeot_poly.h:748
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48