GetFEM++  5.3
gmm_blas.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_blas.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>
34  @date October 13, 2002.
35  @brief Basic linear algebra functions.
36 */
37 
38 #ifndef GMM_BLAS_H__
39 #define GMM_BLAS_H__
40 
41 #include "gmm_scaled.h"
42 #include "gmm_transposed.h"
43 #include "gmm_conjugated.h"
44 
45 namespace gmm {
46 
47  /* ******************************************************************** */
48  /* */
49  /* Generic algorithms */
50  /* */
51  /* ******************************************************************** */
52 
53 
54  /* ******************************************************************** */
55  /* Miscellaneous */
56  /* ******************************************************************** */
57 
58  /** clear (fill with zeros) a vector or matrix. */
59  template <typename L> inline void clear(L &l)
60  { linalg_traits<L>::do_clear(l); }
61  /** @cond DOXY_SHOW_ALL_FUNCTIONS
62  skip all these redundant definitions in doxygen documentation..
63  */
64  template <typename L> inline void clear(const L &l)
65  { linalg_traits<L>::do_clear(linalg_const_cast(l)); }
66 
67  ///@endcond
68  /** count the number of non-zero entries of a vector or matrix. */ template <typename L> inline size_type nnz(const L& l)
69  { return nnz(l, typename linalg_traits<L>::linalg_type()); }
70 
71  ///@cond DOXY_SHOW_ALL_FUNCTIONS
72  template <typename L> inline size_type nnz(const L& l, abstract_vector) {
73  auto it = vect_const_begin(l), ite = vect_const_end(l);
74  size_type res(0);
75  for (; it != ite; ++it) ++res;
76  return res;
77  }
78 
79  template <typename L> inline size_type nnz(const L& l, abstract_matrix) {
80  return nnz(l, typename principal_orientation_type<typename
81  linalg_traits<L>::sub_orientation>::potype());
82  }
83 
84  template <typename L> inline size_type nnz(const L& l, row_major) {
85  size_type res(0);
86  for (size_type i = 0; i < mat_nrows(l); ++i)
87  res += nnz(mat_const_row(l, i));
88  return res;
89  }
90 
91  template <typename L> inline size_type nnz(const L& l, col_major) {
92  size_type res(0);
93  for (size_type i = 0; i < mat_ncols(l); ++i)
94  res += nnz(mat_const_col(l, i));
95  return res;
96  }
97 
98  ///@endcond
99 
100 
101  /** fill a vector or matrix with x. */
102  template <typename L> inline
103  void fill(L& l, typename gmm::linalg_traits<L>::value_type x) {
104  typedef typename gmm::linalg_traits<L>::value_type T;
105  if (x == T(0)) gmm::clear(l);
106  fill(l, x, typename linalg_traits<L>::linalg_type());
107  }
108 
109  template <typename L> inline
110  void fill(const L& l, typename gmm::linalg_traits<L>::value_type x) {
111  fill(linalg_const_cast(l), x);
112  }
113 
114  template <typename L> inline // to be optimized for dense vectors ...
115  void fill(L& l, typename gmm::linalg_traits<L>::value_type x,
116  abstract_vector) {
117  for (size_type i = 0; i < vect_size(l); ++i) l[i] = x;
118  }
119 
120  template <typename L> inline // to be optimized for dense matrices ...
121  void fill(L& l, typename gmm::linalg_traits<L>::value_type x,
122  abstract_matrix) {
123  for (size_type i = 0; i < mat_nrows(l); ++i)
124  for (size_type j = 0; j < mat_ncols(l); ++j)
125  l(i,j) = x;
126  }
127 
128  /** fill a vector or matrix with random value (uniform [-1,1]). */
129  template <typename L> inline void fill_random(L& l)
130  { fill_random(l, typename linalg_traits<L>::linalg_type()); }
131 
132  ///@cond DOXY_SHOW_ALL_FUNCTIONS
133  template <typename L> inline void fill_random(const L& l) {
134  fill_random(linalg_const_cast(l),
135  typename linalg_traits<L>::linalg_type());
136  }
137 
138  template <typename L> inline void fill_random(L& l, abstract_vector) {
139  for (size_type i = 0; i < vect_size(l); ++i)
140  l[i] = gmm::random(typename linalg_traits<L>::value_type());
141  }
142 
143  template <typename L> inline void fill_random(L& l, abstract_matrix) {
144  for (size_type i = 0; i < mat_nrows(l); ++i)
145  for (size_type j = 0; j < mat_ncols(l); ++j)
146  l(i,j) = gmm::random(typename linalg_traits<L>::value_type());
147  }
148 
149  ///@endcond
150  /** fill a vector or matrix with random value.
151  @param l a vector or matrix.
152  @param cfill probability of a non-zero value.
153  */
154  template <typename L> inline void fill_random(L& l, double cfill)
155  { fill_random(l, cfill, typename linalg_traits<L>::linalg_type()); }
156  ///@cond DOXY_SHOW_ALL_FUNCTIONS
157 
158  template <typename L> inline void fill_random(const L& l, double cfill) {
159  fill_random(linalg_const_cast(l), cfill,
160  typename linalg_traits<L>::linalg_type());
161  }
162 
163  template <typename L> inline
164  void fill_random(L& l, double cfill, abstract_vector) {
165  typedef typename linalg_traits<L>::value_type T;
166  size_type ntot = std::min(vect_size(l),
167  size_type(double(vect_size(l))*cfill) + 1);
168  for (size_type nb = 0; nb < ntot;) {
169  size_type i = gmm::irandom(vect_size(l));
170  if (l[i] == T(0)) {
171  l[i] = gmm::random(typename linalg_traits<L>::value_type());
172  ++nb;
173  }
174  }
175  }
176 
177  template <typename L> inline
178  void fill_random(L& l, double cfill, abstract_matrix) {
179  fill_random(l, cfill, typename principal_orientation_type<typename
180  linalg_traits<L>::sub_orientation>::potype());
181  }
182 
183  template <typename L> inline
184  void fill_random(L& l, double cfill, row_major) {
185  for (size_type i=0; i < mat_nrows(l); ++i) fill_random(mat_row(l,i),cfill);
186  }
187 
188  template <typename L> inline
189  void fill_random(L& l, double cfill, col_major) {
190  for (size_type j=0; j < mat_ncols(l); ++j) fill_random(mat_col(l,j),cfill);
191  }
192 
193  /* resize a vector */
194  template <typename V> inline
195  void resize(V &v, size_type n, linalg_false)
196  { linalg_traits<V>::resize(v, n); }
197 
198  template <typename V> inline
199  void resize(V &, size_type , linalg_modifiable)
200  { GMM_ASSERT1(false, "You cannot resize a reference"); }
201 
202  template <typename V> inline
203  void resize(V &, size_type , linalg_const)
204  { GMM_ASSERT1(false, "You cannot resize a reference"); }
205 
206  ///@endcond
207  /** resize a vector. */
208  template <typename V> inline
209  void resize(V &v, size_type n) {
210  resize(v, n, typename linalg_traits<V>::is_reference());
211  }
212  ///@cond DOXY_SHOW_ALL_FUNCTIONS
213 
214  /** resize a matrix **/
215  template <typename M> inline
216  void resize(M &v, size_type m, size_type n, linalg_false) {
217  linalg_traits<M>::resize(v, m, n);
218  }
219 
220  template <typename M> inline
221  void resize(M &, size_type, size_type, linalg_modifiable)
222  { GMM_ASSERT1(false, "You cannot resize a reference"); }
223 
224  template <typename M> inline
225  void resize(M &, size_type, size_type, linalg_const)
226  { GMM_ASSERT1(false, "You cannot resize a reference"); }
227 
228  ///@endcond
229  /** resize a matrix */
230  template <typename M> inline
231  void resize(M &v, size_type m, size_type n)
232  { resize(v, m, n, typename linalg_traits<M>::is_reference()); }
233  ///@cond
234 
235  template <typename M> inline
236  void reshape(M &v, size_type m, size_type n, linalg_false)
237  { linalg_traits<M>::reshape(v, m, n); }
238 
239  template <typename M> inline
240  void reshape(M &, size_type, size_type, linalg_modifiable)
241  { GMM_ASSERT1(false, "You cannot reshape a reference"); }
242 
243  template <typename M> inline
244  void reshape(M &, size_type, size_type, linalg_const)
245  { GMM_ASSERT1(false, "You cannot reshape a reference"); }
246 
247  ///@endcond
248  /** reshape a matrix */
249  template <typename M> inline
250  void reshape(M &v, size_type m, size_type n)
251  { reshape(v, m, n, typename linalg_traits<M>::is_reference()); }
252  ///@cond DOXY_SHOW_ALL_FUNCTIONS
253 
254 
255  /* ******************************************************************** */
256  /* Scalar product */
257  /* ******************************************************************** */
258 
259  ///@endcond
260  /** scalar product between two vectors */
261  template <typename V1, typename V2> inline
262  typename strongest_value_type<V1,V2>::value_type
263  vect_sp(const V1 &v1, const V2 &v2) {
264  GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch, "
265  << vect_size(v1) << " !=" << vect_size(v2));
266  return vect_sp(v1, v2,
267  typename linalg_traits<V1>::storage_type(),
268  typename linalg_traits<V2>::storage_type());
269  }
270 
271  /** scalar product between two vectors, using a matrix.
272  @param ps the matrix of the scalar product.
273  @param v1 the first vector
274  @param v2 the second vector
275  */
276  template <typename MATSP, typename V1, typename V2> inline
277  typename strongest_value_type3<V1,V2,MATSP>::value_type
278  vect_sp(const MATSP &ps, const V1 &v1, const V2 &v2) {
279  return vect_sp_with_mat(ps, v1, v2,
280  typename linalg_traits<MATSP>::sub_orientation());
281  }
282  ///@cond DOXY_SHOW_ALL_FUNCTIONS
283 
284  template <typename MATSP, typename V1, typename V2> inline
285  typename strongest_value_type3<V1,V2,MATSP>::value_type
286  vect_sp_with_mat(const MATSP &ps, const V1 &v1, const V2 &v2, row_major) {
287  return vect_sp_with_matr(ps, v1, v2,
288  typename linalg_traits<V2>::storage_type());
289  }
290 
291  template <typename MATSP, typename V1, typename V2> inline
292  typename strongest_value_type3<V1,V2,MATSP>::value_type
293  vect_sp_with_matr(const MATSP &ps, const V1 &v1, const V2 &v2,
294  abstract_sparse) {
295  GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) &&
296  vect_size(v2) == mat_nrows(ps), "dimensions mismatch");
297  size_type nr = mat_nrows(ps);
298  typename linalg_traits<V2>::const_iterator
299  it = vect_const_begin(v2), ite = vect_const_end(v2);
300  typename strongest_value_type3<V1,V2,MATSP>::value_type res(0);
301  for (; it != ite; ++it)
302  res += vect_sp(mat_const_row(ps, it.index()), v1)* (*it);
303  return res;
304  }
305 
306  template <typename MATSP, typename V1, typename V2> inline
307  typename strongest_value_type3<V1,V2,MATSP>::value_type
308  vect_sp_with_matr(const MATSP &ps, const V1 &v1, const V2 &v2,
309  abstract_skyline)
310  { return vect_sp_with_matr(ps, v1, v2, abstract_sparse()); }
311 
312  template <typename MATSP, typename V1, typename V2> inline
313  typename strongest_value_type3<V1,V2,MATSP>::value_type
314  vect_sp_with_matr(const MATSP &ps, const V1 &v1, const V2 &v2,
315  abstract_dense) {
316  GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) &&
317  vect_size(v2) == mat_nrows(ps), "dimensions mismatch");
318  typename linalg_traits<V2>::const_iterator
319  it = vect_const_begin(v2), ite = vect_const_end(v2);
320  typename strongest_value_type3<V1,V2,MATSP>::value_type res(0);
321  for (size_type i = 0; it != ite; ++i, ++it)
322  res += vect_sp(mat_const_row(ps, i), v1) * (*it);
323  return res;
324  }
325 
326  template <typename MATSP, typename V1, typename V2> inline
327  typename strongest_value_type3<V1,V2,MATSP>::value_type
328  vect_sp_with_mat(const MATSP &ps, const V1 &v1,const V2 &v2,row_and_col)
329  { return vect_sp_with_mat(ps, v1, v2, row_major()); }
330 
331  template <typename MATSP, typename V1, typename V2> inline
332  typename strongest_value_type3<V1,V2,MATSP>::value_type
333  vect_sp_with_mat(const MATSP &ps, const V1 &v1, const V2 &v2,col_major){
334  return vect_sp_with_matc(ps, v1, v2,
335  typename linalg_traits<V1>::storage_type());
336  }
337 
338  template <typename MATSP, typename V1, typename V2> inline
339  typename strongest_value_type3<V1,V2,MATSP>::value_type
340  vect_sp_with_matc(const MATSP &ps, const V1 &v1, const V2 &v2,
341  abstract_sparse) {
342  GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) &&
343  vect_size(v2) == mat_nrows(ps), "dimensions mismatch");
344  typename linalg_traits<V1>::const_iterator
345  it = vect_const_begin(v1), ite = vect_const_end(v1);
346  typename strongest_value_type3<V1,V2,MATSP>::value_type res(0);
347  for (; it != ite; ++it)
348  res += vect_sp(mat_const_col(ps, it.index()), v2) * (*it);
349  return res;
350  }
351 
352  template <typename MATSP, typename V1, typename V2> inline
353  typename strongest_value_type3<V1,V2,MATSP>::value_type
354  vect_sp_with_matc(const MATSP &ps, const V1 &v1, const V2 &v2,
355  abstract_skyline)
356  { return vect_sp_with_matc(ps, v1, v2, abstract_sparse()); }
357 
358  template <typename MATSP, typename V1, typename V2> inline
359  typename strongest_value_type3<V1,V2,MATSP>::value_type
360  vect_sp_with_matc(const MATSP &ps, const V1 &v1, const V2 &v2,
361  abstract_dense) {
362  GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) &&
363  vect_size(v2) == mat_nrows(ps), "dimensions mismatch");
364  typename linalg_traits<V1>::const_iterator
365  it = vect_const_begin(v1), ite = vect_const_end(v1);
366  typename strongest_value_type3<V1,V2,MATSP>::value_type res(0);
367  for (size_type i = 0; it != ite; ++i, ++it)
368  res += vect_sp(mat_const_col(ps, i), v2) * (*it);
369  return res;
370  }
371 
372  template <typename MATSP, typename V1, typename V2> inline
373  typename strongest_value_type3<V1,V2,MATSP>::value_type
374  vect_sp_with_mat(const MATSP &ps, const V1 &v1,const V2 &v2,col_and_row)
375  { return vect_sp_with_mat(ps, v1, v2, col_major()); }
376 
377  template <typename MATSP, typename V1, typename V2> inline
378  typename strongest_value_type3<V1,V2,MATSP>::value_type
379  vect_sp_with_mat(const MATSP &ps, const V1 &v1, const V2 &v2,
380  abstract_null_type) {
381  typename temporary_vector<V1>::vector_type w(mat_nrows(ps));
382  GMM_WARNING2("Warning, a temporary is used in scalar product\n");
383  mult(ps, v1, w);
384  return vect_sp(w, v2);
385  }
386 
387  template <typename IT1, typename IT2> inline
388  typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
389  typename std::iterator_traits<IT2>::value_type>::T
390  vect_sp_dense_(IT1 it, IT1 ite, IT2 it2) {
391  typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
392  typename std::iterator_traits<IT2>::value_type>::T res(0);
393  for (; it != ite; ++it, ++it2) res += (*it) * (*it2);
394  return res;
395  }
396 
397  template <typename IT1, typename V> inline
398  typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
399  typename linalg_traits<V>::value_type>::T
400  vect_sp_sparse_(IT1 it, IT1 ite, const V &v) {
401  typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
402  typename linalg_traits<V>::value_type>::T res(0);
403  for (; it != ite; ++it) res += (*it) * v[it.index()];
404  return res;
405  }
406 
407  template <typename V1, typename V2> inline
408  typename strongest_value_type<V1,V2>::value_type
409  vect_sp(const V1 &v1, const V2 &v2, abstract_dense, abstract_dense) {
410  return vect_sp_dense_(vect_const_begin(v1), vect_const_end(v1),
411  vect_const_begin(v2));
412  }
413 
414  template <typename V1, typename V2> inline
415  typename strongest_value_type<V1,V2>::value_type
416  vect_sp(const V1 &v1, const V2 &v2, abstract_skyline, abstract_dense) {
417  typename linalg_traits<V1>::const_iterator it1 = vect_const_begin(v1),
418  ite = vect_const_end(v1);
419  typename linalg_traits<V2>::const_iterator it2 = vect_const_begin(v2);
420  return vect_sp_dense_(it1, ite, it2 + it1.index());
421  }
422 
423  template <typename V1, typename V2> inline
424  typename strongest_value_type<V1,V2>::value_type
425  vect_sp(const V1 &v1, const V2 &v2, abstract_dense, abstract_skyline) {
426  typename linalg_traits<V2>::const_iterator it1 = vect_const_begin(v2),
427  ite = vect_const_end(v2);
428  typename linalg_traits<V1>::const_iterator it2 = vect_const_begin(v1);
429  return vect_sp_dense_(it1, ite, it2 + it1.index());
430  }
431 
432  template <typename V1, typename V2> inline
433  typename strongest_value_type<V1,V2>::value_type
434  vect_sp(const V1 &v1, const V2 &v2, abstract_skyline, abstract_skyline) {
435  typedef typename strongest_value_type<V1,V2>::value_type T;
436  auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
437  auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
438  size_type n = std::min(ite1.index(), ite2.index());
439  size_type l = std::max(it1.index(), it2.index());
440 
441  if (l < n) {
442  size_type m = l - it1.index(), p = l - it2.index(), q = m + n - l;
443  return vect_sp_dense_(it1+m, it1+q, it2 + p);
444  }
445  return T(0);
446  }
447 
448  template <typename V1, typename V2> inline
449  typename strongest_value_type<V1,V2>::value_type
450  vect_sp(const V1 &v1, const V2 &v2,abstract_sparse,abstract_dense) {
451  return vect_sp_sparse_(vect_const_begin(v1), vect_const_end(v1), v2);
452  }
453 
454  template <typename V1, typename V2> inline
455  typename strongest_value_type<V1,V2>::value_type
456  vect_sp(const V1 &v1, const V2 &v2, abstract_sparse, abstract_skyline) {
457  return vect_sp_sparse_(vect_const_begin(v1), vect_const_end(v1), v2);
458  }
459 
460  template <typename V1, typename V2> inline
461  typename strongest_value_type<V1,V2>::value_type
462  vect_sp(const V1 &v1, const V2 &v2, abstract_skyline, abstract_sparse) {
463  return vect_sp_sparse_(vect_const_begin(v2), vect_const_end(v2), v1);
464  }
465 
466  template <typename V1, typename V2> inline
467  typename strongest_value_type<V1,V2>::value_type
468  vect_sp(const V1 &v1, const V2 &v2, abstract_dense,abstract_sparse) {
469  return vect_sp_sparse_(vect_const_begin(v2), vect_const_end(v2), v1);
470  }
471 
472 
473  template <typename V1, typename V2> inline
474  typename strongest_value_type<V1,V2>::value_type
475  vect_sp_sparse_sparse(const V1 &v1, const V2 &v2, linalg_true) {
476  typename linalg_traits<V1>::const_iterator it1 = vect_const_begin(v1),
477  ite1 = vect_const_end(v1);
478  typename linalg_traits<V2>::const_iterator it2 = vect_const_begin(v2),
479  ite2 = vect_const_end(v2);
480  typename strongest_value_type<V1,V2>::value_type res(0);
481 
482  while (it1 != ite1 && it2 != ite2) {
483  if (it1.index() == it2.index())
484  { res += (*it1) * *it2; ++it1; ++it2; }
485  else if (it1.index() < it2.index()) ++it1; else ++it2;
486  }
487  return res;
488  }
489 
490  template <typename V1, typename V2> inline
491  typename strongest_value_type<V1,V2>::value_type
492  vect_sp_sparse_sparse(const V1 &v1, const V2 &v2, linalg_false) {
493  return vect_sp_sparse_(vect_const_begin(v1), vect_const_end(v1), v2);
494  }
495 
496  template <typename V1, typename V2> inline
497  typename strongest_value_type<V1,V2>::value_type
498  vect_sp(const V1 &v1, const V2 &v2,abstract_sparse,abstract_sparse) {
499  return vect_sp_sparse_sparse(v1, v2,
500  typename linalg_and<typename linalg_traits<V1>::index_sorted,
501  typename linalg_traits<V2>::index_sorted>::bool_type());
502  }
503 
504  /* ******************************************************************** */
505  /* Hermitian product */
506  /* ******************************************************************** */
507  ///@endcond
508  /** Hermitian product. */
509  template <typename V1, typename V2>
510  inline typename strongest_value_type<V1,V2>::value_type
511  vect_hp(const V1 &v1, const V2 &v2)
512  { return vect_sp(v1, conjugated(v2)); }
513 
514  /** Hermitian product with a matrix. */
515  template <typename MATSP, typename V1, typename V2> inline
516  typename strongest_value_type3<V1,V2,MATSP>::value_type
517  vect_hp(const MATSP &ps, const V1 &v1, const V2 &v2) {
518  return vect_sp(ps, v1, gmm::conjugated(v2));
519  }
520 
521  /* ******************************************************************** */
522  /* Trace of a matrix */
523  /* ******************************************************************** */
524 
525  /** Trace of a matrix */
526  template <typename M>
527  typename linalg_traits<M>::value_type
528  mat_trace(const M &m) {
529  typedef typename linalg_traits<M>::value_type T;
530  T res(0);
531  for (size_type i = 0; i < std::min(mat_nrows(m), mat_ncols(m)); ++i)
532  res += m(i,i);
533  return res;
534  }
535 
536  /* ******************************************************************** */
537  /* Euclidean norm */
538  /* ******************************************************************** */
539 
540  /** squared Euclidean norm of a vector. */
541  template <typename V>
542  typename number_traits<typename linalg_traits<V>::value_type>
543  ::magnitude_type
544  vect_norm2_sqr(const V &v) {
545  typedef typename linalg_traits<V>::value_type T;
546  typedef typename number_traits<T>::magnitude_type R;
547  auto it = vect_const_begin(v), ite = vect_const_end(v);
548  R res(0);
549  for (; it != ite; ++it) res += gmm::abs_sqr(*it);
550  return res;
551  }
552 
553  /** Euclidean norm of a vector. */
554  template <typename V> inline
555  typename number_traits<typename linalg_traits<V>::value_type>
556  ::magnitude_type
557  vect_norm2(const V &v)
558  { return sqrt(vect_norm2_sqr(v)); }
559 
560 
561  /** squared Euclidean distance between two vectors */
562  template <typename V1, typename V2> inline
563  typename number_traits<typename linalg_traits<V1>::value_type>
564  ::magnitude_type
565  vect_dist2_sqr(const V1 &v1, const V2 &v2) { // not fully optimized
566  typedef typename linalg_traits<V1>::value_type T;
567  typedef typename number_traits<T>::magnitude_type R;
568  auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
569  auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
570  size_type k1(0), k2(0);
571  R res(0);
572  while (it1 != ite1 && it2 != ite2) {
573  size_type i1 = index_of_it(it1, k1,
574  typename linalg_traits<V1>::storage_type());
575  size_type i2 = index_of_it(it2, k2,
576  typename linalg_traits<V2>::storage_type());
577 
578  if (i1 == i2) {
579  res += gmm::abs_sqr(*it2 - *it1); ++it1; ++k1; ++it2; ++k2;
580  }
581  else if (i1 < i2) {
582  res += gmm::abs_sqr(*it1); ++it1; ++k1;
583  }
584  else {
585  res += gmm::abs_sqr(*it2); ++it2; ++k2;
586  }
587  }
588  while (it1 != ite1) { res += gmm::abs_sqr(*it1); ++it1; }
589  while (it2 != ite2) { res += gmm::abs_sqr(*it2); ++it2; }
590  return res;
591  }
592 
593  /** Euclidean distance between two vectors */
594  template <typename V1, typename V2> inline
595  typename number_traits<typename linalg_traits<V1>::value_type>
596  ::magnitude_type
597  vect_dist2(const V1 &v1, const V2 &v2)
598  { return sqrt(vect_dist2_sqr(v1, v2)); }
599  ///@cond DOXY_SHOW_ALL_FUNCTIONS
600  template <typename M>
601  typename number_traits<typename linalg_traits<M>::value_type>
602  ::magnitude_type
603  mat_euclidean_norm_sqr(const M &m, row_major) {
604  typename number_traits<typename linalg_traits<M>::value_type>
605  ::magnitude_type res(0);
606  for (size_type i = 0; i < mat_nrows(m); ++i)
607  res += vect_norm2_sqr(mat_const_row(m, i));
608  return res;
609  }
610 
611  template <typename M>
612  typename number_traits<typename linalg_traits<M>::value_type>
613  ::magnitude_type
614  mat_euclidean_norm_sqr(const M &m, col_major) {
615  typename number_traits<typename linalg_traits<M>::value_type>
616  ::magnitude_type res(0);
617  for (size_type i = 0; i < mat_ncols(m); ++i)
618  res += vect_norm2_sqr(mat_const_col(m, i));
619  return res;
620  }
621  ///@endcond
622  /** squared Euclidean norm of a matrix. */
623  template <typename M> inline
624  typename number_traits<typename linalg_traits<M>::value_type>
625  ::magnitude_type
626  mat_euclidean_norm_sqr(const M &m) {
627  return mat_euclidean_norm_sqr(m,
628  typename principal_orientation_type<typename
629  linalg_traits<M>::sub_orientation>::potype());
630  }
631 
632  /** Euclidean norm of a matrix. */
633  template <typename M> inline
634  typename number_traits<typename linalg_traits<M>::value_type>
635  ::magnitude_type
636  mat_euclidean_norm(const M &m)
637  { return gmm::sqrt(mat_euclidean_norm_sqr(m)); }
638 
639  /* ******************************************************************** */
640  /* vector norm1 */
641  /* ******************************************************************** */
642  /** 1-norm of a vector */
643  template <typename V>
644  typename number_traits<typename linalg_traits<V>::value_type>
645  ::magnitude_type
646  vect_norm1(const V &v) {
647  auto it = vect_const_begin(v), ite = vect_const_end(v);
648  typename number_traits<typename linalg_traits<V>::value_type>
649  ::magnitude_type res(0);
650  for (; it != ite; ++it) res += gmm::abs(*it);
651  return res;
652  }
653 
654  /** 1-distance between two vectors */
655  template <typename V1, typename V2> inline
656  typename number_traits<typename linalg_traits<V1>::value_type>
657  ::magnitude_type
658  vect_dist1(const V1 &v1, const V2 &v2) { // not fully optimized
659  typedef typename linalg_traits<V1>::value_type T;
660  typedef typename number_traits<T>::magnitude_type R;
661  auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
662  auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
663  size_type k1(0), k2(0);
664  R res(0);
665  while (it1 != ite1 && it2 != ite2) {
666  size_type i1 = index_of_it(it1, k1,
667  typename linalg_traits<V1>::storage_type());
668  size_type i2 = index_of_it(it2, k2,
669  typename linalg_traits<V2>::storage_type());
670 
671  if (i1 == i2) {
672  res += gmm::abs(*it2 - *it1); ++it1; ++k1; ++it2; ++k2;
673  }
674  else if (i1 < i2) {
675  res += gmm::abs(*it1); ++it1; ++k1;
676  }
677  else {
678  res += gmm::abs(*it2); ++it2; ++k2;
679  }
680  }
681  while (it1 != ite1) { res += gmm::abs(*it1); ++it1; }
682  while (it2 != ite2) { res += gmm::abs(*it2); ++it2; }
683  return res;
684  }
685 
686  /* ******************************************************************** */
687  /* vector Infinity norm */
688  /* ******************************************************************** */
689  /** Infinity norm of a vector. */
690  template <typename V>
691  typename number_traits<typename linalg_traits<V>::value_type>
692  ::magnitude_type
693  vect_norminf(const V &v) {
694  auto it = vect_const_begin(v), ite = vect_const_end(v);
695  typename number_traits<typename linalg_traits<V>::value_type>
696  ::magnitude_type res(0);
697  for (; it != ite; ++it) res = std::max(res, gmm::abs(*it));
698  return res;
699  }
700 
701  /** Infinity distance between two vectors */
702  template <typename V1, typename V2> inline
703  typename number_traits<typename linalg_traits<V1>::value_type>
704  ::magnitude_type
705  vect_distinf(const V1 &v1, const V2 &v2) { // not fully optimized
706  typedef typename linalg_traits<V1>::value_type T;
707  typedef typename number_traits<T>::magnitude_type R;
708  auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
709  auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
710  size_type k1(0), k2(0);
711  R res(0);
712  while (it1 != ite1 && it2 != ite2) {
713  size_type i1 = index_of_it(it1, k1,
714  typename linalg_traits<V1>::storage_type());
715  size_type i2 = index_of_it(it2, k2,
716  typename linalg_traits<V2>::storage_type());
717 
718  if (i1 == i2) {
719  res = std::max(res, gmm::abs(*it2 - *it1)); ++it1; ++k1; ++it2; ++k2;
720  }
721  else if (i1 < i2) {
722  res = std::max(res, gmm::abs(*it1)); ++it1; ++k1;
723  }
724  else {
725  res = std::max(res, gmm::abs(*it2)); ++it2; ++k2;
726  }
727  }
728  while (it1 != ite1) { res = std::max(res, gmm::abs(*it1)); ++it1; }
729  while (it2 != ite2) { res = std::max(res, gmm::abs(*it2)); ++it2; }
730  return res;
731  }
732 
733  /* ******************************************************************** */
734  /* matrix norm1 */
735  /* ******************************************************************** */
736  ///@cond DOXY_SHOW_ALL_FUNCTIONS
737  template <typename M>
738  typename number_traits<typename linalg_traits<M>::value_type>
739  ::magnitude_type
740  mat_norm1(const M &m, col_major) {
741  typename number_traits<typename linalg_traits<M>::value_type>
742  ::magnitude_type res(0);
743  for (size_type i = 0; i < mat_ncols(m); ++i)
744  res = std::max(res, vect_norm1(mat_const_col(m,i)));
745  return res;
746  }
747 
748  template <typename M>
749  typename number_traits<typename linalg_traits<M>::value_type>
750  ::magnitude_type
751  mat_norm1(const M &m, row_major) {
752  typedef typename linalg_traits<M>::value_type T;
753  typedef typename number_traits<T>::magnitude_type R;
754  typedef typename linalg_traits<M>::storage_type store_type;
755 
756  std::vector<R> aux(mat_ncols(m));
757  for (size_type i = 0; i < mat_nrows(m); ++i) {
758  typename linalg_traits<M>::const_sub_row_type row = mat_const_row(m, i);
759  auto it = vect_const_begin(row), ite = vect_const_end(row);
760  for (size_type k = 0; it != ite; ++it, ++k)
761  aux[index_of_it(it, k, store_type())] += gmm::abs(*it);
762  }
763  return vect_norminf(aux);
764  }
765 
766  template <typename M>
767  typename number_traits<typename linalg_traits<M>::value_type>
768  ::magnitude_type
769  mat_norm1(const M &m, col_and_row)
770  { return mat_norm1(m, col_major()); }
771 
772  template <typename M>
773  typename number_traits<typename linalg_traits<M>::value_type>
774  ::magnitude_type
775  mat_norm1(const M &m, row_and_col)
776  { return mat_norm1(m, col_major()); }
777  ///@endcond
778  /** 1-norm of a matrix */
779  template <typename M>
780  typename number_traits<typename linalg_traits<M>::value_type>
781  ::magnitude_type
782  mat_norm1(const M &m) {
783  return mat_norm1(m, typename linalg_traits<M>::sub_orientation());
784  }
785 
786 
787  /* ******************************************************************** */
788  /* matrix Infinity norm */
789  /* ******************************************************************** */
790  ///@cond DOXY_SHOW_ALL_FUNCTIONS
791  template <typename M>
792  typename number_traits<typename linalg_traits<M>::value_type>
793  ::magnitude_type
794  mat_norminf(const M &m, row_major) {
795  typename number_traits<typename linalg_traits<M>::value_type>
796  ::magnitude_type res(0);
797  for (size_type i = 0; i < mat_nrows(m); ++i)
798  res = std::max(res, vect_norm1(mat_const_row(m,i)));
799  return res;
800  }
801 
802  template <typename M>
803  typename number_traits<typename linalg_traits<M>::value_type>
804  ::magnitude_type
805  mat_norminf(const M &m, col_major) {
806  typedef typename linalg_traits<M>::value_type T;
807  typedef typename number_traits<T>::magnitude_type R;
808  typedef typename linalg_traits<M>::storage_type store_type;
809 
810  std::vector<R> aux(mat_nrows(m));
811  for (size_type i = 0; i < mat_ncols(m); ++i) {
812  typename linalg_traits<M>::const_sub_col_type col = mat_const_col(m, i);
813  auto it = vect_const_begin(col), ite = vect_const_end(col);
814  for (size_type k = 0; it != ite; ++it, ++k)
815  aux[index_of_it(it, k, store_type())] += gmm::abs(*it);
816  }
817  return vect_norminf(aux);
818  }
819 
820  template <typename M>
821  typename number_traits<typename linalg_traits<M>::value_type>
822  ::magnitude_type
823  mat_norminf(const M &m, col_and_row)
824  { return mat_norminf(m, row_major()); }
825 
826  template <typename M>
827  typename number_traits<typename linalg_traits<M>::value_type>
828  ::magnitude_type
829  mat_norminf(const M &m, row_and_col)
830  { return mat_norminf(m, row_major()); }
831  ///@endcond
832  /** infinity-norm of a matrix.*/
833  template <typename M>
834  typename number_traits<typename linalg_traits<M>::value_type>
835  ::magnitude_type
836  mat_norminf(const M &m) {
837  return mat_norminf(m, typename linalg_traits<M>::sub_orientation());
838  }
839 
840  /* ******************************************************************** */
841  /* Max norm for matrices */
842  /* ******************************************************************** */
843  ///@cond DOXY_SHOW_ALL_FUNCTIONS
844  template <typename M>
845  typename number_traits<typename linalg_traits<M>::value_type>
846  ::magnitude_type
847  mat_maxnorm(const M &m, row_major) {
848  typename number_traits<typename linalg_traits<M>::value_type>
849  ::magnitude_type res(0);
850  for (size_type i = 0; i < mat_nrows(m); ++i)
851  res = std::max(res, vect_norminf(mat_const_row(m,i)));
852  return res;
853  }
854 
855  template <typename M>
856  typename number_traits<typename linalg_traits<M>::value_type>
857  ::magnitude_type
858  mat_maxnorm(const M &m, col_major) {
859  typename number_traits<typename linalg_traits<M>::value_type>
860  ::magnitude_type res(0);
861  for (size_type i = 0; i < mat_ncols(m); ++i)
862  res = std::max(res, vect_norminf(mat_const_col(m,i)));
863  return res;
864  }
865  ///@endcond
866  /** max-norm of a matrix. */
867  template <typename M>
868  typename number_traits<typename linalg_traits<M>::value_type>
869  ::magnitude_type
870  mat_maxnorm(const M &m) {
871  return mat_maxnorm(m,
872  typename principal_orientation_type<typename
873  linalg_traits<M>::sub_orientation>::potype());
874  }
875 
876  /* ******************************************************************** */
877  /* Clean */
878  /* ******************************************************************** */
879  /** Clean a vector or matrix (replace near-zero entries with zeroes). */
880 
881  template <typename L> inline void clean(L &l, double threshold);
882 
883  ///@cond DOXY_SHOW_ALL_FUNCTIONS
884 
885  template <typename L, typename T>
886  void clean(L &l, double threshold, abstract_dense, T) {
887  typedef typename number_traits<T>::magnitude_type R;
888  auto it = vect_begin(l), ite = vect_end(l);
889  for (; it != ite; ++it)
890  if (gmm::abs(*it) < R(threshold)) *it = T(0);
891  }
892 
893  template <typename L, typename T>
894  void clean(L &l, double threshold, abstract_skyline, T)
895  { gmm::clean(l, threshold, abstract_dense(), T()); }
896 
897  template <typename L, typename T>
898  void clean(L &l, double threshold, abstract_sparse, T) {
899  typedef typename number_traits<T>::magnitude_type R;
900  auto it = vect_begin(l), ite = vect_end(l);
901  std::vector<size_type> ind;
902  for (; it != ite; ++it)
903  if (gmm::abs(*it) < R(threshold)) ind.push_back(it.index());
904  for (size_type i = 0; i < ind.size(); ++i) l[ind[i]] = T(0);
905  }
906 
907  template <typename L, typename T>
908  void clean(L &l, double threshold, abstract_dense, std::complex<T>) {
909  auto it = vect_begin(l), ite = vect_end(l);
910  for (; it != ite; ++it){
911  if (gmm::abs((*it).real()) < T(threshold))
912  *it = std::complex<T>(T(0), (*it).imag());
913  if (gmm::abs((*it).imag()) < T(threshold))
914  *it = std::complex<T>((*it).real(), T(0));
915  }
916  }
917 
918  template <typename L, typename T>
919  void clean(L &l, double threshold, abstract_skyline, std::complex<T>)
920  { gmm::clean(l, threshold, abstract_dense(), std::complex<T>()); }
921 
922  template <typename L, typename T>
923  void clean(L &l, double threshold, abstract_sparse, std::complex<T>) {
924  auto it = vect_begin(l), ite = vect_end(l);
925  std::vector<size_type> ind;
926  for (; it != ite; ++it) {
927  bool r = (gmm::abs((*it).real()) < T(threshold));
928  bool i = (gmm::abs((*it).imag()) < T(threshold));
929  if (r && i) ind.push_back(it.index());
930  else if (r) *it = std::complex<T>(T(0), (*it).imag());
931  else if (i) *it = std::complex<T>((*it).real(), T(0));
932  }
933  for (size_type i = 0; i < ind.size(); ++i)
934  l[ind[i]] = std::complex<T>(T(0),T(0));
935  }
936 
937  template <typename L> inline void clean(L &l, double threshold,
938  abstract_vector) {
939  gmm::clean(l, threshold, typename linalg_traits<L>::storage_type(),
940  typename linalg_traits<L>::value_type());
941  }
942 
943  template <typename L> inline void clean(const L &l, double threshold);
944 
945  template <typename L> void clean(L &l, double threshold, row_major) {
946  for (size_type i = 0; i < mat_nrows(l); ++i)
947  gmm::clean(mat_row(l, i), threshold);
948  }
949 
950  template <typename L> void clean(L &l, double threshold, col_major) {
951  for (size_type i = 0; i < mat_ncols(l); ++i)
952  gmm::clean(mat_col(l, i), threshold);
953  }
954 
955  template <typename L> inline void clean(L &l, double threshold,
956  abstract_matrix) {
957  gmm::clean(l, threshold,
958  typename principal_orientation_type<typename
959  linalg_traits<L>::sub_orientation>::potype());
960  }
961 
962  template <typename L> inline void clean(L &l, double threshold)
963  { clean(l, threshold, typename linalg_traits<L>::linalg_type()); }
964 
965  template <typename L> inline void clean(const L &l, double threshold)
966  { gmm::clean(linalg_const_cast(l), threshold); }
967 
968  /* ******************************************************************** */
969  /* Copy */
970  /* ******************************************************************** */
971  ///@endcond
972  /** Copy vectors or matrices.
973  @param l1 source vector or matrix.
974  @param l2 destination.
975  */
976  template <typename L1, typename L2> inline
977  void copy(const L1& l1, L2& l2) {
978  if ((const void *)(&l1) != (const void *)(&l2)) {
979  if (same_origin(l1,l2))
980  GMM_WARNING2("Warning : a conflict is possible in copy\n");
981 
982  copy(l1, l2, typename linalg_traits<L1>::linalg_type(),
983  typename linalg_traits<L2>::linalg_type());
984  }
985  }
986  ///@cond DOXY_SHOW_ALL_FUNCTIONS
987 
988  template <typename L1, typename L2> inline
989  void copy(const L1& l1, const L2& l2) { copy(l1, linalg_const_cast(l2)); }
990 
991  template <typename L1, typename L2> inline
992  void copy(const L1& l1, L2& l2, abstract_vector, abstract_vector) {
993  GMM_ASSERT2(vect_size(l1) == vect_size(l2), "dimensions mismatch, "
994  << vect_size(l1) << " !=" << vect_size(l2));
995  copy_vect(l1, l2, typename linalg_traits<L1>::storage_type(),
996  typename linalg_traits<L2>::storage_type());
997  }
998 
999  template <typename L1, typename L2> inline
1000  void copy(const L1& l1, L2& l2, abstract_matrix, abstract_matrix) {
1001  size_type m = mat_nrows(l1), n = mat_ncols(l1);
1002  if (!m || !n) return;
1003  GMM_ASSERT2(n==mat_ncols(l2) && m==mat_nrows(l2), "dimensions mismatch");
1004  copy_mat(l1, l2, typename linalg_traits<L1>::sub_orientation(),
1005  typename linalg_traits<L2>::sub_orientation());
1006  }
1007 
1008  template <typename V1, typename V2, typename C1, typename C2> inline
1009  void copy_vect(const V1 &v1, const V2 &v2, C1, C2)
1010  { copy_vect(v1, const_cast<V2 &>(v2), C1(), C2()); }
1011 
1012 
1013  template <typename L1, typename L2>
1014  void copy_mat_by_row(const L1& l1, L2& l2) {
1015  size_type nbr = mat_nrows(l1);
1016  for (size_type i = 0; i < nbr; ++i)
1017  copy(mat_const_row(l1, i), mat_row(l2, i));
1018  }
1019 
1020  template <typename L1, typename L2>
1021  void copy_mat_by_col(const L1 &l1, L2 &l2) {
1022  size_type nbc = mat_ncols(l1);
1023  for (size_type i = 0; i < nbc; ++i) {
1024  copy(mat_const_col(l1, i), mat_col(l2, i));
1025  }
1026  }
1027 
1028  template <typename L1, typename L2> inline
1029  void copy_mat(const L1& l1, L2& l2, row_major, row_major)
1030  { copy_mat_by_row(l1, l2); }
1031 
1032  template <typename L1, typename L2> inline
1033  void copy_mat(const L1& l1, L2& l2, row_major, row_and_col)
1034  { copy_mat_by_row(l1, l2); }
1035 
1036  template <typename L1, typename L2> inline
1037  void copy_mat(const L1& l1, L2& l2, row_and_col, row_and_col)
1038  { copy_mat_by_row(l1, l2); }
1039 
1040  template <typename L1, typename L2> inline
1041  void copy_mat(const L1& l1, L2& l2, row_and_col, row_major)
1042  { copy_mat_by_row(l1, l2); }
1043 
1044  template <typename L1, typename L2> inline
1045  void copy_mat(const L1& l1, L2& l2, col_and_row, row_major)
1046  { copy_mat_by_row(l1, l2); }
1047 
1048  template <typename L1, typename L2> inline
1049  void copy_mat(const L1& l1, L2& l2, row_major, col_and_row)
1050  { copy_mat_by_row(l1, l2); }
1051 
1052  template <typename L1, typename L2> inline
1053  void copy_mat(const L1& l1, L2& l2, col_and_row, row_and_col)
1054  { copy_mat_by_row(l1, l2); }
1055 
1056  template <typename L1, typename L2> inline
1057  void copy_mat(const L1& l1, L2& l2, row_and_col, col_and_row)
1058  { copy_mat_by_row(l1, l2); }
1059 
1060  template <typename L1, typename L2> inline
1061  void copy_mat(const L1& l1, L2& l2, col_major, col_major)
1062  { copy_mat_by_col(l1, l2); }
1063 
1064  template <typename L1, typename L2> inline
1065  void copy_mat(const L1& l1, L2& l2, col_major, col_and_row)
1066  { copy_mat_by_col(l1, l2); }
1067 
1068  template <typename L1, typename L2> inline
1069  void copy_mat(const L1& l1, L2& l2, col_major, row_and_col)
1070  { copy_mat_by_col(l1, l2); }
1071 
1072  template <typename L1, typename L2> inline
1073  void copy_mat(const L1& l1, L2& l2, row_and_col, col_major)
1074  { copy_mat_by_col(l1, l2); }
1075 
1076  template <typename L1, typename L2> inline
1077  void copy_mat(const L1& l1, L2& l2, col_and_row, col_major)
1078  { copy_mat_by_col(l1, l2); }
1079 
1080  template <typename L1, typename L2> inline
1081  void copy_mat(const L1& l1, L2& l2, col_and_row, col_and_row)
1082  { copy_mat_by_col(l1, l2); }
1083 
1084  template <typename L1, typename L2> inline
1085  void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i) {
1086  copy_mat_mixed_rc(l1, l2, i, typename linalg_traits<L1>::storage_type());
1087  }
1088 
1089  template <typename L1, typename L2>
1090  void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_sparse) {
1091  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1092  for (; it != ite; ++it)
1093  l2(i, it.index()) = *it;
1094  }
1095 
1096  template <typename L1, typename L2>
1097  void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_skyline) {
1098  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1099  for (; it != ite; ++it)
1100  l2(i, it.index()) = *it;
1101  }
1102 
1103  template <typename L1, typename L2>
1104  void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_dense) {
1105  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1106  for (size_type j = 0; it != ite; ++it, ++j) l2(i, j) = *it;
1107  }
1108 
1109  template <typename L1, typename L2> inline
1110  void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i) {
1111  copy_mat_mixed_cr(l1, l2, i, typename linalg_traits<L1>::storage_type());
1112  }
1113 
1114  template <typename L1, typename L2>
1115  void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_sparse) {
1116  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1117  for (; it != ite; ++it) l2(it.index(), i) = *it;
1118  }
1119 
1120  template <typename L1, typename L2>
1121  void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_skyline) {
1122  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1123  for (; it != ite; ++it) l2(it.index(), i) = *it;
1124  }
1125 
1126  template <typename L1, typename L2>
1127  void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_dense) {
1128  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1129  for (size_type j = 0; it != ite; ++it, ++j) l2(j, i) = *it;
1130  }
1131 
1132  template <typename L1, typename L2>
1133  void copy_mat(const L1& l1, L2& l2, row_major, col_major) {
1134  clear(l2);
1135  size_type nbr = mat_nrows(l1);
1136  for (size_type i = 0; i < nbr; ++i)
1137  copy_mat_mixed_rc(mat_const_row(l1, i), l2, i);
1138  }
1139 
1140  template <typename L1, typename L2>
1141  void copy_mat(const L1& l1, L2& l2, col_major, row_major) {
1142  clear(l2);
1143  size_type nbc = mat_ncols(l1);
1144  for (size_type i = 0; i < nbc; ++i)
1145  copy_mat_mixed_cr(mat_const_col(l1, i), l2, i);
1146  }
1147 
1148  template <typename L1, typename L2> inline
1149  void copy_vect(const L1 &l1, L2 &l2, abstract_dense, abstract_dense) {
1150  std::copy(vect_const_begin(l1), vect_const_end(l1), vect_begin(l2));
1151  }
1152 
1153  template <typename L1, typename L2> inline // to be optimised ?
1154  void copy_vect(const L1 &l1, L2 &l2, abstract_skyline, abstract_skyline) {
1155  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1156  while (it1 != ite1 && *it1 == typename linalg_traits<L1>::value_type(0))
1157  ++it1;
1158 
1159  if (ite1 - it1 > 0) {
1160  clear(l2);
1161  auto it2 = vect_begin(l2), ite2 = vect_end(l2);
1162  while (*(ite1-1) == typename linalg_traits<L1>::value_type(0)) ite1--;
1163 
1164  if (it2 == ite2) {
1165  l2[it1.index()] = *it1; ++it1;
1166  l2[ite1.index()-1] = *(ite1-1); --ite1;
1167  if (it1 < ite1)
1168  { it2 = vect_begin(l2); ++it2; std::copy(it1, ite1, it2); }
1169  }
1170  else {
1171  ptrdiff_t m = it1.index() - it2.index();
1172  if (m >= 0 && ite1.index() <= ite2.index())
1173  std::copy(it1, ite1, it2 + m);
1174  else {
1175  if (m < 0) l2[it1.index()] = *it1;
1176  if (ite1.index() > ite2.index()) l2[ite1.index()-1] = *(ite1-1);
1177  it2 = vect_begin(l2); ite2 = vect_end(l2);
1178  m = it1.index() - it2.index();
1179  std::copy(it1, ite1, it2 + m);
1180  }
1181  }
1182  }
1183  }
1184 
1185  template <typename L1, typename L2>
1186  void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_dense) {
1187  clear(l2);
1188  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1189  for (; it != ite; ++it) { l2[it.index()] = *it; }
1190  }
1191 
1192  template <typename L1, typename L2>
1193  void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_skyline) {
1194  clear(l2);
1195  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1196  for (; it != ite; ++it) l2[it.index()] = *it;
1197  }
1198 
1199  template <typename L1, typename L2>
1200  void copy_vect(const L1& l1, L2& l2, abstract_skyline, abstract_dense) {
1201  typedef typename linalg_traits<L1>::value_type T;
1202  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1203  if (it == ite)
1204  gmm::clear(l2);
1205  else {
1206  auto it2 = vect_begin(l2), ite2 = vect_end(l2);
1207 
1208  size_type i = it.index(), j;
1209  for (j = 0; j < i; ++j, ++it2) *it2 = T(0);
1210  for (; it != ite; ++it, ++it2) *it2 = *it;
1211  for (; it2 != ite2; ++it2) *it2 = T(0);
1212  }
1213  }
1214 
1215  template <typename L1, typename L2>
1216  void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_sparse) {
1217  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1218  clear(l2);
1219  // cout << "copy " << l1 << " of size " << vect_size(l1) << " nnz = " << nnz(l1) << endl;
1220  for (; it != ite; ++it) {
1221  // cout << "*it = " << *it << endl;
1222  // cout << "it.index() = " << it.index() << endl;
1223  if (*it != (typename linalg_traits<L1>::value_type)(0))
1224  l2[it.index()] = *it;
1225  }
1226  }
1227 
1228  template <typename L1, typename L2>
1229  void copy_vect(const L1& l1, L2& l2, abstract_dense, abstract_sparse) {
1230  clear(l2);
1231  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1232  for (size_type i = 0; it != ite; ++it, ++i)
1233  if (*it != (typename linalg_traits<L1>::value_type)(0))
1234  l2[i] = *it;
1235  }
1236 
1237  template <typename L1, typename L2> // to be optimised ...
1238  void copy_vect(const L1& l1, L2& l2, abstract_dense, abstract_skyline) {
1239  clear(l2);
1240  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1241  for (size_type i = 0; it != ite; ++it, ++i)
1242  if (*it != (typename linalg_traits<L1>::value_type)(0))
1243  l2[i] = *it;
1244  }
1245 
1246 
1247  template <typename L1, typename L2>
1248  void copy_vect(const L1& l1, L2& l2, abstract_skyline, abstract_sparse) {
1249  clear(l2);
1250  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1251  for (; it != ite; ++it)
1252  if (*it != (typename linalg_traits<L1>::value_type)(0))
1253  l2[it.index()] = *it;
1254  }
1255 
1256  /* ******************************************************************** */
1257  /* Matrix and vector addition */
1258  /* algorithms are built in order to avoid some conflicts with */
1259  /* repeated arguments or with overlapping part of a same object. */
1260  /* In the latter case, conflicts are still possible. */
1261  /* ******************************************************************** */
1262  ///@endcond
1263  /** Add two vectors or matrices
1264  @param l1
1265  @param l2 contains on output, l2+l1.
1266  */
1267  template <typename L1, typename L2> inline
1268  void add(const L1& l1, L2& l2) {
1269  add_spec(l1, l2, typename linalg_traits<L2>::linalg_type());
1270  }
1271  ///@cond
1272 
1273  template <typename L1, typename L2> inline
1274  void add(const L1& l1, const L2& l2) { add(l1, linalg_const_cast(l2)); }
1275 
1276  template <typename L1, typename L2> inline
1277  void add_spec(const L1& l1, L2& l2, abstract_vector) {
1278  GMM_ASSERT2(vect_size(l1) == vect_size(l2), "dimensions mismatch, "
1279  << vect_size(l1) << " !=" << vect_size(l2));
1280  add(l1, l2, typename linalg_traits<L1>::storage_type(),
1281  typename linalg_traits<L2>::storage_type());
1282  }
1283 
1284  template <typename L1, typename L2> inline
1285  void add_spec(const L1& l1, L2& l2, abstract_matrix) {
1286  GMM_ASSERT2(mat_nrows(l1)==mat_nrows(l2) && mat_ncols(l1)==mat_ncols(l2),
1287  "dimensions mismatch l1 is " << mat_nrows(l1) << "x"
1288  << mat_ncols(l1) << " and l2 is " << mat_nrows(l2)
1289  << "x" << mat_ncols(l2));
1290  add(l1, l2, typename linalg_traits<L1>::sub_orientation(),
1291  typename linalg_traits<L2>::sub_orientation());
1292  }
1293 
1294  template <typename L1, typename L2>
1295  void add(const L1& l1, L2& l2, row_major, row_major) {
1296  auto it1 = mat_row_begin(l1), ite = mat_row_end(l1);
1297  auto it2 = mat_row_begin(l2);
1298  for ( ; it1 != ite; ++it1, ++it2)
1299  add(linalg_traits<L1>::row(it1), linalg_traits<L2>::row(it2));
1300  }
1301 
1302  template <typename L1, typename L2>
1303  void add(const L1& l1, L2& l2, col_major, col_major) {
1304  auto it1 = mat_col_const_begin(l1), ite = mat_col_const_end(l1);
1305  typename linalg_traits<L2>::col_iterator it2 = mat_col_begin(l2);
1306  for ( ; it1 != ite; ++it1, ++it2)
1307  add(linalg_traits<L1>::col(it1), linalg_traits<L2>::col(it2));
1308  }
1309 
1310  template <typename L1, typename L2> inline
1311  void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i) {
1312  add_mat_mixed_rc(l1, l2, i, typename linalg_traits<L1>::storage_type());
1313  }
1314 
1315  template <typename L1, typename L2>
1316  void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_sparse) {
1317  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1318  for (; it != ite; ++it) l2(i, it.index()) += *it;
1319  }
1320 
1321  template <typename L1, typename L2>
1322  void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_skyline) {
1323  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1324  for (; it != ite; ++it) l2(i, it.index()) += *it;
1325  }
1326 
1327  template <typename L1, typename L2>
1328  void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_dense) {
1329  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1330  for (size_type j = 0; it != ite; ++it, ++j) l2(i, j) += *it;
1331  }
1332 
1333  template <typename L1, typename L2> inline
1334  void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i) {
1335  add_mat_mixed_cr(l1, l2, i, typename linalg_traits<L1>::storage_type());
1336  }
1337 
1338  template <typename L1, typename L2>
1339  void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_sparse) {
1340  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1341  for (; it != ite; ++it) l2(it.index(), i) += *it;
1342  }
1343 
1344  template <typename L1, typename L2>
1345  void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_skyline) {
1346  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1347  for (; it != ite; ++it) l2(it.index(), i) += *it;
1348  }
1349 
1350  template <typename L1, typename L2>
1351  void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_dense) {
1352  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1353  for (size_type j = 0; it != ite; ++it, ++j) l2(j, i) += *it;
1354  }
1355 
1356  template <typename L1, typename L2>
1357  void add(const L1& l1, L2& l2, row_major, col_major) {
1358  size_type nbr = mat_nrows(l1);
1359  for (size_type i = 0; i < nbr; ++i)
1360  add_mat_mixed_rc(mat_const_row(l1, i), l2, i);
1361  }
1362 
1363  template <typename L1, typename L2>
1364  void add(const L1& l1, L2& l2, col_major, row_major) {
1365  size_type nbc = mat_ncols(l1);
1366  for (size_type i = 0; i < nbc; ++i)
1367  add_mat_mixed_cr(mat_const_col(l1, i), l2, i);
1368  }
1369 
1370  template <typename L1, typename L2> inline
1371  void add(const L1& l1, L2& l2, row_and_col, row_major)
1372  { add(l1, l2, row_major(), row_major()); }
1373 
1374  template <typename L1, typename L2> inline
1375  void add(const L1& l1, L2& l2, row_and_col, row_and_col)
1376  { add(l1, l2, row_major(), row_major()); }
1377 
1378  template <typename L1, typename L2> inline
1379  void add(const L1& l1, L2& l2, row_and_col, col_and_row)
1380  { add(l1, l2, row_major(), row_major()); }
1381 
1382  template <typename L1, typename L2> inline
1383  void add(const L1& l1, L2& l2, col_and_row, row_and_col)
1384  { add(l1, l2, row_major(), row_major()); }
1385 
1386  template <typename L1, typename L2> inline
1387  void add(const L1& l1, L2& l2, row_major, row_and_col)
1388  { add(l1, l2, row_major(), row_major()); }
1389 
1390  template <typename L1, typename L2> inline
1391  void add(const L1& l1, L2& l2, col_and_row, row_major)
1392  { add(l1, l2, row_major(), row_major()); }
1393 
1394  template <typename L1, typename L2> inline
1395  void add(const L1& l1, L2& l2, row_major, col_and_row)
1396  { add(l1, l2, row_major(), row_major()); }
1397 
1398  template <typename L1, typename L2> inline
1399  void add(const L1& l1, L2& l2, row_and_col, col_major)
1400  { add(l1, l2, col_major(), col_major()); }
1401 
1402  template <typename L1, typename L2> inline
1403  void add(const L1& l1, L2& l2, col_major, row_and_col)
1404  { add(l1, l2, col_major(), col_major()); }
1405 
1406  template <typename L1, typename L2> inline
1407  void add(const L1& l1, L2& l2, col_and_row, col_major)
1408  { add(l1, l2, col_major(), col_major()); }
1409 
1410  template <typename L1, typename L2> inline
1411  void add(const L1& l1, L2& l2, col_and_row, col_and_row)
1412  { add(l1, l2, col_major(), col_major()); }
1413 
1414  template <typename L1, typename L2> inline
1415  void add(const L1& l1, L2& l2, col_major, col_and_row)
1416  { add(l1, l2, col_major(), col_major()); }
1417 
1418  ///@endcond
1419  /** Addition of two vectors/matrices
1420  @param l1
1421  @param l2
1422  @param l3 contains l1+l2 on output
1423  */
1424  template <typename L1, typename L2, typename L3> inline
1425  void add(const L1& l1, const L2& l2, L3& l3) {
1426  add_spec(l1, l2, l3, typename linalg_traits<L2>::linalg_type());
1427  }
1428  ///@cond DOXY_SHOW_ALL_FUNCTIONS
1429 
1430  template <typename L1, typename L2, typename L3> inline
1431  void add(const L1& l1, const L2& l2, const L3& l3)
1432  { add(l1, l2, linalg_const_cast(l3)); }
1433 
1434  template <typename L1, typename L2, typename L3> inline
1435  void add_spec(const L1& l1, const L2& l2, L3& l3, abstract_matrix)
1436  { copy(l2, l3); add(l1, l3); }
1437 
1438  template <typename L1, typename L2, typename L3> inline
1439  void add_spec(const L1& l1, const L2& l2, L3& l3, abstract_vector) {
1440  GMM_ASSERT2(vect_size(l1) == vect_size(l2) &&
1441  vect_size(l1) == vect_size(l3), "dimensions mismatch");
1442  if ((const void *)(&l1) == (const void *)(&l3))
1443  add(l2, l3);
1444  else if ((const void *)(&l2) == (const void *)(&l3))
1445  add(l1, l3);
1446  else
1447  add(l1, l2, l3, typename linalg_traits<L1>::storage_type(),
1448  typename linalg_traits<L2>::storage_type(),
1449  typename linalg_traits<L3>::storage_type());
1450  }
1451 
1452  template <typename IT1, typename IT2, typename IT3>
1453  void add_full_(IT1 it1, IT2 it2, IT3 it3, IT3 ite) {
1454  for (; it3 != ite; ++it3, ++it2, ++it1) *it3 = *it1 + *it2;
1455  }
1456 
1457  template <typename IT1, typename IT2, typename IT3>
1458  void add_almost_full_(IT1 it1, IT1 ite1, IT2 it2, IT3 it3, IT3 ite3) {
1459  IT3 it = it3;
1460  for (; it != ite3; ++it, ++it2) *it = *it2;
1461  for (; it1 != ite1; ++it1)
1462  *(it3 + it1.index()) += *it1;
1463  }
1464 
1465  template <typename IT1, typename IT2, typename IT3>
1466  void add_to_full_(IT1 it1, IT1 ite1, IT2 it2, IT2 ite2,
1467  IT3 it3, IT3 ite3) {
1468  typedef typename std::iterator_traits<IT3>::value_type T;
1469  IT3 it = it3;
1470  for (; it != ite3; ++it) *it = T(0);
1471  for (; it1 != ite1; ++it1) *(it3 + it1.index()) = *it1;
1472  for (; it2 != ite2; ++it2) *(it3 + it2.index()) += *it2;
1473  }
1474 
1475  template <typename L1, typename L2, typename L3> inline
1476  void add(const L1& l1, const L2& l2, L3& l3,
1477  abstract_dense, abstract_dense, abstract_dense) {
1478  add_full_(vect_const_begin(l1), vect_const_begin(l2),
1479  vect_begin(l3), vect_end(l3));
1480  }
1481 
1482  // generic function for add(v1, v2, v3).
1483  // Need to be specialized to optimize particular additions.
1484  template <typename L1, typename L2, typename L3,
1485  typename ST1, typename ST2, typename ST3>
1486  inline void add(const L1& l1, const L2& l2, L3& l3, ST1, ST2, ST3)
1487  { copy(l2, l3); add(l1, l3, ST1(), ST3()); }
1488 
1489  template <typename L1, typename L2, typename L3> inline
1490  void add(const L1& l1, const L2& l2, L3& l3,
1491  abstract_sparse, abstract_dense, abstract_dense) {
1492  add_almost_full_(vect_const_begin(l1), vect_const_end(l1),
1493  vect_const_begin(l2), vect_begin(l3), vect_end(l3));
1494  }
1495 
1496  template <typename L1, typename L2, typename L3> inline
1497  void add(const L1& l1, const L2& l2, L3& l3,
1498  abstract_dense, abstract_sparse, abstract_dense)
1499  { add(l2, l1, l3, abstract_sparse(), abstract_dense(), abstract_dense()); }
1500 
1501  template <typename L1, typename L2, typename L3> inline
1502  void add(const L1& l1, const L2& l2, L3& l3,
1503  abstract_sparse, abstract_sparse, abstract_dense) {
1504  add_to_full_(vect_const_begin(l1), vect_const_end(l1),
1505  vect_const_begin(l2), vect_const_end(l2),
1506  vect_begin(l3), vect_end(l3));
1507  }
1508 
1509 
1510  template <typename L1, typename L2, typename L3>
1511  void add_spspsp(const L1& l1, const L2& l2, L3& l3, linalg_true) {
1512  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1513  auto it2 = vect_const_begin(l2), ite2 = vect_const_end(l2);
1514  clear(l3);
1515  while (it1 != ite1 && it2 != ite2) {
1516  ptrdiff_t d = it1.index() - it2.index();
1517  if (d < 0)
1518  { l3[it1.index()] += *it1; ++it1; }
1519  else if (d > 0)
1520  { l3[it2.index()] += *it2; ++it2; }
1521  else
1522  { l3[it1.index()] = *it1 + *it2; ++it1; ++it2; }
1523  }
1524  for (; it1 != ite1; ++it1) l3[it1.index()] += *it1;
1525  for (; it2 != ite2; ++it2) l3[it2.index()] += *it2;
1526  }
1527 
1528  template <typename L1, typename L2, typename L3>
1529  void add_spspsp(const L1& l1, const L2& l2, L3& l3, linalg_false)
1530  { copy(l2, l3); add(l2, l3); }
1531 
1532  template <typename L1, typename L2, typename L3>
1533  void add(const L1& l1, const L2& l2, L3& l3,
1534  abstract_sparse, abstract_sparse, abstract_sparse) {
1535  add_spspsp(l1, l2, l3, typename linalg_and<typename
1536  linalg_traits<L1>::index_sorted,
1537  typename linalg_traits<L2>::index_sorted>::bool_type());
1538  }
1539 
1540  template <typename L1, typename L2>
1541  void add(const L1& l1, L2& l2, abstract_dense, abstract_dense) {
1542  auto it1 = vect_const_begin(l1);
1543  auto it2 = vect_begin(l2), ite = vect_end(l2);
1544  for (; it2 != ite; ++it2, ++it1) *it2 += *it1;
1545  }
1546 
1547  template <typename L1, typename L2>
1548  void add(const L1& l1, L2& l2, abstract_dense, abstract_skyline) {
1549  typedef typename linalg_traits<L1>::value_type T;
1550 
1551  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1552  size_type i1 = 0, ie1 = vect_size(l1);
1553  while (it1 != ite1 && *it1 == T(0)) { ++it1; ++i1; }
1554  if (it1 != ite1) {
1555  auto it2 = vect_begin(l2), ite2 = vect_end(l2);
1556  while (ie1 && *(ite1-1) == T(0)) { ite1--; --ie1; }
1557 
1558  if (it2 == ite2 || i1 < it2.index()) {
1559  l2[i1] = *it1; ++i1; ++it1;
1560  if (it1 == ite1) return;
1561  it2 = vect_begin(l2); ite2 = vect_end(l2);
1562  }
1563  if (ie1 > ite2.index()) {
1564  --ite1; l2[ie1 - 1] = *ite1;
1565  it2 = vect_begin(l2);
1566  }
1567  it2 += i1 - it2.index();
1568  for (; it1 != ite1; ++it1, ++it2) { *it2 += *it1; }
1569  }
1570  }
1571 
1572 
1573  template <typename L1, typename L2>
1574  void add(const L1& l1, L2& l2, abstract_skyline, abstract_dense) {
1575  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1576  if (it1 != ite1) {
1577  auto it2 = vect_begin(l2);
1578  it2 += it1.index();
1579  for (; it1 != ite1; ++it2, ++it1) *it2 += *it1;
1580  }
1581  }
1582 
1583 
1584  template <typename L1, typename L2>
1585  void add(const L1& l1, L2& l2, abstract_sparse, abstract_dense) {
1586  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1587  for (; it1 != ite1; ++it1) l2[it1.index()] += *it1;
1588  }
1589 
1590  template <typename L1, typename L2>
1591  void add(const L1& l1, L2& l2, abstract_sparse, abstract_sparse) {
1592  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1593  for (; it1 != ite1; ++it1) l2[it1.index()] += *it1;
1594  }
1595 
1596  template <typename L1, typename L2>
1597  void add(const L1& l1, L2& l2, abstract_sparse, abstract_skyline) {
1598  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1599  for (; it1 != ite1; ++it1) l2[it1.index()] += *it1;
1600  }
1601 
1602 
1603  template <typename L1, typename L2>
1604  void add(const L1& l1, L2& l2, abstract_skyline, abstract_sparse) {
1605  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1606  for (; it1 != ite1; ++it1)
1607  if (*it1 != typename linalg_traits<L1>::value_type(0))
1608  l2[it1.index()] += *it1;
1609  }
1610 
1611  template <typename L1, typename L2>
1612  void add(const L1& l1, L2& l2, abstract_skyline, abstract_skyline) {
1613  typedef typename linalg_traits<L1>::value_type T1;
1614  typedef typename linalg_traits<L2>::value_type T2;
1615 
1616  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1617 
1618  while (it1 != ite1 && *it1 == T1(0)) ++it1;
1619  if (ite1 != it1) {
1620  auto it2 = vect_begin(l2), ite2 = vect_end(l2);
1621  while (*(ite1-1) == T1(0)) ite1--;
1622  if (it2 == ite2 || it1.index() < it2.index()) {
1623  l2[it1.index()] = T2(0);
1624  it2 = vect_begin(l2); ite2 = vect_end(l2);
1625  }
1626  if (ite1.index() > ite2.index()) {
1627  l2[ite1.index() - 1] = T2(0);
1628  it2 = vect_begin(l2);
1629  }
1630  it2 += it1.index() - it2.index();
1631  for (; it1 != ite1; ++it1, ++it2) *it2 += *it1;
1632  }
1633  }
1634 
1635  template <typename L1, typename L2>
1636  void add(const L1& l1, L2& l2, abstract_dense, abstract_sparse) {
1637  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1638  for (size_type i = 0; it1 != ite1; ++it1, ++i)
1639  if (*it1 != typename linalg_traits<L1>::value_type(0)) l2[i] += *it1;
1640  }
1641 
1642  /* ******************************************************************** */
1643  /* Matrix-vector mult */
1644  /* ******************************************************************** */
1645  ///@endcond
1646  /** matrix-vector or matrix-matrix product.
1647  @param l1 a matrix.
1648  @param l2 a vector or matrix.
1649  @param l3 the product l1*l2.
1650  */
1651  template <typename L1, typename L2, typename L3> inline
1652  void mult(const L1& l1, const L2& l2, L3& l3) {
1653  mult_dispatch(l1, l2, l3, typename linalg_traits<L2>::linalg_type());
1654  }
1655  ///@cond DOXY_SHOW_ALL_FUNCTIONS
1656 
1657  template <typename L1, typename L2, typename L3> inline
1658  void mult(const L1& l1, const L2& l2, const L3& l3)
1659  { mult(l1, l2, linalg_const_cast(l3)); }
1660 
1661  template <typename L1, typename L2, typename L3> inline
1662  void mult_dispatch(const L1& l1, const L2& l2, L3& l3, abstract_vector) {
1663  size_type m = mat_nrows(l1), n = mat_ncols(l1);
1664  if (!m || !n) { gmm::clear(l3); return; }
1665  GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l3), "dimensions mismatch");
1666  if (!same_origin(l2, l3))
1667  mult_spec(l1, l2, l3, typename principal_orientation_type<typename
1668  linalg_traits<L1>::sub_orientation>::potype());
1669  else {
1670  GMM_WARNING2("Warning, A temporary is used for mult\n");
1671  typename temporary_vector<L3>::vector_type temp(vect_size(l3));
1672  mult_spec(l1, l2, temp, typename principal_orientation_type<typename
1673  linalg_traits<L1>::sub_orientation>::potype());
1674  copy(temp, l3);
1675  }
1676  }
1677 
1678  template <typename L1, typename L2, typename L3>
1679  void mult_by_row(const L1& l1, const L2& l2, L3& l3, abstract_sparse) {
1680  typedef typename linalg_traits<L3>::value_type T;
1681  clear(l3);
1682  size_type nr = mat_nrows(l1);
1683  for (size_type i = 0; i < nr; ++i) {
1684  T aux = vect_sp(mat_const_row(l1, i), l2);
1685  if (aux != T(0)) l3[i] = aux;
1686  }
1687  }
1688 
1689  template <typename L1, typename L2, typename L3>
1690  void mult_by_row(const L1& l1, const L2& l2, L3& l3, abstract_skyline) {
1691  typedef typename linalg_traits<L3>::value_type T;
1692  clear(l3);
1693  size_type nr = mat_nrows(l1);
1694  for (size_type i = 0; i < nr; ++i) {
1695  T aux = vect_sp(mat_const_row(l1, i), l2);
1696  if (aux != T(0)) l3[i] = aux;
1697  }
1698  }
1699 
1700  template <typename L1, typename L2, typename L3>
1701  void mult_by_row(const L1& l1, const L2& l2, L3& l3, abstract_dense) {
1702  typename linalg_traits<L3>::iterator it=vect_begin(l3), ite=vect_end(l3);
1703  auto itr = mat_row_const_begin(l1);
1704  for (; it != ite; ++it, ++itr)
1705  *it = vect_sp(linalg_traits<L1>::row(itr), l2,
1706  typename linalg_traits<L1>::storage_type(),
1707  typename linalg_traits<L2>::storage_type());
1708  }
1709 
1710  template <typename L1, typename L2, typename L3>
1711  void mult_by_col(const L1& l1, const L2& l2, L3& l3, abstract_dense) {
1712  clear(l3);
1713  size_type nc = mat_ncols(l1);
1714  for (size_type i = 0; i < nc; ++i)
1715  add(scaled(mat_const_col(l1, i), l2[i]), l3);
1716  }
1717 
1718  template <typename L1, typename L2, typename L3>
1719  void mult_by_col(const L1& l1, const L2& l2, L3& l3, abstract_sparse) {
1720  typedef typename linalg_traits<L2>::value_type T;
1721  clear(l3);
1722  auto it = vect_const_begin(l2), ite = vect_const_end(l2);
1723  for (; it != ite; ++it)
1724  if (*it != T(0)) add(scaled(mat_const_col(l1, it.index()), *it), l3);
1725  }
1726 
1727  template <typename L1, typename L2, typename L3>
1728  void mult_by_col(const L1& l1, const L2& l2, L3& l3, abstract_skyline) {
1729  typedef typename linalg_traits<L2>::value_type T;
1730  clear(l3);
1731  auto it = vect_const_begin(l2), ite = vect_const_end(l2);
1732  for (; it != ite; ++it)
1733  if (*it != T(0)) add(scaled(mat_const_col(l1, it.index()), *it), l3);
1734  }
1735 
1736  template <typename L1, typename L2, typename L3> inline
1737  void mult_spec(const L1& l1, const L2& l2, L3& l3, row_major)
1738  { mult_by_row(l1, l2, l3, typename linalg_traits<L3>::storage_type()); }
1739 
1740  template <typename L1, typename L2, typename L3> inline
1741  void mult_spec(const L1& l1, const L2& l2, L3& l3, col_major)
1742  { mult_by_col(l1, l2, l3, typename linalg_traits<L2>::storage_type()); }
1743 
1744  template <typename L1, typename L2, typename L3> inline
1745  void mult_spec(const L1& l1, const L2& l2, L3& l3, abstract_null_type)
1746  { mult_ind(l1, l2, l3, typename linalg_traits<L1>::storage_type()); }
1747 
1748  template <typename L1, typename L2, typename L3>
1749  void mult_ind(const L1& l1, const L2& l2, L3& l3, abstract_indirect) {
1750  GMM_ASSERT1(false, "gmm::mult(m, ., .) undefined for this kind of matrix");
1751  }
1752 
1753  template <typename L1, typename L2, typename L3, typename L4> inline
1754  void mult(const L1& l1, const L2& l2, const L3& l3, L4& l4) {
1755  size_type m = mat_nrows(l1), n = mat_ncols(l1);
1756  copy(l3, l4);
1757  if (!m || !n) { gmm::copy(l3, l4); return; }
1758  GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l4), "dimensions mismatch");
1759  if (!same_origin(l2, l4)) {
1760  mult_add_spec(l1, l2, l4, typename principal_orientation_type<typename
1761  linalg_traits<L1>::sub_orientation>::potype());
1762  }
1763  else {
1764  GMM_WARNING2("Warning, A temporary is used for mult\n");
1765  typename temporary_vector<L2>::vector_type temp(vect_size(l2));
1766  copy(l2, temp);
1767  mult_add_spec(l1,temp, l4, typename principal_orientation_type<typename
1768  linalg_traits<L1>::sub_orientation>::potype());
1769  }
1770  }
1771 
1772  template <typename L1, typename L2, typename L3, typename L4> inline
1773  void mult(const L1& l1, const L2& l2, const L3& l3, const L4& l4)
1774  { mult(l1, l2, l3, linalg_const_cast(l4)); }
1775 
1776  ///@endcond
1777  /** Multiply-accumulate. l3 += l1*l2; */
1778  template <typename L1, typename L2, typename L3> inline
1779  void mult_add(const L1& l1, const L2& l2, L3& l3) {
1780  size_type m = mat_nrows(l1), n = mat_ncols(l1);
1781  if (!m || !n) return;
1782  GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l3), "dimensions mismatch");
1783  if (!same_origin(l2, l3)) {
1784  mult_add_spec(l1, l2, l3, typename principal_orientation_type<typename
1785  linalg_traits<L1>::sub_orientation>::potype());
1786  }
1787  else {
1788  GMM_WARNING2("Warning, A temporary is used for mult\n");
1789  typename temporary_vector<L3>::vector_type temp(vect_size(l2));
1790  copy(l2, temp);
1791  mult_add_spec(l1,temp, l3, typename principal_orientation_type<typename
1792  linalg_traits<L1>::sub_orientation>::potype());
1793  }
1794  }
1795  ///@cond DOXY_SHOW_ALL_FUNCTIONS
1796 
1797  template <typename L1, typename L2, typename L3> inline
1798  void mult_add(const L1& l1, const L2& l2, const L3& l3)
1799  { mult_add(l1, l2, linalg_const_cast(l3)); }
1800 
1801  template <typename L1, typename L2, typename L3>
1802  void mult_add_by_row(const L1& l1, const L2& l2, L3& l3, abstract_sparse) {
1803  typedef typename linalg_traits<L3>::value_type T;
1804  size_type nr = mat_nrows(l1);
1805  for (size_type i = 0; i < nr; ++i) {
1806  T aux = vect_sp(mat_const_row(l1, i), l2);
1807  if (aux != T(0)) l3[i] += aux;
1808  }
1809  }
1810 
1811  template <typename L1, typename L2, typename L3>
1812  void mult_add_by_row(const L1& l1, const L2& l2, L3& l3, abstract_skyline) {
1813  typedef typename linalg_traits<L3>::value_type T;
1814  size_type nr = mat_nrows(l1);
1815  for (size_type i = 0; i < nr; ++i) {
1816  T aux = vect_sp(mat_const_row(l1, i), l2);
1817  if (aux != T(0)) l3[i] += aux;
1818  }
1819  }
1820 
1821  template <typename L1, typename L2, typename L3>
1822  void mult_add_by_row(const L1& l1, const L2& l2, L3& l3, abstract_dense) {
1823  auto it=vect_begin(l3), ite=vect_end(l3);
1824  auto itr = mat_row_const_begin(l1);
1825  for (; it != ite; ++it, ++itr)
1826  *it += vect_sp(linalg_traits<L1>::row(itr), l2);
1827  }
1828 
1829  template <typename L1, typename L2, typename L3>
1830  void mult_add_by_col(const L1& l1, const L2& l2, L3& l3, abstract_dense) {
1831  size_type nc = mat_ncols(l1);
1832  for (size_type i = 0; i < nc; ++i)
1833  add(scaled(mat_const_col(l1, i), l2[i]), l3);
1834  }
1835 
1836  template <typename L1, typename L2, typename L3>
1837  void mult_add_by_col(const L1& l1, const L2& l2, L3& l3, abstract_sparse) {
1838  auto it = vect_const_begin(l2), ite = vect_const_end(l2);
1839  for (; it != ite; ++it)
1840  if (*it != typename linalg_traits<L2>::value_type(0))
1841  add(scaled(mat_const_col(l1, it.index()), *it), l3);
1842  }
1843 
1844  template <typename L1, typename L2, typename L3>
1845  void mult_add_by_col(const L1& l1, const L2& l2, L3& l3, abstract_skyline) {
1846  auto it = vect_const_begin(l2), ite = vect_const_end(l2);
1847  for (; it != ite; ++it)
1848  if (*it != typename linalg_traits<L2>::value_type(0))
1849  add(scaled(mat_const_col(l1, it.index()), *it), l3);
1850  }
1851 
1852  template <typename L1, typename L2, typename L3> inline
1853  void mult_add_spec(const L1& l1, const L2& l2, L3& l3, row_major)
1854  { mult_add_by_row(l1, l2, l3, typename linalg_traits<L3>::storage_type()); }
1855 
1856  template <typename L1, typename L2, typename L3> inline
1857  void mult_add_spec(const L1& l1, const L2& l2, L3& l3, col_major)
1858  { mult_add_by_col(l1, l2, l3, typename linalg_traits<L2>::storage_type()); }
1859 
1860  template <typename L1, typename L2, typename L3> inline
1861  void mult_add_spec(const L1& l1, const L2& l2, L3& l3, abstract_null_type)
1862  { mult_ind(l1, l2, l3, typename linalg_traits<L1>::storage_type()); }
1863 
1864  template <typename L1, typename L2, typename L3>
1865  void transposed_mult(const L1& l1, const L2& l2, const L3& l3)
1866  { mult(gmm::transposed(l1), l2, l3); }
1867 
1868 
1869  /* ******************************************************************** */
1870  /* Matrix-matrix mult */
1871  /* ******************************************************************** */
1872 
1873 
1874  struct g_mult {}; // generic mult, less optimized
1875  struct c_mult {}; // col x col -> col mult
1876  struct r_mult {}; // row x row -> row mult
1877  struct rcmult {}; // row x col -> col mult
1878  struct crmult {}; // col x row -> row mult
1879 
1880 
1881  template<typename SO1, typename SO2, typename SO3> struct mult_t;
1882  #define DEFMU__ template<> struct mult_t
1883  DEFMU__<row_major , row_major , row_major > { typedef r_mult t; };
1884  DEFMU__<row_major , row_major , col_major > { typedef g_mult t; };
1885  DEFMU__<row_major , row_major , col_and_row> { typedef r_mult t; };
1886  DEFMU__<row_major , row_major , row_and_col> { typedef r_mult t; };
1887  DEFMU__<row_major , col_major , row_major > { typedef rcmult t; };
1888  DEFMU__<row_major , col_major , col_major > { typedef rcmult t; };
1889  DEFMU__<row_major , col_major , col_and_row> { typedef rcmult t; };
1890  DEFMU__<row_major , col_major , row_and_col> { typedef rcmult t; };
1891  DEFMU__<row_major , col_and_row, row_major > { typedef r_mult t; };
1892  DEFMU__<row_major , col_and_row, col_major > { typedef rcmult t; };
1893  DEFMU__<row_major , col_and_row, col_and_row> { typedef rcmult t; };
1894  DEFMU__<row_major , col_and_row, row_and_col> { typedef rcmult t; };
1895  DEFMU__<row_major , row_and_col, row_major > { typedef r_mult t; };
1896  DEFMU__<row_major , row_and_col, col_major > { typedef rcmult t; };
1897  DEFMU__<row_major , row_and_col, col_and_row> { typedef r_mult t; };
1898  DEFMU__<row_major , row_and_col, row_and_col> { typedef r_mult t; };
1899  DEFMU__<col_major , row_major , row_major > { typedef crmult t; };
1900  DEFMU__<col_major , row_major , col_major > { typedef g_mult t; };
1901  DEFMU__<col_major , row_major , col_and_row> { typedef crmult t; };
1902  DEFMU__<col_major , row_major , row_and_col> { typedef crmult t; };
1903  DEFMU__<col_major , col_major , row_major > { typedef g_mult t; };
1904  DEFMU__<col_major , col_major , col_major > { typedef c_mult t; };
1905  DEFMU__<col_major , col_major , col_and_row> { typedef c_mult t; };
1906  DEFMU__<col_major , col_major , row_and_col> { typedef c_mult t; };
1907  DEFMU__<col_major , col_and_row, row_major > { typedef crmult t; };
1908  DEFMU__<col_major , col_and_row, col_major > { typedef c_mult t; };
1909  DEFMU__<col_major , col_and_row, col_and_row> { typedef c_mult t; };
1910  DEFMU__<col_major , col_and_row, row_and_col> { typedef c_mult t; };
1911  DEFMU__<col_major , row_and_col, row_major > { typedef crmult t; };
1912  DEFMU__<col_major , row_and_col, col_major > { typedef c_mult t; };
1913  DEFMU__<col_major , row_and_col, col_and_row> { typedef c_mult t; };
1914  DEFMU__<col_major , row_and_col, row_and_col> { typedef c_mult t; };
1915  DEFMU__<col_and_row, row_major , row_major > { typedef r_mult t; };
1916  DEFMU__<col_and_row, row_major , col_major > { typedef c_mult t; };
1917  DEFMU__<col_and_row, row_major , col_and_row> { typedef r_mult t; };
1918  DEFMU__<col_and_row, row_major , row_and_col> { typedef r_mult t; };
1919  DEFMU__<col_and_row, col_major , row_major > { typedef rcmult t; };
1920  DEFMU__<col_and_row, col_major , col_major > { typedef c_mult t; };
1921  DEFMU__<col_and_row, col_major , col_and_row> { typedef c_mult t; };
1922  DEFMU__<col_and_row, col_major , row_and_col> { typedef c_mult t; };
1923  DEFMU__<col_and_row, col_and_row, row_major > { typedef r_mult t; };
1924  DEFMU__<col_and_row, col_and_row, col_major > { typedef c_mult t; };
1925  DEFMU__<col_and_row, col_and_row, col_and_row> { typedef c_mult t; };
1926  DEFMU__<col_and_row, col_and_row, row_and_col> { typedef c_mult t; };
1927  DEFMU__<col_and_row, row_and_col, row_major > { typedef r_mult t; };
1928  DEFMU__<col_and_row, row_and_col, col_major > { typedef c_mult t; };
1929  DEFMU__<col_and_row, row_and_col, col_and_row> { typedef c_mult t; };
1930  DEFMU__<col_and_row, row_and_col, row_and_col> { typedef r_mult t; };
1931  DEFMU__<row_and_col, row_major , row_major > { typedef r_mult t; };
1932  DEFMU__<row_and_col, row_major , col_major > { typedef c_mult t; };
1933  DEFMU__<row_and_col, row_major , col_and_row> { typedef r_mult t; };
1934  DEFMU__<row_and_col, row_major , row_and_col> { typedef r_mult t; };
1935  DEFMU__<row_and_col, col_major , row_major > { typedef rcmult t; };
1936  DEFMU__<row_and_col, col_major , col_major > { typedef c_mult t; };
1937  DEFMU__<row_and_col, col_major , col_and_row> { typedef c_mult t; };
1938  DEFMU__<row_and_col, col_major , row_and_col> { typedef c_mult t; };
1939  DEFMU__<row_and_col, col_and_row, row_major > { typedef rcmult t; };
1940  DEFMU__<row_and_col, col_and_row, col_major > { typedef rcmult t; };
1941  DEFMU__<row_and_col, col_and_row, col_and_row> { typedef rcmult t; };
1942  DEFMU__<row_and_col, col_and_row, row_and_col> { typedef rcmult t; };
1943  DEFMU__<row_and_col, row_and_col, row_major > { typedef r_mult t; };
1944  DEFMU__<row_and_col, row_and_col, col_major > { typedef c_mult t; };
1945  DEFMU__<row_and_col, row_and_col, col_and_row> { typedef r_mult t; };
1946  DEFMU__<row_and_col, row_and_col, row_and_col> { typedef r_mult t; };
1947 
1948  template <typename L1, typename L2, typename L3>
1949  void mult_dispatch(const L1& l1, const L2& l2, L3& l3, abstract_matrix) {
1950  typedef typename temporary_matrix<L3>::matrix_type temp_mat_type;
1951  size_type n = mat_ncols(l1);
1952  if (n == 0) { gmm::clear(l3); return; }
1953  GMM_ASSERT2(n == mat_nrows(l2) && mat_nrows(l1) == mat_nrows(l3) &&
1954  mat_ncols(l2) == mat_ncols(l3), "dimensions mismatch");
1955 
1956  if (same_origin(l2, l3) || same_origin(l1, l3)) {
1957  GMM_WARNING2("A temporary is used for mult");
1958  temp_mat_type temp(mat_nrows(l3), mat_ncols(l3));
1959  mult_spec(l1, l2, temp, typename mult_t<
1960  typename linalg_traits<L1>::sub_orientation,
1961  typename linalg_traits<L2>::sub_orientation,
1962  typename linalg_traits<temp_mat_type>::sub_orientation>::t());
1963  copy(temp, l3);
1964  }
1965  else
1966  mult_spec(l1, l2, l3, typename mult_t<
1967  typename linalg_traits<L1>::sub_orientation,
1968  typename linalg_traits<L2>::sub_orientation,
1969  typename linalg_traits<L3>::sub_orientation>::t());
1970  }
1971 
1972  // Completely generic but inefficient
1973 
1974  template <typename L1, typename L2, typename L3>
1975  void mult_spec(const L1& l1, const L2& l2, L3& l3, g_mult) {
1976  typedef typename linalg_traits<L3>::value_type T;
1977  GMM_WARNING2("Inefficient generic matrix-matrix mult is used");
1978  for (size_type i = 0; i < mat_nrows(l3) ; ++i)
1979  for (size_type j = 0; j < mat_ncols(l3) ; ++j) {
1980  T a(0);
1981  for (size_type k = 0; k < mat_nrows(l2) ; ++k) a += l1(i, k)*l2(k, j);
1982  l3(i, j) = a;
1983  }
1984  }
1985 
1986  // row x col matrix-matrix mult
1987 
1988  template <typename L1, typename L2, typename L3>
1989  void mult_row_col_with_temp(const L1& l1, const L2& l2, L3& l3, col_major) {
1990  typedef typename temporary_col_matrix<L1>::matrix_type temp_col_mat;
1991  temp_col_mat temp(mat_nrows(l1), mat_ncols(l1));
1992  copy(l1, temp);
1993  mult(temp, l2, l3);
1994  }
1995 
1996  template <typename L1, typename L2, typename L3>
1997  void mult_row_col_with_temp(const L1& l1, const L2& l2, L3& l3, row_major) {
1998  typedef typename temporary_row_matrix<L2>::matrix_type temp_row_mat;
1999  temp_row_mat temp(mat_nrows(l2), mat_ncols(l2));
2000  copy(l2, temp);
2001  mult(l1, temp, l3);
2002  }
2003 
2004  template <typename L1, typename L2, typename L3>
2005  void mult_spec(const L1& l1, const L2& l2, L3& l3, rcmult) {
2006  if (is_sparse(l1) && is_sparse(l2)) {
2007  GMM_WARNING3("Inefficient row matrix - col matrix mult for "
2008  "sparse matrices, using temporary");
2009  mult_row_col_with_temp(l1, l2, l3,
2010  typename principal_orientation_type<typename
2011  linalg_traits<L3>::sub_orientation>::potype());
2012  }
2013  else {
2014  auto it2b = linalg_traits<L2>::col_begin(l2), it2 = it2b,
2015  ite = linalg_traits<L2>::col_end(l2);
2016  size_type i,j, k = mat_nrows(l1);
2017 
2018  for (i = 0; i < k; ++i) {
2019  typename linalg_traits<L1>::const_sub_row_type r1=mat_const_row(l1, i);
2020  for (it2 = it2b, j = 0; it2 != ite; ++it2, ++j)
2021  l3(i,j) = vect_sp(r1, linalg_traits<L2>::col(it2));
2022  }
2023  }
2024  }
2025 
2026  // row - row matrix-matrix mult
2027 
2028  template <typename L1, typename L2, typename L3> inline
2029  void mult_spec(const L1& l1, const L2& l2, L3& l3, r_mult) {
2030  mult_spec(l1, l2, l3,r_mult(),typename linalg_traits<L1>::storage_type());
2031  }
2032 
2033  template <typename L1, typename L2, typename L3>
2034  void mult_spec(const L1& l1, const L2& l2, L3& l3, r_mult, abstract_dense) {
2035  // optimizable
2036  clear(l3);
2037  size_type nn = mat_nrows(l3), mm = mat_nrows(l2);
2038  for (size_type i = 0; i < nn; ++i) {
2039  for (size_type j = 0; j < mm; ++j) {
2040  add(scaled(mat_const_row(l2, j), l1(i, j)), mat_row(l3, i));
2041  }
2042  }
2043  }
2044 
2045  template <typename L1, typename L2, typename L3>
2046  void mult_spec(const L1& l1, const L2& l2, L3& l3, r_mult, abstract_sparse) {
2047  // optimizable
2048  clear(l3);
2049  size_type nn = mat_nrows(l3);
2050  for (size_type i = 0; i < nn; ++i) {
2051  typename linalg_traits<L1>::const_sub_row_type rl1=mat_const_row(l1, i);
2052  auto it = vect_const_begin(rl1), ite = vect_const_end(rl1);
2053  for (; it != ite; ++it)
2054  add(scaled(mat_const_row(l2, it.index()), *it), mat_row(l3, i));
2055  }
2056  }
2057 
2058  template <typename L1, typename L2, typename L3>
2059  void mult_spec(const L1& l1, const L2& l2, L3& l3, r_mult, abstract_skyline)
2060  { mult_spec(l1, l2, l3, r_mult(), abstract_sparse()); }
2061 
2062  // col - col matrix-matrix mult
2063 
2064  template <typename L1, typename L2, typename L3> inline
2065  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult) {
2066  mult_spec(l1, l2,l3,c_mult(),typename linalg_traits<L2>::storage_type(),
2067  typename linalg_traits<L2>::sub_orientation());
2068  }
2069 
2070 
2071  template <typename L1, typename L2, typename L3, typename ORIEN>
2072  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult,
2073  abstract_dense, ORIEN) {
2074  typedef typename linalg_traits<L2>::value_type T;
2075  size_type nn = mat_ncols(l3), mm = mat_ncols(l1);
2076 
2077  for (size_type i = 0; i < nn; ++i) {
2078  clear(mat_col(l3, i));
2079  for (size_type j = 0; j < mm; ++j) {
2080  T b = l2(j, i);
2081  if (b != T(0)) add(scaled(mat_const_col(l1, j), b), mat_col(l3, i));
2082  }
2083  }
2084  }
2085 
2086  template <typename L1, typename L2, typename L3, typename ORIEN>
2087  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult,
2088  abstract_sparse, ORIEN) {
2089  // optimizable
2090  clear(l3);
2091  size_type nn = mat_ncols(l3);
2092  for (size_type i = 0; i < nn; ++i) {
2093  typename linalg_traits<L2>::const_sub_col_type rc2 = mat_const_col(l2, i);
2094  auto it = vect_const_begin(rc2), ite = vect_const_end(rc2);
2095  for (; it != ite; ++it)
2096  add(scaled(mat_const_col(l1, it.index()), *it), mat_col(l3, i));
2097  }
2098  }
2099 
2100  template <typename L1, typename L2, typename L3>
2101  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult,
2102  abstract_sparse, row_major) {
2103  typedef typename linalg_traits<L2>::value_type T;
2104  GMM_WARNING3("Inefficient matrix-matrix mult for sparse matrices");
2105  clear(l3);
2106  size_type mm = mat_nrows(l2), nn = mat_ncols(l3);
2107  for (size_type i = 0; i < nn; ++i)
2108  for (size_type j = 0; j < mm; ++j) {
2109  T a = l2(i,j);
2110  if (a != T(0)) add(scaled(mat_const_col(l1, j), a), mat_col(l3, i));
2111  }
2112  }
2113 
2114  template <typename L1, typename L2, typename L3, typename ORIEN> inline
2115  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult,
2116  abstract_skyline, ORIEN)
2117  { mult_spec(l1, l2, l3, c_mult(), abstract_sparse(), ORIEN()); }
2118 
2119 
2120  // col - row matrix-matrix mult
2121 
2122  template <typename L1, typename L2, typename L3> inline
2123  void mult_spec(const L1& l1, const L2& l2, L3& l3, crmult)
2124  { mult_spec(l1,l2,l3,crmult(), typename linalg_traits<L1>::storage_type()); }
2125 
2126 
2127  template <typename L1, typename L2, typename L3>
2128  void mult_spec(const L1& l1, const L2& l2, L3& l3, crmult, abstract_dense) {
2129  // optimizable
2130  clear(l3);
2131  size_type nn = mat_ncols(l1), mm = mat_nrows(l1);
2132  for (size_type i = 0; i < nn; ++i) {
2133  for (size_type j = 0; j < mm; ++j)
2134  add(scaled(mat_const_row(l2, i), l1(j, i)), mat_row(l3, j));
2135  }
2136  }
2137 
2138  template <typename L1, typename L2, typename L3>
2139  void mult_spec(const L1& l1, const L2& l2, L3& l3, crmult, abstract_sparse) {
2140  // optimizable
2141  clear(l3);
2142  size_type nn = mat_ncols(l1);
2143  for (size_type i = 0; i < nn; ++i) {
2144  typename linalg_traits<L1>::const_sub_col_type rc1 = mat_const_col(l1, i);
2145  auto it = vect_const_begin(rc1), ite = vect_const_end(rc1);
2146  for (; it != ite; ++it)
2147  add(scaled(mat_const_row(l2, i), *it), mat_row(l3, it.index()));
2148  }
2149  }
2150 
2151  template <typename L1, typename L2, typename L3> inline
2152  void mult_spec(const L1& l1, const L2& l2, L3& l3, crmult, abstract_skyline)
2153  { mult_spec(l1, l2, l3, crmult(), abstract_sparse()); }
2154 
2155 
2156  /* ******************************************************************** */
2157  /* Symmetry test. */
2158  /* ******************************************************************** */
2159 
2160  ///@endcond
2161  /** test if A is symmetric.
2162  @param A a matrix.
2163  @param tol a threshold.
2164  */
2165  template <typename MAT> inline
2166  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol
2167  = magnitude_of_linalg(MAT)(-1)) {
2168  typedef magnitude_of_linalg(MAT) R;
2169  if (tol < R(0)) tol = default_tol(R()) * mat_maxnorm(A);
2170  if (mat_nrows(A) != mat_ncols(A)) return false;
2171  return is_symmetric(A, tol, typename linalg_traits<MAT>::storage_type());
2172  }
2173  ///@cond DOXY_SHOW_ALL_FUNCTIONS
2174 
2175  template <typename MAT>
2176  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2177  abstract_dense) {
2178  size_type m = mat_nrows(A);
2179  for (size_type i = 1; i < m; ++i)
2180  for (size_type j = 0; j < i; ++j)
2181  if (gmm::abs(A(i, j)-A(j, i)) > tol) return false;
2182  return true;
2183  }
2184 
2185  template <typename MAT>
2186  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2187  abstract_sparse) {
2188  return is_symmetric(A, tol, typename principal_orientation_type<typename
2189  linalg_traits<MAT>::sub_orientation>::potype());
2190  }
2191 
2192  template <typename MAT>
2193  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2194  row_major) {
2195  for (size_type i = 0; i < mat_nrows(A); ++i) {
2196  typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, i);
2197  auto it = vect_const_begin(row), ite = vect_const_end(row);
2198  for (; it != ite; ++it)
2199  if (gmm::abs(*it - A(it.index(), i)) > tol) return false;
2200  }
2201  return true;
2202  }
2203 
2204  template <typename MAT>
2205  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2206  col_major) {
2207  for (size_type i = 0; i < mat_ncols(A); ++i) {
2208  typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, i);
2209  auto it = vect_const_begin(col), ite = vect_const_end(col);
2210  for (; it != ite; ++it)
2211  if (gmm::abs(*it - A(i, it.index())) > tol) return false;
2212  }
2213  return true;
2214  }
2215 
2216  template <typename MAT>
2217  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2218  abstract_skyline)
2219  { return is_symmetric(A, tol, abstract_sparse()); }
2220 
2221  ///@endcond
2222  /** test if A is Hermitian.
2223  @param A a matrix.
2224  @param tol a threshold.
2225  */
2226  template <typename MAT> inline
2227  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol
2228  = magnitude_of_linalg(MAT)(-1)) {
2229  typedef magnitude_of_linalg(MAT) R;
2230  if (tol < R(0)) tol = default_tol(R()) * mat_maxnorm(A);
2231  if (mat_nrows(A) != mat_ncols(A)) return false;
2232  return is_hermitian(A, tol, typename linalg_traits<MAT>::storage_type());
2233  }
2234  ///@cond DOXY_SHOW_ALL_FUNCTIONS
2235 
2236  template <typename MAT>
2237  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2238  abstract_dense) {
2239  size_type m = mat_nrows(A);
2240  for (size_type i = 1; i < m; ++i)
2241  for (size_type j = 0; j < i; ++j)
2242  if (gmm::abs(A(i, j)-gmm::conj(A(j, i))) > tol) return false;
2243  return true;
2244  }
2245 
2246  template <typename MAT>
2247  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2248  abstract_sparse) {
2249  return is_hermitian(A, tol, typename principal_orientation_type<typename
2250  linalg_traits<MAT>::sub_orientation>::potype());
2251  }
2252 
2253  template <typename MAT>
2254  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2255  row_major) {
2256  for (size_type i = 0; i < mat_nrows(A); ++i) {
2257  typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, i);
2258  auto it = vect_const_begin(row), ite = vect_const_end(row);
2259  for (; it != ite; ++it)
2260  if (gmm::abs(gmm::conj(*it) - A(it.index(), i)) > tol) return false;
2261  }
2262  return true;
2263  }
2264 
2265  template <typename MAT>
2266  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2267  col_major) {
2268  for (size_type i = 0; i < mat_ncols(A); ++i) {
2269  typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, i);
2270  auto it = vect_const_begin(col), ite = vect_const_end(col);
2271  for (; it != ite; ++it)
2272  if (gmm::abs(gmm::conj(*it) - A(i, it.index())) > tol) return false;
2273  }
2274  return true;
2275  }
2276 
2277  template <typename MAT>
2278  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2279  abstract_skyline)
2280  { return is_hermitian(A, tol, abstract_sparse()); }
2281  ///@endcond
2282 }
2283 
2284 
2285 #endif // GMM_BLAS_H__
gmm::fill_random
void fill_random(L &l)
fill a vector or matrix with random value (uniform [-1,1]).
Definition: gmm_blas.h:129
gmm::mat_norminf
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_norminf(const M &m)
*/
Definition: gmm_blas.h:836
gmm::mat_maxnorm
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_maxnorm(const M &m)
*/
Definition: gmm_blas.h:870
gmm::mat_euclidean_norm
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_euclidean_norm(const M &m)
Euclidean norm of a matrix.
Definition: gmm_blas.h:636
gmm::vect_norm2_sqr
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2_sqr(const V &v)
squared Euclidean norm of a vector.
Definition: gmm_blas.h:544
gmm::mat_norm1
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_norm1(const M &m)
*/
Definition: gmm_blas.h:782
gmm_scaled.h
get a scaled view of a vector/matrix.
gmm_conjugated.h
handle conjugation of complex matrices/vectors.
gmm::is_hermitian
bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol=magnitude_of_linalg(MAT)(-1))
*/
Definition: gmm_blas.h:2227
gmm::vect_dist2
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist2(const V1 &v1, const V2 &v2)
Euclidean distance between two vectors.
Definition: gmm_blas.h:597
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
gmm::vect_norm1
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm1(const V &v)
1-norm of a vector
Definition: gmm_blas.h:646
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::conjugated
conjugated_return< L >::return_type conjugated(const L &v)
return a conjugated view of the input matrix or vector.
Definition: gmm_conjugated.h:294
gmm::copy
void copy(const L1 &l1, L2 &l2)
*/
Definition: gmm_blas.h:977
gmm::vect_norminf
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norminf(const V &v)
Infinity norm of a vector.
Definition: gmm_blas.h:693
gmm::is_symmetric
bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol=magnitude_of_linalg(MAT)(-1))
*/
Definition: gmm_blas.h:2166
gmm::resize
void resize(V &v, size_type n)
*/
Definition: gmm_blas.h:209
gmm::vect_dist1
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist1(const V1 &v1, const V2 &v2)
1-distance between two vectors
Definition: gmm_blas.h:658
gmm::fill
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*/
Definition: gmm_blas.h:103
gmm::vect_dist2_sqr
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist2_sqr(const V1 &v1, const V2 &v2)
squared Euclidean distance between two vectors
Definition: gmm_blas.h:565
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::clean
void clean(L &l, double threshold)
Clean a vector or matrix (replace near-zero entries with zeroes).
gmm::vect_distinf
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_distinf(const V1 &v1, const V2 &v2)
Infinity distance between two vectors.
Definition: gmm_blas.h:705
gmm::mat_trace
linalg_traits< M >::value_type mat_trace(const M &m)
Trace of a matrix.
Definition: gmm_blas.h:528
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
gmm::mat_euclidean_norm_sqr
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_euclidean_norm_sqr(const M &m)
*/
Definition: gmm_blas.h:626