GetFEM  5.5
gmm_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_matrix.h
32  @author Yves Renard <Yves.Renard@insa-lyon.fr>
33  @date October 13, 2002.
34  @brief Declaration of some matrix types (gmm::dense_matrix,
35  gmm::row_matrix, gmm::col_matrix, gmm::csc_matrix, etc.)
36 */
37 
38 #ifndef GMM_MATRIX_H__
39 #define GMM_MATRIX_H__
40 
41 #include "gmm_vector.h"
42 #include "gmm_sub_vector.h"
43 #include "gmm_sub_matrix.h"
44 #include "gmm_transposed.h"
45 
46 namespace gmm
47 {
48 
49  /* ******************************************************************** */
50  /* */
51  /* Identity matrix */
52  /* */
53  /* ******************************************************************** */
54 
55  struct identity_matrix {
56  template <class MAT> void build_with(const MAT &) {}
57  };
58 
59  template <typename M> inline
60  void add(const identity_matrix&, M &v1) {
61  size_type n = std::min(gmm::mat_nrows(v1), gmm::mat_ncols(v1));
62  for (size_type i = 0; i < n; ++i)
63  v1(i,i) += typename linalg_traits<M>::value_type(1);
64  }
65  template <typename M> inline
66  void add(const identity_matrix &II, const M &v1)
67  { add(II, linalg_const_cast(v1)); }
68 
69  template <typename V1, typename V2> inline
70  void mult(const identity_matrix&, const V1 &v1, V2 &v2)
71  { copy(v1, v2); }
72  template <typename V1, typename V2> inline
73  void mult(const identity_matrix&, const V1 &v1, const V2 &v2)
74  { copy(v1, v2); }
75  template <typename V1, typename V2, typename V3> inline
76  void mult(const identity_matrix&, const V1 &v1, const V2 &v2, V3 &v3)
77  { add(v1, v2, v3); }
78  template <typename V1, typename V2, typename V3> inline
79  void mult(const identity_matrix&, const V1 &v1, const V2 &v2, const V3 &v3)
80  { add(v1, v2, v3); }
81  template <typename V1, typename V2> inline
82  void left_mult(const identity_matrix&, const V1 &v1, V2 &v2)
83  { copy(v1, v2); }
84  template <typename V1, typename V2> inline
85  void left_mult(const identity_matrix&, const V1 &v1, const V2 &v2)
86  { copy(v1, v2); }
87  template <typename V1, typename V2> inline
88  void right_mult(const identity_matrix&, const V1 &v1, V2 &v2)
89  { copy(v1, v2); }
90  template <typename V1, typename V2> inline
91  void right_mult(const identity_matrix&, const V1 &v1, const V2 &v2)
92  { copy(v1, v2); }
93  template <typename V1, typename V2> inline
94  void transposed_left_mult(const identity_matrix&, const V1 &v1, V2 &v2)
95  { copy(v1, v2); }
96  template <typename V1, typename V2> inline
97  void transposed_left_mult(const identity_matrix&, const V1 &v1,const V2 &v2)
98  { copy(v1, v2); }
99  template <typename V1, typename V2> inline
100  void transposed_right_mult(const identity_matrix&, const V1 &v1, V2 &v2)
101  { copy(v1, v2); }
102  template <typename V1, typename V2> inline
103  void transposed_right_mult(const identity_matrix&,const V1 &v1,const V2 &v2)
104  { copy(v1, v2); }
105  template <typename M> void copy_ident(const identity_matrix&, M &m) {
106  size_type i = 0, n = std::min(mat_nrows(m), mat_ncols(m));
107  clear(m);
108  for (; i < n; ++i) m(i,i) = typename linalg_traits<M>::value_type(1);
109  }
110  template <typename M> inline void copy(const identity_matrix&, M &m)
111  { copy_ident(identity_matrix(), m); }
112  template <typename M> inline void copy(const identity_matrix &, const M &m)
113  { copy_ident(identity_matrix(), linalg_const_cast(m)); }
114  template <typename V1, typename V2> inline
115  typename linalg_traits<V1>::value_type
116  vect_sp(const identity_matrix &, const V1 &v1, const V2 &v2)
117  { return vect_sp(v1, v2); }
118  template <typename V1, typename V2> inline
119  typename linalg_traits<V1>::value_type
120  vect_hp(const identity_matrix &, const V1 &v1, const V2 &v2)
121  { return vect_hp(v1, v2); }
122  template<typename M> inline bool is_identity(const M&) { return false; }
123  inline bool is_identity(const identity_matrix&) { return true; }
124 
125  /* ******************************************************************** */
126  /* */
127  /* Row matrix */
128  /* */
129  /* ******************************************************************** */
130 
131  template<typename V> class row_matrix {
132  protected :
133  std::vector<V> li; /* array of rows. */
134  size_type nc;
135 
136  public :
137 
138  typedef typename linalg_traits<V>::reference reference;
139  typedef typename linalg_traits<V>::value_type value_type;
140 
141  row_matrix(size_type r, size_type c) : li(r, V(c)), nc(c) {}
142  row_matrix() : nc(0) {}
143  reference operator ()(size_type l, size_type c)
144  { return li[l][c]; }
145  value_type operator ()(size_type l, size_type c) const
146  { return li[l][c]; }
147 
148  void clear_mat();
149  void resize(size_type m, size_type n);
150 
151  typename std::vector<V>::iterator begin()
152  { return li.begin(); }
153  typename std::vector<V>::iterator end()
154  { return li.end(); }
155  typename std::vector<V>::const_iterator begin() const
156  { return li.begin(); }
157  typename std::vector<V>::const_iterator end() const
158  { return li.end(); }
159 
160 
161  V& row(size_type i) { return li[i]; }
162  const V& row(size_type i) const { return li[i]; }
163  V& operator[](size_type i) { return li[i]; }
164  const V& operator[](size_type i) const { return li[i]; }
165 
166  inline size_type nrows() const { return li.size(); }
167  inline size_type ncols() const { return nc; }
168 
169  void swap(row_matrix<V> &m) { std::swap(li, m.li); std::swap(nc, m.nc); }
170  void swap_row(size_type i, size_type j) { std::swap(li[i], li[j]); }
171  };
172 
173  template<typename V>
174  void row_matrix<V>::resize(size_type m, size_type n) {
175  size_type nr = std::min(nrows(), m);
176  li.resize(m);
177  for (size_type i=nr; i < m; ++i) gmm::resize(li[i], n);
178  if (n != nc) {
179  for (size_type i=0; i < nr; ++i) gmm::resize(li[i], n);
180  nc = n;
181  }
182  }
183 
184 
185  template<typename V>
186  void row_matrix<V>::clear_mat()
187  { for (size_type i=0; i < nrows(); ++i) clear(li[i]); }
188 
189  template <typename V>
190  struct linalg_traits<row_matrix<V> > {
191  typedef row_matrix<V> this_type;
192  typedef this_type origin_type;
193  typedef linalg_false is_reference;
194  typedef abstract_matrix linalg_type;
195  typedef typename linalg_traits<V>::value_type value_type;
196  typedef typename linalg_traits<V>::reference reference;
197  typedef typename linalg_traits<V>::storage_type storage_type;
198  typedef V & sub_row_type;
199  typedef const V & const_sub_row_type;
200  typedef typename std::vector<V>::iterator row_iterator;
201  typedef typename std::vector<V>::const_iterator const_row_iterator;
202  typedef abstract_null_type sub_col_type;
203  typedef abstract_null_type const_sub_col_type;
204  typedef abstract_null_type col_iterator;
205  typedef abstract_null_type const_col_iterator;
206  typedef row_major sub_orientation;
207  typedef linalg_true index_sorted;
208  static size_type nrows(const this_type &m) { return m.nrows(); }
209  static size_type ncols(const this_type &m) { return m.ncols(); }
210  static row_iterator row_begin(this_type &m) { return m.begin(); }
211  static row_iterator row_end(this_type &m) { return m.end(); }
212  static const_row_iterator row_begin(const this_type &m)
213  { return m.begin(); }
214  static const_row_iterator row_end(const this_type &m)
215  { return m.end(); }
216  static const_sub_row_type row(const const_row_iterator &it)
217  { return const_sub_row_type(*it); }
218  static sub_row_type row(const row_iterator &it)
219  { return sub_row_type(*it); }
220  static origin_type* origin(this_type &m) { return &m; }
221  static const origin_type* origin(const this_type &m) { return &m; }
222  static void do_clear(this_type &m) { m.clear_mat(); }
223  static value_type access(const const_row_iterator &itrow, size_type j)
224  { return (*itrow)[j]; }
225  static reference access(const row_iterator &itrow, size_type j)
226  { return (*itrow)[j]; }
227  static void resize(this_type &v, size_type m, size_type n)
228  { v.resize(m, n); }
229  static void reshape(this_type &, size_type, size_type)
230  { GMM_ASSERT1(false, "Sorry, to be done"); }
231  };
232 
233  template<typename V> std::ostream &operator <<
234  (std::ostream &o, const row_matrix<V>& m) { gmm::write(o,m); return o; }
235 
236  /* ******************************************************************** */
237  /* */
238  /* Column matrix */
239  /* */
240  /* ******************************************************************** */
241 
242  template<typename V> class col_matrix {
243  protected :
244  std::vector<V> li; /* array of columns. */
245  size_type nr;
246 
247  public :
248 
249  typedef typename linalg_traits<V>::reference reference;
250  typedef typename linalg_traits<V>::value_type value_type;
251 
252  col_matrix(size_type r, size_type c) : li(c, V(r)), nr(r) { }
253  col_matrix() : nr(0) {}
254  reference operator ()(size_type l, size_type c)
255  { return li[c][l]; }
256  value_type operator ()(size_type l, size_type c) const
257  { return li[c][l]; }
258 
259  void clear_mat();
260  void resize(size_type, size_type);
261 
262  V& col(size_type i) { return li[i]; }
263  const V& col(size_type i) const { return li[i]; }
264  V& operator[](size_type i) { return li[i]; }
265  const V& operator[](size_type i) const { return li[i]; }
266 
267  typename std::vector<V>::iterator begin()
268  { return li.begin(); }
269  typename std::vector<V>::iterator end()
270  { return li.end(); }
271  typename std::vector<V>::const_iterator begin() const
272  { return li.begin(); }
273  typename std::vector<V>::const_iterator end() const
274  { return li.end(); }
275 
276  inline size_type ncols() const { return li.size(); }
277  inline size_type nrows() const { return nr; }
278 
279  void swap(col_matrix<V> &m) { std::swap(li, m.li); std::swap(nr, m.nr); }
280  void swap_col(size_type i, size_type j) { std::swap(li[i], li[j]); }
281  };
282 
283  template<typename V> void col_matrix<V>::resize(size_type m, size_type n) {
284  size_type nc = std::min(ncols(), n);
285  li.resize(n);
286  for (size_type i=nc; i < n; ++i) gmm::resize(li[i], m);
287  if (m != nr) {
288  for (size_type i=0; i < nc; ++i) gmm::resize(li[i], m);
289  nr = m;
290  }
291  }
292 
293  template<typename V> void col_matrix<V>::clear_mat()
294  { for (size_type i=0; i < ncols(); ++i) clear(li[i]); }
295 
296  template <typename V> struct linalg_traits<col_matrix<V> > {
297  typedef col_matrix<V> this_type;
298  typedef this_type origin_type;
299  typedef linalg_false is_reference;
300  typedef abstract_matrix linalg_type;
301  typedef typename linalg_traits<V>::value_type value_type;
302  typedef typename linalg_traits<V>::reference reference;
303  typedef typename linalg_traits<V>::storage_type storage_type;
304  typedef V &sub_col_type;
305  typedef const V &const_sub_col_type;
306  typedef typename std::vector<V>::iterator col_iterator;
307  typedef typename std::vector<V>::const_iterator const_col_iterator;
308  typedef abstract_null_type sub_row_type;
309  typedef abstract_null_type const_sub_row_type;
310  typedef abstract_null_type row_iterator;
311  typedef abstract_null_type const_row_iterator;
312  typedef col_major sub_orientation;
313  typedef linalg_true index_sorted;
314  static size_type nrows(const this_type &m) { return m.nrows(); }
315  static size_type ncols(const this_type &m) { return m.ncols(); }
316  static col_iterator col_begin(this_type &m) { return m.begin(); }
317  static col_iterator col_end(this_type &m) { return m.end(); }
318  static const_col_iterator col_begin(const this_type &m)
319  { return m.begin(); }
320  static const_col_iterator col_end(const this_type &m)
321  { return m.end(); }
322  static const_sub_col_type col(const const_col_iterator &it)
323  { return *it; }
324  static sub_col_type col(const col_iterator &it)
325  { return *it; }
326  static origin_type* origin(this_type &m) { return &m; }
327  static const origin_type* origin(const this_type &m) { return &m; }
328  static void do_clear(this_type &m) { m.clear_mat(); }
329  static value_type access(const const_col_iterator &itcol, size_type j)
330  { return (*itcol)[j]; }
331  static reference access(const col_iterator &itcol, size_type j)
332  { return (*itcol)[j]; }
333  static void resize(this_type &v, size_type m, size_type n)
334  { v.resize(m,n); }
335  static void reshape(this_type &, size_type, size_type)
336  { GMM_ASSERT1(false, "Sorry, to be done"); }
337  };
338 
339  template<typename V> std::ostream &operator <<
340  (std::ostream &o, const col_matrix<V>& m) { gmm::write(o,m); return o; }
341 
342  /* ******************************************************************** */
343  /* */
344  /* Dense matrix */
345  /* */
346  /* ******************************************************************** */
347 
348  template<typename T> class dense_matrix : public std::vector<T> {
349  public:
350  typedef typename std::vector<T>::size_type size_type;
351  typedef typename std::vector<T>::iterator iterator;
352  typedef typename std::vector<T>::const_iterator const_iterator;
353  typedef typename std::vector<T>::reference reference;
354  typedef typename std::vector<T>::const_reference const_reference;
355 
356  protected:
357  size_type nbc, nbl;
358 
359  public:
360 
361  inline const_reference operator ()(size_type l, size_type c) const {
362  GMM_ASSERT2(l < nbl && c < nbc, "out of range");
363  return *(this->begin() + c*nbl+l);
364  }
365  inline reference operator ()(size_type l, size_type c) {
366  GMM_ASSERT2(l < nbl && c < nbc, "out of range");
367  return *(this->begin() + c*nbl+l);
368  }
369 
370  std::vector<T> &as_vector() { return *this; }
371  const std::vector<T> &as_vector() const { return *this; }
372 
373  void resize(size_type, size_type);
374  void base_resize(size_type, size_type);
375  void reshape(size_type, size_type);
376 
377  void fill(T a, T b = T(0));
378  inline size_type nrows() const { return nbl; }
379  inline size_type ncols() const { return nbc; }
380  void swap(dense_matrix<T> &m)
381  { std::vector<T>::swap(m); std::swap(nbc, m.nbc); std::swap(nbl, m.nbl); }
382 
383  dense_matrix(size_type l, size_type c)
384  : std::vector<T>(c*l), nbc(c), nbl(l) {}
385  dense_matrix() { nbl = nbc = 0; }
386  };
387 
388  template<typename T>
389  void dense_matrix<T>::reshape(size_type m,size_type n) {
390  GMM_ASSERT2(n*m == nbl*nbc, "dimensions mismatch");
391  nbl = m; nbc = n;
392  }
393 
394  template<typename T>
395  void dense_matrix<T>::base_resize(size_type m, size_type n)
396  { std::vector<T>::resize(n*m); nbl = m; nbc = n; }
397 
398  template<typename T>
399  void dense_matrix<T>::resize(size_type m, size_type n) {
400  if (n*m > nbc*nbl) std::vector<T>::resize(n*m);
401  if (m < nbl) {
402  for (size_type i = 1; i < std::min(nbc, n); ++i)
403  std::copy(this->begin()+i*nbl, this->begin()+(i*nbl+m),
404  this->begin()+i*m);
405  for (size_type i = std::min(nbc, n); i < n; ++i)
406  std::fill(this->begin()+(i*m), this->begin()+(i+1)*m, T(0));
407  }
408  else if (m > nbl) { /* do nothing when the nb of rows does not change */
409  for (size_type i = std::min(nbc, n); i > 1; --i)
410  std::copy(this->begin()+(i-1)*nbl, this->begin()+i*nbl,
411  this->begin()+(i-1)*m);
412  for (size_type i = 0; i < std::min(nbc, n); ++i)
413  std::fill(this->begin()+(i*m+nbl), this->begin()+(i+1)*m, T(0));
414  }
415  if (n*m < nbc*nbl) std::vector<T>::resize(n*m);
416  nbl = m; nbc = n;
417  }
418 
419  template<typename T>
420  void dense_matrix<T>::fill(T a, T b) {
421  std::fill(this->begin(), this->end(), b);
422  size_type n = std::min(nbl, nbc);
423  if (a != b) for (size_type i = 0; i < n; ++i) (*this)(i,i) = a;
424  }
425 
426  template <typename T>
427  struct linalg_traits<dense_matrix<T> > {
428  typedef dense_matrix<T> this_type;
429  typedef this_type origin_type;
430  typedef linalg_false is_reference;
431  typedef abstract_matrix linalg_type;
432  typedef T value_type;
433  typedef T& reference;
434  typedef abstract_dense storage_type;
435  typedef tab_ref_reg_spaced_with_origin<typename this_type::iterator,
436  this_type> sub_row_type;
437  typedef tab_ref_reg_spaced_with_origin<typename this_type::const_iterator,
438  this_type> const_sub_row_type;
439  typedef dense_compressed_iterator<typename this_type::iterator,
440  typename this_type::iterator,
441  this_type *> row_iterator;
442  typedef dense_compressed_iterator<typename this_type::const_iterator,
443  typename this_type::iterator,
444  const this_type *> const_row_iterator;
445  typedef tab_ref_with_origin<typename this_type::iterator,
446  this_type> sub_col_type;
447  typedef tab_ref_with_origin<typename this_type::const_iterator,
448  this_type> const_sub_col_type;
449  typedef dense_compressed_iterator<typename this_type::iterator,
450  typename this_type::iterator,
451  this_type *> col_iterator;
452  typedef dense_compressed_iterator<typename this_type::const_iterator,
453  typename this_type::iterator,
454  const this_type *> const_col_iterator;
455  typedef col_and_row sub_orientation;
456  typedef linalg_true index_sorted;
457  static size_type nrows(const this_type &m) { return m.nrows(); }
458  static size_type ncols(const this_type &m) { return m.ncols(); }
459  static const_sub_row_type row(const const_row_iterator &it)
460  { return const_sub_row_type(*it, it.nrows, it.ncols, it.origin); }
461  static const_sub_col_type col(const const_col_iterator &it)
462  { return const_sub_col_type(*it, *it + it.nrows, it.origin); }
463  static sub_row_type row(const row_iterator &it)
464  { return sub_row_type(*it, it.nrows, it.ncols, it.origin); }
465  static sub_col_type col(const col_iterator &it)
466  { return sub_col_type(*it, *it + it.nrows, it.origin); }
467  static row_iterator row_begin(this_type &m)
468  { return row_iterator(m.begin(), m.size() ? 1 : 0, m.nrows(), m.ncols(), 0, &m); }
469  static row_iterator row_end(this_type &m)
470  { return row_iterator(m.begin(), m.size() ? 1 : 0, m.nrows(), m.ncols(), m.nrows(), &m); }
471  static const_row_iterator row_begin(const this_type &m)
472  { return const_row_iterator(m.begin(), m.size() ? 1 : 0, m.nrows(), m.ncols(), 0, &m); }
473  static const_row_iterator row_end(const this_type &m)
474  { return const_row_iterator(m.begin(), m.size() ? 1 : 0, m.nrows(), m.ncols(), m.nrows(), &m); }
475  static col_iterator col_begin(this_type &m)
476  { return col_iterator(m.begin(), m.nrows(), m.nrows(), m.ncols(), 0, &m); }
477  static col_iterator col_end(this_type &m)
478  { return col_iterator(m.begin(), m.nrows(), m.nrows(), m.ncols(), m.ncols(), &m); }
479  static const_col_iterator col_begin(const this_type &m)
480  { return const_col_iterator(m.begin(), m.nrows(), m.nrows(), m.ncols(), 0, &m); }
481  static const_col_iterator col_end(const this_type &m)
482  { return const_col_iterator(m.begin(),m.nrows(),m.nrows(),m.ncols(),m.ncols(), &m); }
483  static origin_type* origin(this_type &m) { return &m; }
484  static const origin_type* origin(const this_type &m) { return &m; }
485  static void do_clear(this_type &m) { m.fill(value_type(0)); }
486  static value_type access(const const_col_iterator &itcol, size_type j)
487  { return (*itcol)[j]; }
488  static reference access(const col_iterator &itcol, size_type j)
489  { return (*itcol)[j]; }
490  static void resize(this_type &v, size_type m, size_type n)
491  { v.resize(m,n); }
492  static void reshape(this_type &v, size_type m, size_type n)
493  { v.reshape(m, n); }
494  };
495 
496  template<typename T> std::ostream &operator <<
497  (std::ostream &o, const dense_matrix<T>& m) { gmm::write(o,m); return o; }
498 
499 
500  /* ******************************************************************** */
501  /* */
502  /* Read only compressed sparse column matrix */
503  /* */
504  /* ******************************************************************** */
505 
506  template <typename T, typename IND_TYPE = unsigned int, int shift = 0>
507  struct csc_matrix {
508 
509  std::vector<T> pr;
510  std::vector<IND_TYPE> ir;
511  std::vector<IND_TYPE> jc;
512  size_type nc, nr;
513 
514  typedef T value_type;
515  typedef T& access_type;
516 
517  template <typename Matrix> void init_with_good_format(const Matrix &B);
518  template <typename Matrix> void init_with(const Matrix &A);
519  void init_with(const col_matrix<gmm::rsvector<T> > &B)
520  { init_with_good_format(B); }
521  void init_with(const col_matrix<wsvector<T> > &B)
522  { init_with_good_format(B); }
523  template <typename PT1, typename PT2, typename PT3, int cshift>
524  void init_with(const csc_matrix_ref<PT1,PT2,PT3,cshift>& B)
525  { init_with_good_format(B); }
526  template <typename U, int cshift>
527  void init_with(const csc_matrix<U, IND_TYPE, cshift>& B)
528  { init_with_good_format(B); }
529 
530  void init_with_identity(size_type n);
531 
532  csc_matrix() : nc(0), nr(0) {}
533  csc_matrix(size_type nnr, size_type nnc);
534 
535  size_type nrows() const { return nr; }
536  size_type ncols() const { return nc; }
537  void swap(csc_matrix<T, IND_TYPE, shift> &m) {
538  std::swap(pr, m.pr);
539  std::swap(ir, m.ir); std::swap(jc, m.jc);
540  std::swap(nc, m.nc); std::swap(nr, m.nr);
541  }
542  value_type operator()(size_type i, size_type j) const
543  { return mat_col(*this, j)[i]; }
544  };
545 
546  template <typename T, typename IND_TYPE, int shift> template<typename Matrix>
547  void csc_matrix<T, IND_TYPE, shift>::init_with_good_format(const Matrix &B) {
548  typedef typename linalg_traits<Matrix>::const_sub_col_type col_type;
549  nc = mat_ncols(B); nr = mat_nrows(B);
550  jc.resize(nc+1);
551  jc[0] = shift;
552  for (size_type j = 0; j < nc; ++j) {
553  jc[j+1] = IND_TYPE(jc[j] + nnz(mat_const_col(B, j)));
554  }
555  pr.resize(jc[nc]);
556  ir.resize(jc[nc]);
557  for (size_type j = 0; j < nc; ++j) {
558  col_type col = mat_const_col(B, j);
559  typename linalg_traits<typename org_type<col_type>::t>::const_iterator
560  it = vect_const_begin(col), ite = vect_const_end(col);
561  for (size_type k = 0; it != ite; ++it, ++k) {
562  pr[jc[j]-shift+k] = *it;
563  ir[jc[j]-shift+k] = IND_TYPE(it.index() + shift);
564  }
565  }
566  }
567 
568  template <typename T, typename IND_TYPE, int shift>
569  template <typename Matrix>
570  void csc_matrix<T, IND_TYPE, shift>::init_with(const Matrix &A) {
571  col_matrix<wsvector<T> > B(mat_nrows(A), mat_ncols(A));
572  copy(A, B);
573  init_with_good_format(B);
574  }
575 
576  template <typename T, typename IND_TYPE, int shift>
577  void csc_matrix<T, IND_TYPE, shift>::init_with_identity(size_type n) {
578  nc = nr = n;
579  pr.resize(nc); ir.resize(nc); jc.resize(nc+1);
580  for (size_type j = 0; j < nc; ++j)
581  { ir[j] = jc[j] = shift + j; pr[j] = T(1); }
582  jc[nc] = shift + nc;
583  }
584 
585  template <typename T, typename IND_TYPE, int shift>
586  csc_matrix<T, IND_TYPE, shift>::csc_matrix(size_type nnr, size_type nnc)
587  : nc(nnc), nr(nnr) {
588  pr.resize(1); ir.resize(1); jc.resize(nc+1);
589  for (size_type j = 0; j <= nc; ++j) jc[j] = shift;
590  }
591 
592  template <typename T, typename IND_TYPE, int shift>
593  struct linalg_traits<csc_matrix<T, IND_TYPE, shift> > {
594  typedef csc_matrix<T, IND_TYPE, shift> this_type;
595  typedef linalg_const is_reference;
596  typedef abstract_matrix linalg_type;
597  typedef T value_type;
598  typedef T origin_type;
599  typedef T reference;
600  typedef abstract_sparse storage_type;
601  typedef abstract_null_type sub_row_type;
602  typedef abstract_null_type const_sub_row_type;
603  typedef abstract_null_type row_iterator;
604  typedef abstract_null_type const_row_iterator;
605  typedef abstract_null_type sub_col_type;
606  typedef cs_vector_ref<const T *, const IND_TYPE *, shift>
607  const_sub_col_type;
608  typedef sparse_compressed_iterator<const T *, const IND_TYPE *,
609  const IND_TYPE *, shift>
610  const_col_iterator;
611  typedef abstract_null_type col_iterator;
612  typedef col_major sub_orientation;
613  typedef linalg_true index_sorted;
614  static size_type nrows(const this_type &m) { return m.nrows(); }
615  static size_type ncols(const this_type &m) { return m.ncols(); }
616  static const_col_iterator col_begin(const this_type &m)
617  { return const_col_iterator(&m.pr[0],&m.ir[0],&m.jc[0], m.nr, &m.pr[0]); }
618  static const_col_iterator col_end(const this_type &m) {
619  return const_col_iterator(&m.pr[0],&m.ir[0],&m.jc[0]+m.nc,
620  m.nr,&m.pr[0]);
621  }
622  static const_sub_col_type col(const const_col_iterator &it) {
623  return const_sub_col_type(it.pr + *(it.jc) - shift,
624  it.ir + *(it.jc) - shift,
625  *(it.jc + 1) - *(it.jc), it.n);
626  }
627  static const origin_type* origin(const this_type &m) { return &m.pr[0]; }
628  static void do_clear(this_type &m) { m.do_clear(); }
629  static value_type access(const const_col_iterator &itcol, size_type j)
630  { return col(itcol)[j]; }
631  };
632 
633  template <typename T, typename IND_TYPE, int shift>
634  std::ostream &operator <<
635  (std::ostream &o, const csc_matrix<T, IND_TYPE, shift>& m)
636  { gmm::write(o,m); return o; }
637 
638  template <typename T, typename IND_TYPE, int shift>
639  inline void copy(const identity_matrix &, csc_matrix<T, IND_TYPE, shift>& M)
640  { M.init_with_identity(mat_nrows(M)); }
641 
642  template <typename Matrix, typename T, typename IND_TYPE, int shift>
643  inline void copy(const Matrix &A, csc_matrix<T, IND_TYPE, shift>& M)
644  { M.init_with(A); }
645 
646  /* ******************************************************************** */
647  /* */
648  /* Read only compressed sparse row matrix */
649  /* */
650  /* ******************************************************************** */
651 
652  template <typename T, typename IND_TYPE = unsigned int, int shift = 0>
653  struct csr_matrix {
654 
655  std::vector<T> pr; // values.
656  std::vector<IND_TYPE> ir; // col indices.
657  std::vector<IND_TYPE> jc; // row repartition on pr and ir.
658  size_type nc, nr;
659 
660  typedef T value_type;
661  typedef T& access_type;
662 
663 
664  template <typename Matrix> void init_with_good_format(const Matrix &B);
665  void init_with(const row_matrix<wsvector<T> > &B)
666  { init_with_good_format(B); }
667  void init_with(const row_matrix<rsvector<T> > &B)
668  { init_with_good_format(B); }
669  template <typename PT1, typename PT2, typename PT3, int cshift>
670  void init_with(const csr_matrix_ref<PT1,PT2,PT3,cshift>& B)
671  { init_with_good_format(B); }
672  template <typename U, int cshift>
673  void init_with(const csr_matrix<U, IND_TYPE, cshift>& B)
674  { init_with_good_format(B); }
675 
676  template <typename Matrix> void init_with(const Matrix &A);
677  void init_with_identity(size_type n);
678 
679  csr_matrix() : nc(0), nr(0) {}
680  csr_matrix(size_type nnr, size_type nnc);
681 
682  size_type nrows() const { return nr; }
683  size_type ncols() const { return nc; }
684  void swap(csr_matrix<T, IND_TYPE, shift> &m) {
685  std::swap(pr, m.pr);
686  std::swap(ir,m.ir); std::swap(jc, m.jc);
687  std::swap(nc, m.nc); std::swap(nr,m.nr);
688  }
689 
690  value_type operator()(size_type i, size_type j) const
691  { return mat_row(*this, i)[j]; }
692  };
693 
694  template <typename T, typename IND_TYPE, int shift> template <typename Matrix>
695  void csr_matrix<T, IND_TYPE, shift>::init_with_good_format(const Matrix &B) {
696  typedef typename linalg_traits<Matrix>::const_sub_row_type row_type;
697  nc = mat_ncols(B); nr = mat_nrows(B);
698  jc.resize(nr+1);
699  jc[0] = shift;
700  for (size_type j = 0; j < nr; ++j) {
701  jc[j+1] = IND_TYPE(jc[j] + nnz(mat_const_row(B, j)));
702  }
703  pr.resize(jc[nr]);
704  ir.resize(jc[nr]);
705  for (size_type j = 0; j < nr; ++j) {
706  row_type row = mat_const_row(B, j);
707  typename linalg_traits<typename org_type<row_type>::t>::const_iterator
708  it = vect_const_begin(row), ite = vect_const_end(row);
709  for (size_type k = 0; it != ite; ++it, ++k) {
710  pr[jc[j]-shift+k] = *it;
711  ir[jc[j]-shift+k] = IND_TYPE(it.index()+shift);
712  }
713  }
714  }
715 
716  template <typename T, typename IND_TYPE, int shift> template <typename Matrix>
717  void csr_matrix<T, IND_TYPE, shift>::init_with(const Matrix &A) {
718  row_matrix<wsvector<T> > B(mat_nrows(A), mat_ncols(A));
719  copy(A, B);
720  init_with_good_format(B);
721  }
722 
723  template <typename T, typename IND_TYPE, int shift>
724  void csr_matrix<T, IND_TYPE, shift>::init_with_identity(size_type n) {
725  nc = nr = n;
726  pr.resize(nr); ir.resize(nr); jc.resize(nr+1);
727  for (size_type j = 0; j < nr; ++j)
728  { ir[j] = jc[j] = shift + j; pr[j] = T(1); }
729  jc[nr] = shift + nr;
730  }
731 
732  template <typename T, typename IND_TYPE, int shift>
733  csr_matrix<T, IND_TYPE, shift>::csr_matrix(size_type nnr, size_type nnc)
734  : nc(nnc), nr(nnr) {
735  pr.resize(1); ir.resize(1); jc.resize(nr+1);
736  for (size_type j = 0; j < nr; ++j) jc[j] = shift;
737  jc[nr] = shift;
738  }
739 
740 
741  template <typename T, typename IND_TYPE, int shift>
742  struct linalg_traits<csr_matrix<T, IND_TYPE, shift> > {
743  typedef csr_matrix<T, IND_TYPE, shift> this_type;
744  typedef linalg_const is_reference;
745  typedef abstract_matrix linalg_type;
746  typedef T value_type;
747  typedef T origin_type;
748  typedef T reference;
749  typedef abstract_sparse storage_type;
750  typedef abstract_null_type sub_col_type;
751  typedef abstract_null_type const_sub_col_type;
752  typedef abstract_null_type col_iterator;
753  typedef abstract_null_type const_col_iterator;
754  typedef abstract_null_type sub_row_type;
755  typedef cs_vector_ref<const T *, const IND_TYPE *, shift>
756  const_sub_row_type;
757  typedef sparse_compressed_iterator<const T *, const IND_TYPE *,
758  const IND_TYPE *, shift>
759  const_row_iterator;
760  typedef abstract_null_type row_iterator;
761  typedef row_major sub_orientation;
762  typedef linalg_true index_sorted;
763  static size_type nrows(const this_type &m) { return m.nrows(); }
764  static size_type ncols(const this_type &m) { return m.ncols(); }
765  static const_row_iterator row_begin(const this_type &m)
766  { return const_row_iterator(&m.pr[0], &m.ir[0], &m.jc[0], m.nc, &m.pr[0]); }
767  static const_row_iterator row_end(const this_type &m)
768  { return const_row_iterator(&m.pr[0], &m.ir[0], &m.jc[0] + m.nr, m.nc, &m.pr[0]); }
769  static const_sub_row_type row(const const_row_iterator &it) {
770  return const_sub_row_type(it.pr + *(it.jc) - shift,
771  it.ir + *(it.jc) - shift,
772  *(it.jc + 1) - *(it.jc), it.n);
773  }
774  static const origin_type* origin(const this_type &m) { return &m.pr[0]; }
775  static void do_clear(this_type &m) { m.do_clear(); }
776  static value_type access(const const_row_iterator &itrow, size_type j)
777  { return row(itrow)[j]; }
778  };
779 
780  template <typename T, typename IND_TYPE, int shift>
781  std::ostream &operator <<
782  (std::ostream &o, const csr_matrix<T, IND_TYPE, shift>& m)
783  { gmm::write(o,m); return o; }
784 
785  template <typename T, typename IND_TYPE, int shift>
786  inline void copy(const identity_matrix &, csr_matrix<T, IND_TYPE, shift>& M)
787  { M.init_with_identity(mat_nrows(M)); }
788 
789  template <typename Matrix, typename T, typename IND_TYPE, int shift>
790  inline void copy(const Matrix &A, csr_matrix<T, IND_TYPE, shift>& M)
791  { M.init_with(A); }
792 
793  /* ******************************************************************** */
794  /* */
795  /* Block matrix */
796  /* */
797  /* ******************************************************************** */
798 
799  template <typename MAT> class block_matrix {
800  protected :
801  std::vector<MAT> blocks;
802  size_type nrowblocks_;
803  size_type ncolblocks_;
804  std::vector<sub_interval> introw, intcol;
805 
806  public :
807  typedef typename linalg_traits<MAT>::value_type value_type;
808  typedef typename linalg_traits<MAT>::reference reference;
809 
810  size_type nrows() const { return introw[nrowblocks_-1].max; }
811  size_type ncols() const { return intcol[ncolblocks_-1].max; }
812  size_type nrowblocks() const { return nrowblocks_; }
813  size_type ncolblocks() const { return ncolblocks_; }
814  const sub_interval &subrowinterval(size_type i) const { return introw[i]; }
815  const sub_interval &subcolinterval(size_type i) const { return intcol[i]; }
816  const MAT &block(size_type i, size_type j) const
817  { return blocks[j*ncolblocks_+i]; }
818  MAT &block(size_type i, size_type j)
819  { return blocks[j*ncolblocks_+i]; }
820  void do_clear();
821  // to be done : read and write access to a component
822  value_type operator() (size_type i, size_type j) const {
823  size_type k, l;
824  for (k = 0; k < nrowblocks_; ++k)
825  if (i >= introw[k].min && i < introw[k].max) break;
826  for (l = 0; l < nrowblocks_; ++l)
827  if (j >= introw[l].min && j < introw[l].max) break;
828  return (block(k, l))(i - introw[k].min, j - introw[l].min);
829  }
830  reference operator() (size_type i, size_type j) {
831  size_type k, l;
832  for (k = 0; k < nrowblocks_; ++k)
833  if (i >= introw[k].min && i < introw[k].max) break;
834  for (l = 0; l < nrowblocks_; ++l)
835  if (j >= introw[l].min && j < introw[l].max) break;
836  return (block(k, l))(i - introw[k].min, j - introw[l].min);
837  }
838 
839  template <typename CONT> void resize(const CONT &c1, const CONT &c2);
840  template <typename CONT> block_matrix(const CONT &c1, const CONT &c2)
841  { resize(c1, c2); }
842  block_matrix() {}
843 
844  };
845 
846  template <typename MAT> struct linalg_traits<block_matrix<MAT> > {
847  typedef block_matrix<MAT> this_type;
848  typedef linalg_false is_reference;
849  typedef abstract_matrix linalg_type;
850  typedef this_type origin_type;
851  typedef typename linalg_traits<MAT>::value_type value_type;
852  typedef typename linalg_traits<MAT>::reference reference;
853  typedef typename linalg_traits<MAT>::storage_type storage_type;
854  typedef abstract_null_type sub_row_type; // to be done ...
855  typedef abstract_null_type const_sub_row_type; // to be done ...
856  typedef abstract_null_type row_iterator; // to be done ...
857  typedef abstract_null_type const_row_iterator; // to be done ...
858  typedef abstract_null_type sub_col_type; // to be done ...
859  typedef abstract_null_type const_sub_col_type; // to be done ...
860  typedef abstract_null_type col_iterator; // to be done ...
861  typedef abstract_null_type const_col_iterator; // to be done ...
862  typedef abstract_null_type sub_orientation; // to be done ...
863  typedef linalg_true index_sorted;
864  static size_type nrows(const this_type &m) { return m.nrows(); }
865  static size_type ncols(const this_type &m) { return m.ncols(); }
866  static origin_type* origin(this_type &m) { return &m; }
867  static const origin_type* origin(const this_type &m) { return &m; }
868  static void do_clear(this_type &m) { m.do_clear(); }
869  // access to be done ...
870  static void resize(this_type &, size_type , size_type)
871  { GMM_ASSERT1(false, "Sorry, to be done"); }
872  static void reshape(this_type &, size_type , size_type)
873  { GMM_ASSERT1(false, "Sorry, to be done"); }
874  };
875 
876  template <typename MAT>
877  void block_matrix<MAT>::do_clear() {
878  for (size_type j = 0; j < ncolblocks_; ++j)
879  for (size_type i = 0; i < nrowblocks_; ++i)
880  clear(block(i,j));
881  }
882 
883  template <typename MAT> template <typename CONT>
884  void block_matrix<MAT>::resize(const CONT &c1, const CONT &c2) {
885  nrowblocks_ = c1.size(); ncolblocks_ = c2.size();
886  blocks.resize(nrowblocks_ * ncolblocks_);
887  intcol.resize(ncolblocks_);
888  introw.resize(nrowblocks_);
889  for (size_type j = 0, l = 0; j < ncolblocks_; ++j) {
890  intcol[j] = sub_interval(l, c2[j]); l += c2[j];
891  for (size_type i = 0, k = 0; i < nrowblocks_; ++i) {
892  if (j == 0) { introw[i] = sub_interval(k, c1[i]); k += c1[i]; }
893  block(i, j) = MAT(c1[i], c2[j]);
894  }
895  }
896  }
897 
898  template <typename M1, typename M2>
899  void copy(const block_matrix<M1> &m1, M2 &m2) {
900  for (size_type j = 0; j < m1.ncolblocks(); ++j)
901  for (size_type i = 0; i < m1.nrowblocks(); ++i)
902  copy(m1.block(i,j), sub_matrix(m2, m1.subrowinterval(i),
903  m1.subcolinterval(j)));
904  }
905 
906  template <typename M1, typename M2>
907  void copy(const block_matrix<M1> &m1, const M2 &m2)
908  { copy(m1, linalg_const_cast(m2)); }
909 
910 
911  template <typename MAT, typename V1, typename V2>
912  void mult(const block_matrix<MAT> &m, const V1 &v1, V2 &v2) {
913  clear(v2);
914  typename sub_vector_type<V2 *, sub_interval>::vector_type sv;
915  for (size_type i = 0; i < m.nrowblocks() ; ++i)
916  for (size_type j = 0; j < m.ncolblocks() ; ++j) {
917  sv = sub_vector(v2, m.subrowinterval(i));
918  mult(m.block(i,j),
919  sub_vector(v1, m.subcolinterval(j)), sv, sv);
920  }
921  }
922 
923  template <typename MAT, typename V1, typename V2, typename V3>
924  void mult(const block_matrix<MAT> &m, const V1 &v1, const V2 &v2, V3 &v3) {
925  typename sub_vector_type<V3 *, sub_interval>::vector_type sv;
926  for (size_type i = 0; i < m.nrowblocks() ; ++i)
927  for (size_type j = 0; j < m.ncolblocks() ; ++j) {
928  sv = sub_vector(v3, m.subrowinterval(i));
929  if (j == 0)
930  mult(m.block(i,j),
931  sub_vector(v1, m.subcolinterval(j)),
932  sub_vector(v2, m.subrowinterval(i)), sv);
933  else
934  mult(m.block(i,j),
935  sub_vector(v1, m.subcolinterval(j)), sv, sv);
936  }
937 
938  }
939 
940  template <typename MAT, typename V1, typename V2>
941  void mult(const block_matrix<MAT> &m, const V1 &v1, const V2 &v2)
942  { mult(m, v1, linalg_const_cast(v2)); }
943 
944  template <typename MAT, typename V1, typename V2, typename V3>
945  void mult(const block_matrix<MAT> &m, const V1 &v1, const V2 &v2,
946  const V3 &v3)
947  { mult_const(m, v1, v2, linalg_const_cast(v3)); }
948 
949 }
950  /* ******************************************************************** */
951  /* */
952  /* Distributed matrices */
953  /* */
954  /* ******************************************************************** */
955 
956 #ifdef GMM_USES_MPI
957 # include <mpi.h>
958 
959 namespace gmm {
960 
961 
962 
963  template <typename T> inline MPI_Datatype mpi_type(T)
964  { GMM_ASSERT1(false, "Sorry unsupported type"); return MPI_FLOAT; }
965  inline MPI_Datatype mpi_type(double) { return MPI_DOUBLE; }
966  inline MPI_Datatype mpi_type(float) { return MPI_FLOAT; }
967  inline MPI_Datatype mpi_type(long double) { return MPI_LONG_DOUBLE; }
968 #ifndef LAM_MPI
969  inline MPI_Datatype mpi_type(std::complex<float>) { return MPI_COMPLEX; }
970  inline MPI_Datatype mpi_type(std::complex<double>) { return MPI_DOUBLE_COMPLEX; }
971 #endif
972  inline MPI_Datatype mpi_type(int) { return MPI_INT; }
973  inline MPI_Datatype mpi_type(unsigned int) { return MPI_UNSIGNED; }
974  inline MPI_Datatype mpi_type(long) { return MPI_LONG; }
975  inline MPI_Datatype mpi_type(unsigned long) { return MPI_UNSIGNED_LONG; }
976 
977  template <typename MAT> struct mpi_distributed_matrix {
978  MAT M;
979 
980  mpi_distributed_matrix(size_type n, size_type m) : M(n, m) {}
981  mpi_distributed_matrix() {}
982 
983  const MAT &local_matrix() const { return M; }
984  MAT &local_matrix() { return M; }
985  };
986 
987  template <typename MAT> inline MAT &eff_matrix(MAT &m) { return m; }
988  template <typename MAT> inline
989  const MAT &eff_matrix(const MAT &m) { return m; }
990  template <typename MAT> inline
991  MAT &eff_matrix(mpi_distributed_matrix<MAT> &m) { return m.M; }
992  template <typename MAT> inline
993  const MAT &eff_matrix(const mpi_distributed_matrix<MAT> &m) { return m.M; }
994 
995 
996  template <typename MAT1, typename MAT2>
997  inline void copy(const mpi_distributed_matrix<MAT1> &m1,
998  mpi_distributed_matrix<MAT2> &m2)
999  { copy(eff_matrix(m1), eff_matrix(m2)); }
1000  template <typename MAT1, typename MAT2>
1001  inline void copy(const mpi_distributed_matrix<MAT1> &m1,
1002  const mpi_distributed_matrix<MAT2> &m2)
1003  { copy(m1.M, m2.M); }
1004 
1005  template <typename MAT1, typename MAT2>
1006  inline void copy(const mpi_distributed_matrix<MAT1> &m1, MAT2 &m2)
1007  { copy(m1.M, m2); }
1008  template <typename MAT1, typename MAT2>
1009  inline void copy(const mpi_distributed_matrix<MAT1> &m1, const MAT2 &m2)
1010  { copy(m1.M, m2); }
1011 
1012 
1013  template <typename MATSP, typename V1, typename V2> inline
1014  typename strongest_value_type3<V1,V2,MATSP>::value_type
1015  vect_sp(const mpi_distributed_matrix<MATSP> &ps, const V1 &v1,
1016  const V2 &v2) {
1017  typedef typename strongest_value_type3<V1,V2,MATSP>::value_type T;
1018  T res = vect_sp(ps.M, v1, v2), rest;
1019  MPI_Allreduce(&res, &rest, 1, mpi_type(T()), MPI_SUM,MPI_COMM_WORLD);
1020  return rest;
1021  }
1022 
1023  template <typename MAT, typename V1, typename V2>
1024  inline void mult_add(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1025  V2 &v2) {
1026  typedef typename linalg_traits<V2>::value_type T;
1027  std::vector<T> v3(vect_size(v2)), v4(vect_size(v2));
1028  static double tmult_tot = 0.0;
1029  static double tmult_tot2 = 0.0;
1030  double t_ref = MPI_Wtime();
1031  gmm::mult(m.M, v1, v3);
1032  if (is_sparse(v2)) GMM_WARNING2("Using a plain temporary, here.");
1033  double t_ref2 = MPI_Wtime();
1034  MPI_Allreduce(&(v3[0]), &(v4[0]),gmm::vect_size(v2), mpi_type(T()),
1035  MPI_SUM,MPI_COMM_WORLD);
1036  tmult_tot2 = MPI_Wtime()-t_ref2;
1037  cout << "reduce mult mpi = " << tmult_tot2 << endl;
1038  gmm::add(v4, v2);
1039  tmult_tot = MPI_Wtime()-t_ref;
1040  cout << "tmult mpi = " << tmult_tot << endl;
1041  }
1042 
1043  template <typename MAT, typename V1, typename V2>
1044  void mult_add(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1045  const V2 &v2_)
1046  { mult_add(m, v1, const_cast<V2 &>(v2_)); }
1047 
1048  template <typename MAT, typename V1, typename V2>
1049  inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1050  const V2 &v2_)
1051  { V2 &v2 = const_cast<V2 &>(v2_); clear(v2); mult_add(m, v1, v2); }
1052 
1053  template <typename MAT, typename V1, typename V2>
1054  inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1055  V2 &v2)
1056  { clear(v2); mult_add(m, v1, v2); }
1057 
1058  template <typename MAT, typename V1, typename V2, typename V3>
1059  inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1060  const V2 &v2, const V3 &v3_)
1061  { V3 &v3 = const_cast<V3 &>(v3_); gmm::copy(v2, v3); mult_add(m, v1, v3); }
1062 
1063  template <typename MAT, typename V1, typename V2, typename V3>
1064  inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1065  const V2 &v2, V3 &v3)
1066  { gmm::copy(v2, v3); mult_add(m, v1, v3); }
1067 
1068 
1069  template <typename MAT> inline
1070  size_type mat_nrows(const mpi_distributed_matrix<MAT> &M)
1071  { return mat_nrows(M.M); }
1072  template <typename MAT> inline
1073  size_type mat_ncols(const mpi_distributed_matrix<MAT> &M)
1074  { return mat_nrows(M.M); }
1075  template <typename MAT> inline
1076  void resize(mpi_distributed_matrix<MAT> &M, size_type m, size_type n)
1077  { resize(M.M, m, n); }
1078  template <typename MAT> inline void clear(mpi_distributed_matrix<MAT> &M)
1079  { clear(M.M); }
1080 
1081 
1082  // For compute reduced system
1083  template <typename MAT1, typename MAT2> inline
1084  void mult(const MAT1 &M1, const mpi_distributed_matrix<MAT2> &M2,
1085  mpi_distributed_matrix<MAT2> &M3)
1086  { mult(M1, M2.M, M3.M); }
1087  template <typename MAT1, typename MAT2> inline
1088  void mult(const mpi_distributed_matrix<MAT2> &M2,
1089  const MAT1 &M1, mpi_distributed_matrix<MAT2> &M3)
1090  { mult(M2.M, M1, M3.M); }
1091  template <typename MAT1, typename MAT2, typename MAT3> inline
1092  void mult(const MAT1 &M1, const mpi_distributed_matrix<MAT2> &M2,
1093  MAT3 &M3)
1094  { mult(M1, M2.M, M3); }
1095  template <typename MAT1, typename MAT2, typename MAT3> inline
1096  void mult(const MAT1 &M1, const mpi_distributed_matrix<MAT2> &M2,
1097  const MAT3 &M3)
1098  { mult(M1, M2.M, M3); }
1099 
1100  template <typename M, typename SUBI1, typename SUBI2>
1101  struct sub_matrix_type<const mpi_distributed_matrix<M> *, SUBI1, SUBI2>
1102  { typedef abstract_null_type matrix_type; };
1103 
1104  template <typename M, typename SUBI1, typename SUBI2>
1105  struct sub_matrix_type<mpi_distributed_matrix<M> *, SUBI1, SUBI2>
1106  { typedef abstract_null_type matrix_type; };
1107 
1108  template <typename M, typename SUBI1, typename SUBI2> inline
1109  typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI2>
1110  ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI2>::matrix_type,
1111  M *>::return_type
1112  sub_matrix(mpi_distributed_matrix<M> &m, const SUBI1 &si1, const SUBI2 &si2)
1113  { return sub_matrix(m.M, si1, si2); }
1114 
1115  template <typename MAT, typename SUBI1, typename SUBI2> inline
1116  typename select_return<typename sub_matrix_type<const MAT *, SUBI1, SUBI2>
1117  ::matrix_type, typename sub_matrix_type<MAT *, SUBI1, SUBI2>::matrix_type,
1118  const MAT *>::return_type
1119  sub_matrix(const mpi_distributed_matrix<MAT> &m, const SUBI1 &si1,
1120  const SUBI2 &si2)
1121  { return sub_matrix(m.M, si1, si2); }
1122 
1123  template <typename M, typename SUBI1> inline
1124  typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI1>
1125  ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type,
1126  M *>::return_type
1127  sub_matrix(mpi_distributed_matrix<M> &m, const SUBI1 &si1)
1128  { return sub_matrix(m.M, si1, si1); }
1129 
1130  template <typename M, typename SUBI1> inline
1131  typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI1>
1132  ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type,
1133  const M *>::return_type
1134  sub_matrix(const mpi_distributed_matrix<M> &m, const SUBI1 &si1)
1135  { return sub_matrix(m.M, si1, si1); }
1136 
1137 
1138  template <typename L> struct transposed_return<const mpi_distributed_matrix<L> *>
1139  { typedef abstract_null_type return_type; };
1140  template <typename L> struct transposed_return<mpi_distributed_matrix<L> *>
1141  { typedef abstract_null_type return_type; };
1142 
1143  template <typename L> inline typename transposed_return<const L *>::return_type
1144  transposed(const mpi_distributed_matrix<L> &l)
1145  { return transposed(l.M); }
1146 
1147  template <typename L> inline typename transposed_return<L *>::return_type
1148  transposed(mpi_distributed_matrix<L> &l)
1149  { return transposed(l.M); }
1150 
1151 
1152  template <typename MAT>
1153  struct linalg_traits<mpi_distributed_matrix<MAT> > {
1154  typedef mpi_distributed_matrix<MAT> this_type;
1155  typedef MAT origin_type;
1156  typedef linalg_false is_reference;
1157  typedef abstract_matrix linalg_type;
1158  typedef typename linalg_traits<MAT>::value_type value_type;
1159  typedef typename linalg_traits<MAT>::reference reference;
1160  typedef typename linalg_traits<MAT>::storage_type storage_type;
1161  typedef abstract_null_type sub_row_type;
1162  typedef abstract_null_type const_sub_row_type;
1163  typedef abstract_null_type row_iterator;
1164  typedef abstract_null_type const_row_iterator;
1165  typedef abstract_null_type sub_col_type;
1166  typedef abstract_null_type const_sub_col_type;
1167  typedef abstract_null_type col_iterator;
1168  typedef abstract_null_type const_col_iterator;
1169  typedef abstract_null_type sub_orientation;
1170  typedef abstract_null_type index_sorted;
1171  static size_type nrows(const this_type &m) { return nrows(m.M); }
1172  static size_type ncols(const this_type &m) { return ncols(m.M); }
1173  static void do_clear(this_type &m) { clear(m.M); }
1174  };
1175 
1176 }
1177 
1178 
1179 #endif // GMM_USES_MPI
1180 
1181 namespace std {
1182  template <typename V>
1183  void swap(gmm::row_matrix<V> &m1, gmm::row_matrix<V> &m2)
1184  { m1.swap(m2); }
1185  template <typename V>
1186  void swap(gmm::col_matrix<V> &m1, gmm::col_matrix<V> &m2)
1187  { m1.swap(m2); }
1188  template <typename T>
1189  void swap(gmm::dense_matrix<T> &m1, gmm::dense_matrix<T> &m2)
1190  { m1.swap(m2); }
1191  template <typename T, typename IND_TYPE, int shift> void
1192  swap(gmm::csc_matrix<T, IND_TYPE, shift> &m1, gmm::csc_matrix<T, IND_TYPE, shift> &m2)
1193  { m1.swap(m2); }
1194  template <typename T, typename IND_TYPE, int shift> void
1195  swap(gmm::csr_matrix<T, IND_TYPE, shift> &m1, gmm::csr_matrix<T, IND_TYPE, shift> &m2)
1196  { m1.swap(m2); }
1197 }
1198 
1199 
1200 
1201 
1202 #endif /* GMM_MATRIX_H__ */
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_vector.h:995
gmm::nnz
size_type nnz(const L &l)
count the number of non-zero entries of a vector or matrix.
Definition: gmm_blas.h:68
gmm::mult_add
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1790
gmm::reshape
void reshape(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:250
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
gmm::fill
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*‍/
Definition: gmm_blas.h:103
gmm::vect_hp
strongest_value_type< V1, V2 >::value_type vect_hp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:510
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:209
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:263
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1275
gmm_sub_matrix.h
Generic sub-matrices.
gmm_sub_vector.h
Generic sub-vectors.
gmm_transposed.h
Generic transposed matrices.
gmm_vector.h
Declaration of the vector types (gmm::rsvector, gmm::wsvector, gmm::slvector ,..)
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48