GetFEM++  5.3
gmm_matrix.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-2017 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, int shift = 0>
501  struct csc_matrix {
502  typedef unsigned int IND_TYPE;
503 
504  std::vector<T> pr;
505  std::vector<IND_TYPE> ir;
506  std::vector<IND_TYPE> jc;
507  size_type nc, nr;
508 
509  typedef T value_type;
510  typedef T& access_type;
511 
512  template <typename Matrix> void init_with_good_format(const Matrix &B);
513  template <typename Matrix> void init_with(const Matrix &A);
514  void init_with(const col_matrix<gmm::rsvector<T> > &B)
515  { init_with_good_format(B); }
516  void init_with(const col_matrix<wsvector<T> > &B)
517  { init_with_good_format(B); }
518  template <typename PT1, typename PT2, typename PT3, int cshift>
519  void init_with(const csc_matrix_ref<PT1,PT2,PT3,cshift>& B)
520  { init_with_good_format(B); }
521  template <typename U, int cshift>
522  void init_with(const csc_matrix<U, cshift>& B)
523  { init_with_good_format(B); }
524 
525  void init_with_identity(size_type n);
526 
527  csc_matrix(void) : nc(0), nr(0) {}
528  csc_matrix(size_type nnr, size_type nnc);
529 
530  size_type nrows(void) const { return nr; }
531  size_type ncols(void) const { return nc; }
532  void swap(csc_matrix<T, shift> &m) {
533  std::swap(pr, m.pr);
534  std::swap(ir, m.ir); std::swap(jc, m.jc);
535  std::swap(nc, m.nc); std::swap(nr, m.nr);
536  }
537  value_type operator()(size_type i, size_type j) const
538  { return mat_col(*this, j)[i]; }
539  };
540 
541  template <typename T, int shift> template<typename Matrix>
542  void csc_matrix<T, shift>::init_with_good_format(const Matrix &B) {
543  typedef typename linalg_traits<Matrix>::const_sub_col_type col_type;
544  nc = mat_ncols(B); nr = mat_nrows(B);
545  jc.resize(nc+1);
546  jc[0] = shift;
547  for (size_type j = 0; j < nc; ++j) {
548  jc[j+1] = IND_TYPE(jc[j] + nnz(mat_const_col(B, j)));
549  }
550  pr.resize(jc[nc]);
551  ir.resize(jc[nc]);
552  for (size_type j = 0; j < nc; ++j) {
553  col_type col = mat_const_col(B, j);
554  typename linalg_traits<typename org_type<col_type>::t>::const_iterator
555  it = vect_const_begin(col), ite = vect_const_end(col);
556  for (size_type k = 0; it != ite; ++it, ++k) {
557  pr[jc[j]-shift+k] = *it;
558  ir[jc[j]-shift+k] = IND_TYPE(it.index() + shift);
559  }
560  }
561  }
562 
563  template <typename T, int shift> template <typename Matrix>
564  void csc_matrix<T, 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, int shift>
571  void csc_matrix<T, 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, int shift>
580  csc_matrix<T, 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, int shift>
587  struct linalg_traits<csc_matrix<T, shift> > {
588  typedef csc_matrix<T, shift> this_type;
589  typedef typename this_type::IND_TYPE IND_TYPE;
590  typedef linalg_const is_reference;
591  typedef abstract_matrix linalg_type;
592  typedef T value_type;
593  typedef T origin_type;
594  typedef T reference;
595  typedef abstract_sparse storage_type;
596  typedef abstract_null_type sub_row_type;
597  typedef abstract_null_type const_sub_row_type;
598  typedef abstract_null_type row_iterator;
599  typedef abstract_null_type const_row_iterator;
600  typedef abstract_null_type sub_col_type;
601  typedef cs_vector_ref<const T *, const IND_TYPE *, shift>
602  const_sub_col_type;
603  typedef sparse_compressed_iterator<const T *, const IND_TYPE *,
604  const IND_TYPE *, shift>
605  const_col_iterator;
606  typedef abstract_null_type col_iterator;
607  typedef col_major sub_orientation;
608  typedef linalg_true index_sorted;
609  static size_type nrows(const this_type &m) { return m.nrows(); }
610  static size_type ncols(const this_type &m) { return m.ncols(); }
611  static const_col_iterator col_begin(const this_type &m)
612  { return const_col_iterator(&m.pr[0],&m.ir[0],&m.jc[0], m.nr, &m.pr[0]); }
613  static const_col_iterator col_end(const this_type &m) {
614  return const_col_iterator(&m.pr[0],&m.ir[0],&m.jc[0]+m.nc,
615  m.nr,&m.pr[0]);
616  }
617  static const_sub_col_type col(const const_col_iterator &it) {
618  return const_sub_col_type(it.pr + *(it.jc) - shift,
619  it.ir + *(it.jc) - shift,
620  *(it.jc + 1) - *(it.jc), it.n);
621  }
622  static const origin_type* origin(const this_type &m) { return &m.pr[0]; }
623  static void do_clear(this_type &m) { m.do_clear(); }
624  static value_type access(const const_col_iterator &itcol, size_type j)
625  { return col(itcol)[j]; }
626  };
627 
628  template <typename T, int shift>
629  std::ostream &operator <<
630  (std::ostream &o, const csc_matrix<T, shift>& m)
631  { gmm::write(o,m); return o; }
632 
633  template <typename T, int shift>
634  inline void copy(const identity_matrix &, csc_matrix<T, shift>& M)
635  { M.init_with_identity(mat_nrows(M)); }
636 
637  template <typename Matrix, typename T, int shift>
638  inline void copy(const Matrix &A, csc_matrix<T, shift>& M)
639  { M.init_with(A); }
640 
641  /* ******************************************************************** */
642  /* */
643  /* Read only compressed sparse row matrix */
644  /* */
645  /* ******************************************************************** */
646 
647  template <typename T, int shift = 0>
648  struct csr_matrix {
649 
650  typedef unsigned int IND_TYPE;
651 
652  std::vector<T> pr; // values.
653  std::vector<IND_TYPE> ir; // col indices.
654  std::vector<IND_TYPE> jc; // row repartition on pr and ir.
655  size_type nc, nr;
656 
657  typedef T value_type;
658  typedef T& access_type;
659 
660 
661  template <typename Matrix> void init_with_good_format(const Matrix &B);
662  void init_with(const row_matrix<wsvector<T> > &B)
663  { init_with_good_format(B); }
664  void init_with(const row_matrix<rsvector<T> > &B)
665  { init_with_good_format(B); }
666  template <typename PT1, typename PT2, typename PT3, int cshift>
667  void init_with(const csr_matrix_ref<PT1,PT2,PT3,cshift>& B)
668  { init_with_good_format(B); }
669  template <typename U, int cshift>
670  void init_with(const csr_matrix<U, cshift>& B)
671  { init_with_good_format(B); }
672 
673  template <typename Matrix> void init_with(const Matrix &A);
674  void init_with_identity(size_type n);
675 
676  csr_matrix(void) : nc(0), nr(0) {}
677  csr_matrix(size_type nnr, size_type nnc);
678 
679  size_type nrows(void) const { return nr; }
680  size_type ncols(void) const { return nc; }
681  void swap(csr_matrix<T, shift> &m) {
682  std::swap(pr, m.pr);
683  std::swap(ir,m.ir); std::swap(jc, m.jc);
684  std::swap(nc, m.nc); std::swap(nr,m.nr);
685  }
686 
687  value_type operator()(size_type i, size_type j) const
688  { return mat_row(*this, i)[j]; }
689  };
690 
691  template <typename T, int shift> template <typename Matrix>
692  void csr_matrix<T, shift>::init_with_good_format(const Matrix &B) {
693  typedef typename linalg_traits<Matrix>::const_sub_row_type row_type;
694  nc = mat_ncols(B); nr = mat_nrows(B);
695  jc.resize(nr+1);
696  jc[0] = shift;
697  for (size_type j = 0; j < nr; ++j) {
698  jc[j+1] = IND_TYPE(jc[j] + nnz(mat_const_row(B, j)));
699  }
700  pr.resize(jc[nr]);
701  ir.resize(jc[nr]);
702  for (size_type j = 0; j < nr; ++j) {
703  row_type row = mat_const_row(B, j);
704  typename linalg_traits<typename org_type<row_type>::t>::const_iterator
705  it = vect_const_begin(row), ite = vect_const_end(row);
706  for (size_type k = 0; it != ite; ++it, ++k) {
707  pr[jc[j]-shift+k] = *it;
708  ir[jc[j]-shift+k] = IND_TYPE(it.index()+shift);
709  }
710  }
711  }
712 
713  template <typename T, int shift> template <typename Matrix>
714  void csr_matrix<T, shift>::init_with(const Matrix &A) {
715  row_matrix<wsvector<T> > B(mat_nrows(A), mat_ncols(A));
716  copy(A, B);
717  init_with_good_format(B);
718  }
719 
720  template <typename T, int shift>
721  void csr_matrix<T, shift>::init_with_identity(size_type n) {
722  nc = nr = n;
723  pr.resize(nr); ir.resize(nr); jc.resize(nr+1);
724  for (size_type j = 0; j < nr; ++j)
725  { ir[j] = jc[j] = shift + j; pr[j] = T(1); }
726  jc[nr] = shift + nr;
727  }
728 
729  template <typename T, int shift>
730  csr_matrix<T, shift>::csr_matrix(size_type nnr, size_type nnc)
731  : nc(nnc), nr(nnr) {
732  pr.resize(1); ir.resize(1); jc.resize(nr+1);
733  for (size_type j = 0; j < nr; ++j) jc[j] = shift;
734  jc[nr] = shift;
735  }
736 
737 
738  template <typename T, int shift>
739  struct linalg_traits<csr_matrix<T, shift> > {
740  typedef csr_matrix<T, shift> this_type;
741  typedef typename this_type::IND_TYPE IND_TYPE;
742  typedef linalg_const is_reference;
743  typedef abstract_matrix linalg_type;
744  typedef T value_type;
745  typedef T origin_type;
746  typedef T reference;
747  typedef abstract_sparse storage_type;
748  typedef abstract_null_type sub_col_type;
749  typedef abstract_null_type const_sub_col_type;
750  typedef abstract_null_type col_iterator;
751  typedef abstract_null_type const_col_iterator;
752  typedef abstract_null_type sub_row_type;
753  typedef cs_vector_ref<const T *, const IND_TYPE *, shift>
754  const_sub_row_type;
755  typedef sparse_compressed_iterator<const T *, const IND_TYPE *,
756  const IND_TYPE *, shift>
757  const_row_iterator;
758  typedef abstract_null_type row_iterator;
759  typedef row_major sub_orientation;
760  typedef linalg_true index_sorted;
761  static size_type nrows(const this_type &m) { return m.nrows(); }
762  static size_type ncols(const this_type &m) { return m.ncols(); }
763  static const_row_iterator row_begin(const this_type &m)
764  { return const_row_iterator(&m.pr[0], &m.ir[0], &m.jc[0], m.nc, &m.pr[0]); }
765  static const_row_iterator row_end(const this_type &m)
766  { return const_row_iterator(&m.pr[0], &m.ir[0], &m.jc[0] + m.nr, m.nc, &m.pr[0]); }
767  static const_sub_row_type row(const const_row_iterator &it) {
768  return const_sub_row_type(it.pr + *(it.jc) - shift,
769  it.ir + *(it.jc) - shift,
770  *(it.jc + 1) - *(it.jc), it.n);
771  }
772  static const origin_type* origin(const this_type &m) { return &m.pr[0]; }
773  static void do_clear(this_type &m) { m.do_clear(); }
774  static value_type access(const const_row_iterator &itrow, size_type j)
775  { return row(itrow)[j]; }
776  };
777 
778  template <typename T, int shift>
779  std::ostream &operator <<
780  (std::ostream &o, const csr_matrix<T, shift>& m)
781  { gmm::write(o,m); return o; }
782 
783  template <typename T, int shift>
784  inline void copy(const identity_matrix &, csr_matrix<T, shift>& M)
785  { M.init_with_identity(mat_nrows(M)); }
786 
787  template <typename Matrix, typename T, int shift>
788  inline void copy(const Matrix &A, csr_matrix<T, shift>& M)
789  { M.init_with(A); }
790 
791  /* ******************************************************************** */
792  /* */
793  /* Block matrix */
794  /* */
795  /* ******************************************************************** */
796 
797  template <typename MAT> class block_matrix {
798  protected :
799  std::vector<MAT> blocks;
800  size_type nrowblocks_;
801  size_type ncolblocks_;
802  std::vector<sub_interval> introw, intcol;
803 
804  public :
805  typedef typename linalg_traits<MAT>::value_type value_type;
806  typedef typename linalg_traits<MAT>::reference reference;
807 
808  size_type nrows(void) const { return introw[nrowblocks_-1].max; }
809  size_type ncols(void) const { return intcol[ncolblocks_-1].max; }
810  size_type nrowblocks(void) const { return nrowblocks_; }
811  size_type ncolblocks(void) const { return ncolblocks_; }
812  const sub_interval &subrowinterval(size_type i) const { return introw[i]; }
813  const sub_interval &subcolinterval(size_type i) const { return intcol[i]; }
814  const MAT &block(size_type i, size_type j) const
815  { return blocks[j*ncolblocks_+i]; }
816  MAT &block(size_type i, size_type j)
817  { return blocks[j*ncolblocks_+i]; }
818  void do_clear(void);
819  // to be done : read and write access to a component
820  value_type operator() (size_type i, size_type j) const {
821  size_type k, l;
822  for (k = 0; k < nrowblocks_; ++k)
823  if (i >= introw[k].min && i < introw[k].max) break;
824  for (l = 0; l < nrowblocks_; ++l)
825  if (j >= introw[l].min && j < introw[l].max) break;
826  return (block(k, l))(i - introw[k].min, j - introw[l].min);
827  }
828  reference operator() (size_type i, size_type j) {
829  size_type k, l;
830  for (k = 0; k < nrowblocks_; ++k)
831  if (i >= introw[k].min && i < introw[k].max) break;
832  for (l = 0; l < nrowblocks_; ++l)
833  if (j >= introw[l].min && j < introw[l].max) break;
834  return (block(k, l))(i - introw[k].min, j - introw[l].min);
835  }
836 
837  template <typename CONT> void resize(const CONT &c1, const CONT &c2);
838  template <typename CONT> block_matrix(const CONT &c1, const CONT &c2)
839  { resize(c1, c2); }
840  block_matrix(void) {}
841 
842  };
843 
844  template <typename MAT> struct linalg_traits<block_matrix<MAT> > {
845  typedef block_matrix<MAT> this_type;
846  typedef linalg_false is_reference;
847  typedef abstract_matrix linalg_type;
848  typedef this_type origin_type;
849  typedef typename linalg_traits<MAT>::value_type value_type;
850  typedef typename linalg_traits<MAT>::reference reference;
851  typedef typename linalg_traits<MAT>::storage_type storage_type;
852  typedef abstract_null_type sub_row_type; // to be done ...
853  typedef abstract_null_type const_sub_row_type; // to be done ...
854  typedef abstract_null_type row_iterator; // to be done ...
855  typedef abstract_null_type const_row_iterator; // to be done ...
856  typedef abstract_null_type sub_col_type; // to be done ...
857  typedef abstract_null_type const_sub_col_type; // to be done ...
858  typedef abstract_null_type col_iterator; // to be done ...
859  typedef abstract_null_type const_col_iterator; // to be done ...
860  typedef abstract_null_type sub_orientation; // to be done ...
861  typedef linalg_true index_sorted;
862  static size_type nrows(const this_type &m) { return m.nrows(); }
863  static size_type ncols(const this_type &m) { return m.ncols(); }
864  static origin_type* origin(this_type &m) { return &m; }
865  static const origin_type* origin(const this_type &m) { return &m; }
866  static void do_clear(this_type &m) { m.do_clear(); }
867  // access to be done ...
868  static void resize(this_type &, size_type , size_type)
869  { GMM_ASSERT1(false, "Sorry, to be done"); }
870  static void reshape(this_type &, size_type , size_type)
871  { GMM_ASSERT1(false, "Sorry, to be done"); }
872  };
873 
874  template <typename MAT> void block_matrix<MAT>::do_clear(void) {
875  for (size_type j = 0, l = 0; j < ncolblocks_; ++j)
876  for (size_type i = 0, k = 0; i < nrowblocks_; ++i)
877  clear(block(i,j));
878  }
879 
880  template <typename MAT> template <typename CONT>
881  void block_matrix<MAT>::resize(const CONT &c1, const CONT &c2) {
882  nrowblocks_ = c1.size(); ncolblocks_ = c2.size();
883  blocks.resize(nrowblocks_ * ncolblocks_);
884  intcol.resize(ncolblocks_);
885  introw.resize(nrowblocks_);
886  for (size_type j = 0, l = 0; j < ncolblocks_; ++j) {
887  intcol[j] = sub_interval(l, c2[j]); l += c2[j];
888  for (size_type i = 0, k = 0; i < nrowblocks_; ++i) {
889  if (j == 0) { introw[i] = sub_interval(k, c1[i]); k += c1[i]; }
890  block(i, j) = MAT(c1[i], c2[j]);
891  }
892  }
893  }
894 
895  template <typename M1, typename M2>
896  void copy(const block_matrix<M1> &m1, M2 &m2) {
897  for (size_type j = 0; j < m1.ncolblocks(); ++j)
898  for (size_type i = 0; i < m1.nrowblocks(); ++i)
899  copy(m1.block(i,j), sub_matrix(m2, m1.subrowinterval(i),
900  m1.subcolinterval(j)));
901  }
902 
903  template <typename M1, typename M2>
904  void copy(const block_matrix<M1> &m1, const M2 &m2)
905  { copy(m1, linalg_const_cast(m2)); }
906 
907 
908  template <typename MAT, typename V1, typename V2>
909  void mult(const block_matrix<MAT> &m, const V1 &v1, V2 &v2) {
910  clear(v2);
911  typename sub_vector_type<V2 *, sub_interval>::vector_type sv;
912  for (size_type i = 0; i < m.nrowblocks() ; ++i)
913  for (size_type j = 0; j < m.ncolblocks() ; ++j) {
914  sv = sub_vector(v2, m.subrowinterval(i));
915  mult(m.block(i,j),
916  sub_vector(v1, m.subcolinterval(j)), sv, sv);
917  }
918  }
919 
920  template <typename MAT, typename V1, typename V2, typename V3>
921  void mult(const block_matrix<MAT> &m, const V1 &v1, const V2 &v2, V3 &v3) {
922  typename sub_vector_type<V3 *, sub_interval>::vector_type sv;
923  for (size_type i = 0; i < m.nrowblocks() ; ++i)
924  for (size_type j = 0; j < m.ncolblocks() ; ++j) {
925  sv = sub_vector(v3, m.subrowinterval(i));
926  if (j == 0)
927  mult(m.block(i,j),
928  sub_vector(v1, m.subcolinterval(j)),
929  sub_vector(v2, m.subrowinterval(i)), sv);
930  else
931  mult(m.block(i,j),
932  sub_vector(v1, m.subcolinterval(j)), sv, sv);
933  }
934 
935  }
936 
937  template <typename MAT, typename V1, typename V2>
938  void mult(const block_matrix<MAT> &m, const V1 &v1, const V2 &v2)
939  { mult(m, v1, linalg_const_cast(v2)); }
940 
941  template <typename MAT, typename V1, typename V2, typename V3>
942  void mult(const block_matrix<MAT> &m, const V1 &v1, const V2 &v2,
943  const V3 &v3)
944  { mult_const(m, v1, v2, linalg_const_cast(v3)); }
945 
946 }
947  /* ******************************************************************** */
948  /* */
949  /* Distributed matrices */
950  /* */
951  /* ******************************************************************** */
952 
953 #ifdef GMM_USES_MPI
954 # include <mpi.h>
955 
956 namespace gmm {
957 
958 
959 
960  template <typename T> inline MPI_Datatype mpi_type(T)
961  { GMM_ASSERT1(false, "Sorry unsupported type"); return MPI_FLOAT; }
962  inline MPI_Datatype mpi_type(double) { return MPI_DOUBLE; }
963  inline MPI_Datatype mpi_type(float) { return MPI_FLOAT; }
964  inline MPI_Datatype mpi_type(long double) { return MPI_LONG_DOUBLE; }
965 #ifndef LAM_MPI
966  inline MPI_Datatype mpi_type(std::complex<float>) { return MPI_COMPLEX; }
967  inline MPI_Datatype mpi_type(std::complex<double>) { return MPI_DOUBLE_COMPLEX; }
968 #endif
969  inline MPI_Datatype mpi_type(int) { return MPI_INT; }
970  inline MPI_Datatype mpi_type(unsigned int) { return MPI_UNSIGNED; }
971  inline MPI_Datatype mpi_type(long) { return MPI_LONG; }
972  inline MPI_Datatype mpi_type(unsigned long) { return MPI_UNSIGNED_LONG; }
973 
974  template <typename MAT> struct mpi_distributed_matrix {
975  MAT M;
976 
977  mpi_distributed_matrix(size_type n, size_type m) : M(n, m) {}
978  mpi_distributed_matrix() {}
979 
980  const MAT &local_matrix(void) const { return M; }
981  MAT &local_matrix(void) { return M; }
982  };
983 
984  template <typename MAT> inline MAT &eff_matrix(MAT &m) { return m; }
985  template <typename MAT> inline
986  const MAT &eff_matrix(const MAT &m) { return m; }
987  template <typename MAT> inline
988  MAT &eff_matrix(mpi_distributed_matrix<MAT> &m) { return m.M; }
989  template <typename MAT> inline
990  const MAT &eff_matrix(const mpi_distributed_matrix<MAT> &m) { return m.M; }
991 
992 
993  template <typename MAT1, typename MAT2>
994  inline void copy(const mpi_distributed_matrix<MAT1> &m1,
995  mpi_distributed_matrix<MAT2> &m2)
996  { copy(eff_matrix(m1), eff_matrix(m2)); }
997  template <typename MAT1, typename MAT2>
998  inline void copy(const mpi_distributed_matrix<MAT1> &m1,
999  const mpi_distributed_matrix<MAT2> &m2)
1000  { copy(m1.M, m2.M); }
1001 
1002  template <typename MAT1, typename MAT2>
1003  inline void copy(const mpi_distributed_matrix<MAT1> &m1, MAT2 &m2)
1004  { copy(m1.M, m2); }
1005  template <typename MAT1, typename MAT2>
1006  inline void copy(const mpi_distributed_matrix<MAT1> &m1, const MAT2 &m2)
1007  { copy(m1.M, m2); }
1008 
1009 
1010  template <typename MATSP, typename V1, typename V2> inline
1011  typename strongest_value_type3<V1,V2,MATSP>::value_type
1012  vect_sp(const mpi_distributed_matrix<MATSP> &ps, const V1 &v1,
1013  const V2 &v2) {
1014  typedef typename strongest_value_type3<V1,V2,MATSP>::value_type T;
1015  T res = vect_sp(ps.M, v1, v2), rest;
1016  MPI_Allreduce(&res, &rest, 1, mpi_type(T()), MPI_SUM,MPI_COMM_WORLD);
1017  return rest;
1018  }
1019 
1020  template <typename MAT, typename V1, typename V2>
1021  inline void mult_add(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1022  V2 &v2) {
1023  typedef typename linalg_traits<V2>::value_type T;
1024  std::vector<T> v3(vect_size(v2)), v4(vect_size(v2));
1025  static double tmult_tot = 0.0;
1026  static double tmult_tot2 = 0.0;
1027  double t_ref = MPI_Wtime();
1028  gmm::mult(m.M, v1, v3);
1029  if (is_sparse(v2)) GMM_WARNING2("Using a plain temporary, here.");
1030  double t_ref2 = MPI_Wtime();
1031  MPI_Allreduce(&(v3[0]), &(v4[0]),gmm::vect_size(v2), mpi_type(T()),
1032  MPI_SUM,MPI_COMM_WORLD);
1033  tmult_tot2 = MPI_Wtime()-t_ref2;
1034  cout << "reduce mult mpi = " << tmult_tot2 << endl;
1035  gmm::add(v4, v2);
1036  tmult_tot = MPI_Wtime()-t_ref;
1037  cout << "tmult mpi = " << tmult_tot << endl;
1038  }
1039 
1040  template <typename MAT, typename V1, typename V2>
1041  void mult_add(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1042  const V2 &v2_)
1043  { mult_add(m, v1, const_cast<V2 &>(v2_)); }
1044 
1045  template <typename MAT, typename V1, typename V2>
1046  inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1047  const V2 &v2_)
1048  { V2 &v2 = const_cast<V2 &>(v2_); clear(v2); mult_add(m, v1, v2); }
1049 
1050  template <typename MAT, typename V1, typename V2>
1051  inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1052  V2 &v2)
1053  { clear(v2); mult_add(m, v1, v2); }
1054 
1055  template <typename MAT, typename V1, typename V2, typename V3>
1056  inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1057  const V2 &v2, const V3 &v3_)
1058  { V3 &v3 = const_cast<V3 &>(v3_); gmm::copy(v2, v3); mult_add(m, v1, v3); }
1059 
1060  template <typename MAT, typename V1, typename V2, typename V3>
1061  inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1,
1062  const V2 &v2, V3 &v3)
1063  { gmm::copy(v2, v3); mult_add(m, v1, v3); }
1064 
1065 
1066  template <typename MAT> inline
1067  size_type mat_nrows(const mpi_distributed_matrix<MAT> &M)
1068  { return mat_nrows(M.M); }
1069  template <typename MAT> inline
1070  size_type mat_ncols(const mpi_distributed_matrix<MAT> &M)
1071  { return mat_nrows(M.M); }
1072  template <typename MAT> inline
1073  void resize(mpi_distributed_matrix<MAT> &M, size_type m, size_type n)
1074  { resize(M.M, m, n); }
1075  template <typename MAT> inline void clear(mpi_distributed_matrix<MAT> &M)
1076  { clear(M.M); }
1077 
1078 
1079  // For compute reduced system
1080  template <typename MAT1, typename MAT2> inline
1081  void mult(const MAT1 &M1, const mpi_distributed_matrix<MAT2> &M2,
1082  mpi_distributed_matrix<MAT2> &M3)
1083  { mult(M1, M2.M, M3.M); }
1084  template <typename MAT1, typename MAT2> inline
1085  void mult(const mpi_distributed_matrix<MAT2> &M2,
1086  const MAT1 &M1, mpi_distributed_matrix<MAT2> &M3)
1087  { mult(M2.M, M1, M3.M); }
1088  template <typename MAT1, typename MAT2, typename MAT3> inline
1089  void mult(const MAT1 &M1, const mpi_distributed_matrix<MAT2> &M2,
1090  MAT3 &M3)
1091  { mult(M1, M2.M, M3); }
1092  template <typename MAT1, typename MAT2, typename MAT3> inline
1093  void mult(const MAT1 &M1, const mpi_distributed_matrix<MAT2> &M2,
1094  const MAT3 &M3)
1095  { mult(M1, M2.M, M3); }
1096 
1097  template <typename M, typename SUBI1, typename SUBI2>
1098  struct sub_matrix_type<const mpi_distributed_matrix<M> *, SUBI1, SUBI2>
1099  { typedef abstract_null_type matrix_type; };
1100 
1101  template <typename M, typename SUBI1, typename SUBI2>
1102  struct sub_matrix_type<mpi_distributed_matrix<M> *, SUBI1, SUBI2>
1103  { typedef abstract_null_type matrix_type; };
1104 
1105  template <typename M, typename SUBI1, typename SUBI2> inline
1106  typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI2>
1107  ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI2>::matrix_type,
1108  M *>::return_type
1109  sub_matrix(mpi_distributed_matrix<M> &m, const SUBI1 &si1, const SUBI2 &si2)
1110  { return sub_matrix(m.M, si1, si2); }
1111 
1112  template <typename MAT, typename SUBI1, typename SUBI2> inline
1113  typename select_return<typename sub_matrix_type<const MAT *, SUBI1, SUBI2>
1114  ::matrix_type, typename sub_matrix_type<MAT *, SUBI1, SUBI2>::matrix_type,
1115  const MAT *>::return_type
1116  sub_matrix(const mpi_distributed_matrix<MAT> &m, const SUBI1 &si1,
1117  const SUBI2 &si2)
1118  { return sub_matrix(m.M, si1, si2); }
1119 
1120  template <typename M, typename SUBI1> inline
1121  typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI1>
1122  ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type,
1123  M *>::return_type
1124  sub_matrix(mpi_distributed_matrix<M> &m, const SUBI1 &si1)
1125  { return sub_matrix(m.M, si1, si1); }
1126 
1127  template <typename M, typename SUBI1> inline
1128  typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI1>
1129  ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type,
1130  const M *>::return_type
1131  sub_matrix(const mpi_distributed_matrix<M> &m, const SUBI1 &si1)
1132  { return sub_matrix(m.M, si1, si1); }
1133 
1134 
1135  template <typename L> struct transposed_return<const mpi_distributed_matrix<L> *>
1136  { typedef abstract_null_type return_type; };
1137  template <typename L> struct transposed_return<mpi_distributed_matrix<L> *>
1138  { typedef abstract_null_type return_type; };
1139 
1140  template <typename L> inline typename transposed_return<const L *>::return_type
1141  transposed(const mpi_distributed_matrix<L> &l)
1142  { return transposed(l.M); }
1143 
1144  template <typename L> inline typename transposed_return<L *>::return_type
1145  transposed(mpi_distributed_matrix<L> &l)
1146  { return transposed(l.M); }
1147 
1148 
1149  template <typename MAT>
1150  struct linalg_traits<mpi_distributed_matrix<MAT> > {
1151  typedef mpi_distributed_matrix<MAT> this_type;
1152  typedef MAT origin_type;
1153  typedef linalg_false is_reference;
1154  typedef abstract_matrix linalg_type;
1155  typedef typename linalg_traits<MAT>::value_type value_type;
1156  typedef typename linalg_traits<MAT>::reference reference;
1157  typedef typename linalg_traits<MAT>::storage_type storage_type;
1158  typedef abstract_null_type sub_row_type;
1159  typedef abstract_null_type const_sub_row_type;
1160  typedef abstract_null_type row_iterator;
1161  typedef abstract_null_type const_row_iterator;
1162  typedef abstract_null_type sub_col_type;
1163  typedef abstract_null_type const_sub_col_type;
1164  typedef abstract_null_type col_iterator;
1165  typedef abstract_null_type const_col_iterator;
1166  typedef abstract_null_type sub_orientation;
1167  typedef abstract_null_type index_sorted;
1168  static size_type nrows(const this_type &m) { return nrows(m.M); }
1169  static size_type ncols(const this_type &m) { return ncols(m.M); }
1170  static void do_clear(this_type &m) { clear(m.M); }
1171  };
1172 
1173 }
1174 
1175 
1176 #endif // GMM_USES_MPI
1177 
1178 namespace std {
1179  template <typename V>
1180  void swap(gmm::row_matrix<V> &m1, gmm::row_matrix<V> &m2)
1181  { m1.swap(m2); }
1182  template <typename V>
1183  void swap(gmm::col_matrix<V> &m1, gmm::col_matrix<V> &m2)
1184  { m1.swap(m2); }
1185  template <typename T>
1186  void swap(gmm::dense_matrix<T> &m1, gmm::dense_matrix<T> &m2)
1187  { m1.swap(m2); }
1188  template <typename T, int shift> void
1189  swap(gmm::csc_matrix<T,shift> &m1, gmm::csc_matrix<T,shift> &m2)
1190  { m1.swap(m2); }
1191  template <typename T, int shift> void
1192  swap(gmm::csr_matrix<T,shift> &m1, gmm::csr_matrix<T,shift> &m2)
1193  { m1.swap(m2); }
1194 }
1195 
1196 
1197 
1198 
1199 #endif /* GMM_MATRIX_H__ */
gmm_sub_matrix.h
Generic sub-matrices.
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::copy
void copy(const L1 &l1, L2 &l2)
*/
Definition: gmm_blas.h:977
gmm::resize
void resize(V &v, size_type n)
*/
Definition: gmm_blas.h:209
gmm_vector.h
Declaration of the vector types (gmm::rsvector, gmm::wsvector, gmm::slvector ,..) ...
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_def.h:488
gmm_sub_vector.h
Generic sub-vectors.
gmm::fill
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*/
Definition: gmm_blas.h:103
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*/
Definition: gmm_blas.h:263
gmm::reshape
void reshape(M &v, size_type m, size_type n)
*/
Definition: gmm_blas.h:250
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::mult_add
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*/
Definition: gmm_blas.h:1779
gmm::vect_hp
strongest_value_type< V1, V2 >::value_type vect_hp(const V1 &v1, const V2 &v2)
*/
Definition: gmm_blas.h:511
gmm::nnz
size_type nnz(const L &l)
count the number of non-zero entries of a vector or matrix.
Definition: gmm_blas.h:68
gmm_transposed.h
Generic transposed matrices.
gmm::add
void add(const L1 &l1, L2 &l2)
*/
Definition: gmm_blas.h:1268