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