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