GetFEM  5.4.2
gmm_blas.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-2020 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file gmm_blas.h
33  @author Yves Renard <[email protected]>
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), "dimensions mismatch, "
1441  << vect_size(l1) << " !=" << vect_size(l2));
1442  GMM_ASSERT2(vect_size(l1) == vect_size(l3), "dimensions mismatch, "
1443  << vect_size(l1) << " !=" << vect_size(l3));
1444  if ((const void *)(&l1) == (const void *)(&l3))
1445  add(l2, l3);
1446  else if ((const void *)(&l2) == (const void *)(&l3))
1447  add(l1, l3);
1448  else
1449  add(l1, l2, l3, typename linalg_traits<L1>::storage_type(),
1450  typename linalg_traits<L2>::storage_type(),
1451  typename linalg_traits<L3>::storage_type());
1452  }
1453 
1454  template <typename IT1, typename IT2, typename IT3>
1455  void add_full_(IT1 it1, IT2 it2, IT3 it3, IT3 ite) {
1456  for (; it3 != ite; ++it3, ++it2, ++it1) *it3 = *it1 + *it2;
1457  }
1458 
1459  template <typename IT1, typename IT2, typename IT3>
1460  void add_almost_full_(IT1 it1, IT1 ite1, IT2 it2, IT3 it3, IT3 ite3) {
1461  IT3 it = it3;
1462  for (; it != ite3; ++it, ++it2) *it = *it2;
1463  for (; it1 != ite1; ++it1)
1464  *(it3 + it1.index()) += *it1;
1465  }
1466 
1467  template <typename IT1, typename IT2, typename IT3>
1468  void add_to_full_(IT1 it1, IT1 ite1, IT2 it2, IT2 ite2,
1469  IT3 it3, IT3 ite3) {
1470  typedef typename std::iterator_traits<IT3>::value_type T;
1471  IT3 it = it3;
1472  for (; it != ite3; ++it) *it = T(0);
1473  for (; it1 != ite1; ++it1) *(it3 + it1.index()) = *it1;
1474  for (; it2 != ite2; ++it2) *(it3 + it2.index()) += *it2;
1475  }
1476 
1477  template <typename L1, typename L2, typename L3> inline
1478  void add(const L1& l1, const L2& l2, L3& l3,
1479  abstract_dense, abstract_dense, abstract_dense) {
1480  add_full_(vect_const_begin(l1), vect_const_begin(l2),
1481  vect_begin(l3), vect_end(l3));
1482  }
1483 
1484  // generic function for add(v1, v2, v3).
1485  // Need to be specialized to optimize particular additions.
1486  template <typename L1, typename L2, typename L3,
1487  typename ST1, typename ST2, typename ST3>
1488  inline void add(const L1& l1, const L2& l2, L3& l3, ST1, ST2, ST3)
1489  { copy(l2, l3); add(l1, l3, ST1(), ST3()); }
1490 
1491  template <typename L1, typename L2, typename L3> inline
1492  void add(const L1& l1, const L2& l2, L3& l3,
1493  abstract_sparse, abstract_dense, abstract_dense) {
1494  add_almost_full_(vect_const_begin(l1), vect_const_end(l1),
1495  vect_const_begin(l2), vect_begin(l3), vect_end(l3));
1496  }
1497 
1498  template <typename L1, typename L2, typename L3> inline
1499  void add(const L1& l1, const L2& l2, L3& l3,
1500  abstract_dense, abstract_sparse, abstract_dense)
1501  { add(l2, l1, l3, abstract_sparse(), abstract_dense(), abstract_dense()); }
1502 
1503  template <typename L1, typename L2, typename L3> inline
1504  void add(const L1& l1, const L2& l2, L3& l3,
1505  abstract_sparse, abstract_sparse, abstract_dense) {
1506  add_to_full_(vect_const_begin(l1), vect_const_end(l1),
1507  vect_const_begin(l2), vect_const_end(l2),
1508  vect_begin(l3), vect_end(l3));
1509  }
1510 
1511 
1512  template <typename L1, typename L2, typename L3>
1513  void add_spspsp(const L1& l1, const L2& l2, L3& l3, linalg_true) {
1514  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1515  auto it2 = vect_const_begin(l2), ite2 = vect_const_end(l2);
1516  clear(l3);
1517  while (it1 != ite1 && it2 != ite2) {
1518  ptrdiff_t d = it1.index() - it2.index();
1519  if (d < 0)
1520  { l3[it1.index()] += *it1; ++it1; }
1521  else if (d > 0)
1522  { l3[it2.index()] += *it2; ++it2; }
1523  else
1524  { l3[it1.index()] = *it1 + *it2; ++it1; ++it2; }
1525  }
1526  for (; it1 != ite1; ++it1) l3[it1.index()] += *it1;
1527  for (; it2 != ite2; ++it2) l3[it2.index()] += *it2;
1528  }
1529 
1530  template <typename L1, typename L2, typename L3>
1531  void add_spspsp(const L1& l1, const L2& l2, L3& l3, linalg_false)
1532  { copy(l2, l3); add(l2, l3); }
1533 
1534  template <typename L1, typename L2, typename L3>
1535  void add(const L1& l1, const L2& l2, L3& l3,
1536  abstract_sparse, abstract_sparse, abstract_sparse) {
1537  add_spspsp(l1, l2, l3, typename linalg_and<typename
1538  linalg_traits<L1>::index_sorted,
1539  typename linalg_traits<L2>::index_sorted>::bool_type());
1540  }
1541 
1542  template <typename L1, typename L2>
1543  void add(const L1& l1, L2& l2, abstract_dense, abstract_dense) {
1544  auto it1 = vect_const_begin(l1);
1545  auto it2 = vect_begin(l2), ite = vect_end(l2);
1546  for (; it2 != ite; ++it2, ++it1) *it2 += *it1;
1547  }
1548 
1549  template <typename L1, typename L2>
1550  void add(const L1& l1, L2& l2, abstract_dense, abstract_skyline) {
1551  typedef typename linalg_traits<L1>::value_type T;
1552 
1553  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1554  size_type i1 = 0, ie1 = vect_size(l1);
1555  while (it1 != ite1 && *it1 == T(0)) { ++it1; ++i1; }
1556  if (it1 != ite1) {
1557  auto it2 = vect_begin(l2), ite2 = vect_end(l2);
1558  while (ie1 && *(ite1-1) == T(0)) { ite1--; --ie1; }
1559 
1560  if (it2 == ite2 || i1 < it2.index()) {
1561  l2[i1] = *it1; ++i1; ++it1;
1562  if (it1 == ite1) return;
1563  it2 = vect_begin(l2); ite2 = vect_end(l2);
1564  }
1565  if (ie1 > ite2.index()) {
1566  --ite1; l2[ie1 - 1] = *ite1;
1567  it2 = vect_begin(l2);
1568  }
1569  it2 += i1 - it2.index();
1570  for (; it1 != ite1; ++it1, ++it2) { *it2 += *it1; }
1571  }
1572  }
1573 
1574 
1575  template <typename L1, typename L2>
1576  void add(const L1& l1, L2& l2, abstract_skyline, abstract_dense) {
1577  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1578  if (it1 != ite1) {
1579  auto it2 = vect_begin(l2);
1580  it2 += it1.index();
1581  for (; it1 != ite1; ++it2, ++it1) *it2 += *it1;
1582  }
1583  }
1584 
1585 
1586  template <typename L1, typename L2>
1587  void add(const L1& l1, L2& l2, abstract_sparse, abstract_dense) {
1588  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1589  for (; it1 != ite1; ++it1) l2[it1.index()] += *it1;
1590  }
1591 
1592  template <typename L1, typename L2>
1593  void add(const L1& l1, L2& l2, abstract_sparse, abstract_sparse) {
1594  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1595  for (; it1 != ite1; ++it1) l2[it1.index()] += *it1;
1596  }
1597 
1598  template <typename L1, typename L2>
1599  void add(const L1& l1, L2& l2, abstract_sparse, abstract_skyline) {
1600  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1601  for (; it1 != ite1; ++it1) l2[it1.index()] += *it1;
1602  }
1603 
1604 
1605  template <typename L1, typename L2>
1606  void add(const L1& l1, L2& l2, abstract_skyline, abstract_sparse) {
1607  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1608  for (; it1 != ite1; ++it1)
1609  if (*it1 != typename linalg_traits<L1>::value_type(0))
1610  l2[it1.index()] += *it1;
1611  }
1612 
1613  template <typename L1, typename L2>
1614  void add(const L1& l1, L2& l2, abstract_skyline, abstract_skyline) {
1615  typedef typename linalg_traits<L1>::value_type T1;
1616  typedef typename linalg_traits<L2>::value_type T2;
1617 
1618  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1619 
1620  while (it1 != ite1 && *it1 == T1(0)) ++it1;
1621  if (ite1 != it1) {
1622  auto it2 = vect_begin(l2), ite2 = vect_end(l2);
1623  while (*(ite1-1) == T1(0)) ite1--;
1624  if (it2 == ite2 || it1.index() < it2.index()) {
1625  l2[it1.index()] = T2(0);
1626  it2 = vect_begin(l2); ite2 = vect_end(l2);
1627  }
1628  if (ite1.index() > ite2.index()) {
1629  l2[ite1.index() - 1] = T2(0);
1630  it2 = vect_begin(l2);
1631  }
1632  it2 += it1.index() - it2.index();
1633  for (; it1 != ite1; ++it1, ++it2) *it2 += *it1;
1634  }
1635  }
1636 
1637  template <typename L1, typename L2>
1638  void add(const L1& l1, L2& l2, abstract_dense, abstract_sparse) {
1639  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1640  for (size_type i = 0; it1 != ite1; ++it1, ++i)
1641  if (*it1 != typename linalg_traits<L1>::value_type(0)) l2[i] += *it1;
1642  }
1643 
1644  /* ******************************************************************** */
1645  /* Matrix-vector mult */
1646  /* ******************************************************************** */
1647  ///@endcond
1648  /** matrix-vector or matrix-matrix product.
1649  @param l1 a matrix.
1650  @param l2 a vector or matrix.
1651  @param l3 the product l1*l2.
1652  */
1653  template <typename L1, typename L2, typename L3> inline
1654  void mult(const L1& l1, const L2& l2, L3& l3) {
1655  mult_dispatch(l1, l2, l3, typename linalg_traits<L2>::linalg_type());
1656  }
1657  ///@cond DOXY_SHOW_ALL_FUNCTIONS
1658 
1659  template <typename L1, typename L2, typename L3> inline
1660  void mult(const L1& l1, const L2& l2, const L3& l3)
1661  { mult(l1, l2, linalg_const_cast(l3)); }
1662 
1663  template <typename L1, typename L2, typename L3> inline
1664  void mult_dispatch(const L1& l1, const L2& l2, L3& l3, abstract_vector) {
1665  size_type m = mat_nrows(l1), n = mat_ncols(l1);
1666  if (!m || !n) { gmm::clear(l3); return; }
1667  GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l3), "dimensions mismatch");
1668  if (!same_origin(l2, l3))
1669  mult_spec(l1, l2, l3, typename principal_orientation_type<typename
1670  linalg_traits<L1>::sub_orientation>::potype());
1671  else {
1672  GMM_WARNING2("Warning, A temporary is used for mult\n");
1673  typename temporary_vector<L3>::vector_type temp(vect_size(l3));
1674  mult_spec(l1, l2, temp, typename principal_orientation_type<typename
1675  linalg_traits<L1>::sub_orientation>::potype());
1676  copy(temp, l3);
1677  }
1678  }
1679 
1680  template <typename L1, typename L2, typename L3>
1681  void mult_by_row(const L1& l1, const L2& l2, L3& l3, abstract_sparse) {
1682  typedef typename linalg_traits<L3>::value_type T;
1683  clear(l3);
1684  size_type nr = mat_nrows(l1);
1685  for (size_type i = 0; i < nr; ++i) {
1686  T aux = vect_sp(mat_const_row(l1, i), l2);
1687  if (aux != T(0)) l3[i] = aux;
1688  }
1689  }
1690 
1691  template <typename L1, typename L2, typename L3>
1692  void mult_by_row(const L1& l1, const L2& l2, L3& l3, abstract_skyline) {
1693  typedef typename linalg_traits<L3>::value_type T;
1694  clear(l3);
1695  size_type nr = mat_nrows(l1);
1696  for (size_type i = 0; i < nr; ++i) {
1697  T aux = vect_sp(mat_const_row(l1, i), l2);
1698  if (aux != T(0)) l3[i] = aux;
1699  }
1700  }
1701 
1702  template <typename L1, typename L2, typename L3>
1703  void mult_by_row(const L1& l1, const L2& l2, L3& l3, abstract_dense) {
1704  typename linalg_traits<L3>::iterator it=vect_begin(l3), ite=vect_end(l3);
1705  auto itr = mat_row_const_begin(l1);
1706  for (; it != ite; ++it, ++itr)
1707  *it = vect_sp(linalg_traits<L1>::row(itr), l2,
1708  typename linalg_traits<L1>::storage_type(),
1709  typename linalg_traits<L2>::storage_type());
1710  }
1711 
1712  template <typename L1, typename L2, typename L3>
1713  void mult_by_col(const L1& l1, const L2& l2, L3& l3, abstract_dense) {
1714  clear(l3);
1715  size_type nc = mat_ncols(l1);
1716  for (size_type i = 0; i < nc; ++i)
1717  add(scaled(mat_const_col(l1, i), l2[i]), l3);
1718  }
1719 
1720  template <typename L1, typename L2, typename L3>
1721  void mult_by_col(const L1& l1, const L2& l2, L3& l3, abstract_sparse) {
1722  typedef typename linalg_traits<L2>::value_type T;
1723  clear(l3);
1724  auto it = vect_const_begin(l2), ite = vect_const_end(l2);
1725  for (; it != ite; ++it)
1726  if (*it != T(0)) add(scaled(mat_const_col(l1, it.index()), *it), l3);
1727  }
1728 
1729  template <typename L1, typename L2, typename L3>
1730  void mult_by_col(const L1& l1, const L2& l2, L3& l3, abstract_skyline) {
1731  typedef typename linalg_traits<L2>::value_type T;
1732  clear(l3);
1733  auto it = vect_const_begin(l2), ite = vect_const_end(l2);
1734  for (; it != ite; ++it)
1735  if (*it != T(0)) add(scaled(mat_const_col(l1, it.index()), *it), l3);
1736  }
1737 
1738  template <typename L1, typename L2, typename L3> inline
1739  void mult_spec(const L1& l1, const L2& l2, L3& l3, row_major)
1740  { mult_by_row(l1, l2, l3, typename linalg_traits<L3>::storage_type()); }
1741 
1742  template <typename L1, typename L2, typename L3> inline
1743  void mult_spec(const L1& l1, const L2& l2, L3& l3, col_major)
1744  { mult_by_col(l1, l2, l3, typename linalg_traits<L2>::storage_type()); }
1745 
1746  template <typename L1, typename L2, typename L3> inline
1747  void mult_spec(const L1& l1, const L2& l2, L3& l3, abstract_null_type)
1748  { mult_ind(l1, l2, l3, typename linalg_traits<L1>::storage_type()); }
1749 
1750  template <typename L1, typename L2, typename L3>
1751  void mult_ind(const L1& l1, const L2& l2, L3& l3, abstract_indirect) {
1752  GMM_ASSERT1(false, "gmm::mult(m, ., .) undefined for this kind of matrix");
1753  }
1754 
1755  template <typename L1, typename L2, typename L3, typename L4> inline
1756  void mult(const L1& l1, const L2& l2, const L3& l3, L4& l4) {
1757  size_type m = mat_nrows(l1), n = mat_ncols(l1);
1758  copy(l3, l4);
1759  if (!m || !n) { gmm::copy(l3, l4); return; }
1760  GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l4), "dimensions mismatch");
1761  if (!same_origin(l2, l4)) {
1762  mult_add_spec(l1, l2, l4, typename principal_orientation_type<typename
1763  linalg_traits<L1>::sub_orientation>::potype());
1764  }
1765  else {
1766  GMM_WARNING2("Warning, A temporary is used for mult\n");
1767  typename temporary_vector<L2>::vector_type temp(vect_size(l2));
1768  copy(l2, temp);
1769  mult_add_spec(l1,temp, l4, typename principal_orientation_type<typename
1770  linalg_traits<L1>::sub_orientation>::potype());
1771  }
1772  }
1773 
1774  template <typename L1, typename L2, typename L3, typename L4> inline
1775  void mult(const L1& l1, const L2& l2, const L3& l3, const L4& l4)
1776  { mult(l1, l2, l3, linalg_const_cast(l4)); }
1777 
1778  ///@endcond
1779  /** Multiply-accumulate. l3 += l1*l2; */
1780  template <typename L1, typename L2, typename L3> inline
1781  void mult_add(const L1& l1, const L2& l2, L3& l3) {
1782  size_type m = mat_nrows(l1), n = mat_ncols(l1);
1783  if (!m || !n) return;
1784  GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l3), "dimensions mismatch");
1785  if (!same_origin(l2, l3)) {
1786  mult_add_spec(l1, l2, l3, typename principal_orientation_type<typename
1787  linalg_traits<L1>::sub_orientation>::potype());
1788  }
1789  else {
1790  GMM_WARNING2("Warning, A temporary is used for mult\n");
1791  typename temporary_vector<L3>::vector_type temp(vect_size(l2));
1792  copy(l2, temp);
1793  mult_add_spec(l1,temp, l3, typename principal_orientation_type<typename
1794  linalg_traits<L1>::sub_orientation>::potype());
1795  }
1796  }
1797  ///@cond DOXY_SHOW_ALL_FUNCTIONS
1798 
1799  template <typename L1, typename L2, typename L3> inline
1800  void mult_add(const L1& l1, const L2& l2, const L3& l3)
1801  { mult_add(l1, l2, linalg_const_cast(l3)); }
1802 
1803  template <typename L1, typename L2, typename L3>
1804  void mult_add_by_row(const L1& l1, const L2& l2, L3& l3, abstract_sparse) {
1805  typedef typename linalg_traits<L3>::value_type T;
1806  size_type nr = mat_nrows(l1);
1807  for (size_type i = 0; i < nr; ++i) {
1808  T aux = vect_sp(mat_const_row(l1, i), l2);
1809  if (aux != T(0)) l3[i] += aux;
1810  }
1811  }
1812 
1813  template <typename L1, typename L2, typename L3>
1814  void mult_add_by_row(const L1& l1, const L2& l2, L3& l3, abstract_skyline) {
1815  typedef typename linalg_traits<L3>::value_type T;
1816  size_type nr = mat_nrows(l1);
1817  for (size_type i = 0; i < nr; ++i) {
1818  T aux = vect_sp(mat_const_row(l1, i), l2);
1819  if (aux != T(0)) l3[i] += aux;
1820  }
1821  }
1822 
1823  template <typename L1, typename L2, typename L3>
1824  void mult_add_by_row(const L1& l1, const L2& l2, L3& l3, abstract_dense) {
1825  auto it=vect_begin(l3), ite=vect_end(l3);
1826  auto itr = mat_row_const_begin(l1);
1827  for (; it != ite; ++it, ++itr)
1828  *it += vect_sp(linalg_traits<L1>::row(itr), l2);
1829  }
1830 
1831  template <typename L1, typename L2, typename L3>
1832  void mult_add_by_col(const L1& l1, const L2& l2, L3& l3, abstract_dense) {
1833  size_type nc = mat_ncols(l1);
1834  for (size_type i = 0; i < nc; ++i)
1835  add(scaled(mat_const_col(l1, i), l2[i]), l3);
1836  }
1837 
1838  template <typename L1, typename L2, typename L3>
1839  void mult_add_by_col(const L1& l1, const L2& l2, L3& l3, abstract_sparse) {
1840  auto it = vect_const_begin(l2), ite = vect_const_end(l2);
1841  for (; it != ite; ++it)
1842  if (*it != typename linalg_traits<L2>::value_type(0))
1843  add(scaled(mat_const_col(l1, it.index()), *it), l3);
1844  }
1845 
1846  template <typename L1, typename L2, typename L3>
1847  void mult_add_by_col(const L1& l1, const L2& l2, L3& l3, abstract_skyline) {
1848  auto it = vect_const_begin(l2), ite = vect_const_end(l2);
1849  for (; it != ite; ++it)
1850  if (*it != typename linalg_traits<L2>::value_type(0))
1851  add(scaled(mat_const_col(l1, it.index()), *it), l3);
1852  }
1853 
1854  template <typename L1, typename L2, typename L3> inline
1855  void mult_add_spec(const L1& l1, const L2& l2, L3& l3, row_major)
1856  { mult_add_by_row(l1, l2, l3, typename linalg_traits<L3>::storage_type()); }
1857 
1858  template <typename L1, typename L2, typename L3> inline
1859  void mult_add_spec(const L1& l1, const L2& l2, L3& l3, col_major)
1860  { mult_add_by_col(l1, l2, l3, typename linalg_traits<L2>::storage_type()); }
1861 
1862  template <typename L1, typename L2, typename L3> inline
1863  void mult_add_spec(const L1& l1, const L2& l2, L3& l3, abstract_null_type)
1864  { mult_ind(l1, l2, l3, typename linalg_traits<L1>::storage_type()); }
1865 
1866  template <typename L1, typename L2, typename L3>
1867  void transposed_mult(const L1& l1, const L2& l2, const L3& l3)
1868  { mult(gmm::transposed(l1), l2, l3); }
1869 
1870 
1871  /* ******************************************************************** */
1872  /* Matrix-matrix mult */
1873  /* ******************************************************************** */
1874 
1875 
1876  struct g_mult {}; // generic mult, less optimized
1877  struct c_mult {}; // col x col -> col mult
1878  struct r_mult {}; // row x row -> row mult
1879  struct rcmult {}; // row x col -> col mult
1880  struct crmult {}; // col x row -> row mult
1881 
1882 
1883  template<typename SO1, typename SO2, typename SO3> struct mult_t;
1884  #define DEFMU__ template<> struct mult_t
1885  DEFMU__<row_major , row_major , row_major > { typedef r_mult t; };
1886  DEFMU__<row_major , row_major , col_major > { typedef g_mult t; };
1887  DEFMU__<row_major , row_major , col_and_row> { typedef r_mult t; };
1888  DEFMU__<row_major , row_major , row_and_col> { typedef r_mult t; };
1889  DEFMU__<row_major , col_major , row_major > { typedef rcmult t; };
1890  DEFMU__<row_major , col_major , col_major > { typedef rcmult t; };
1891  DEFMU__<row_major , col_major , col_and_row> { typedef rcmult t; };
1892  DEFMU__<row_major , col_major , row_and_col> { typedef rcmult t; };
1893  DEFMU__<row_major , col_and_row, row_major > { typedef r_mult t; };
1894  DEFMU__<row_major , col_and_row, col_major > { typedef rcmult t; };
1895  DEFMU__<row_major , col_and_row, col_and_row> { typedef rcmult t; };
1896  DEFMU__<row_major , col_and_row, row_and_col> { typedef rcmult t; };
1897  DEFMU__<row_major , row_and_col, row_major > { typedef r_mult t; };
1898  DEFMU__<row_major , row_and_col, col_major > { typedef rcmult t; };
1899  DEFMU__<row_major , row_and_col, col_and_row> { typedef r_mult t; };
1900  DEFMU__<row_major , row_and_col, row_and_col> { typedef r_mult t; };
1901  DEFMU__<col_major , row_major , row_major > { typedef crmult t; };
1902  DEFMU__<col_major , row_major , col_major > { typedef g_mult t; };
1903  DEFMU__<col_major , row_major , col_and_row> { typedef crmult t; };
1904  DEFMU__<col_major , row_major , row_and_col> { typedef crmult t; };
1905  DEFMU__<col_major , col_major , row_major > { typedef g_mult t; };
1906  DEFMU__<col_major , col_major , col_major > { typedef c_mult t; };
1907  DEFMU__<col_major , col_major , col_and_row> { typedef c_mult t; };
1908  DEFMU__<col_major , col_major , row_and_col> { typedef c_mult t; };
1909  DEFMU__<col_major , col_and_row, row_major > { typedef crmult t; };
1910  DEFMU__<col_major , col_and_row, col_major > { typedef c_mult t; };
1911  DEFMU__<col_major , col_and_row, col_and_row> { typedef c_mult t; };
1912  DEFMU__<col_major , col_and_row, row_and_col> { typedef c_mult t; };
1913  DEFMU__<col_major , row_and_col, row_major > { typedef crmult t; };
1914  DEFMU__<col_major , row_and_col, col_major > { typedef c_mult t; };
1915  DEFMU__<col_major , row_and_col, col_and_row> { typedef c_mult t; };
1916  DEFMU__<col_major , row_and_col, row_and_col> { typedef c_mult t; };
1917  DEFMU__<col_and_row, row_major , row_major > { typedef r_mult t; };
1918  DEFMU__<col_and_row, row_major , col_major > { typedef c_mult t; };
1919  DEFMU__<col_and_row, row_major , col_and_row> { typedef r_mult t; };
1920  DEFMU__<col_and_row, row_major , row_and_col> { typedef r_mult t; };
1921  DEFMU__<col_and_row, col_major , row_major > { typedef rcmult t; };
1922  DEFMU__<col_and_row, col_major , col_major > { typedef c_mult t; };
1923  DEFMU__<col_and_row, col_major , col_and_row> { typedef c_mult t; };
1924  DEFMU__<col_and_row, col_major , row_and_col> { typedef c_mult t; };
1925  DEFMU__<col_and_row, col_and_row, row_major > { typedef r_mult t; };
1926  DEFMU__<col_and_row, col_and_row, col_major > { typedef c_mult t; };
1927  DEFMU__<col_and_row, col_and_row, col_and_row> { typedef c_mult t; };
1928  DEFMU__<col_and_row, col_and_row, row_and_col> { typedef c_mult t; };
1929  DEFMU__<col_and_row, row_and_col, row_major > { typedef r_mult t; };
1930  DEFMU__<col_and_row, row_and_col, col_major > { typedef c_mult t; };
1931  DEFMU__<col_and_row, row_and_col, col_and_row> { typedef c_mult t; };
1932  DEFMU__<col_and_row, row_and_col, row_and_col> { typedef r_mult t; };
1933  DEFMU__<row_and_col, row_major , row_major > { typedef r_mult t; };
1934  DEFMU__<row_and_col, row_major , col_major > { typedef c_mult t; };
1935  DEFMU__<row_and_col, row_major , col_and_row> { typedef r_mult t; };
1936  DEFMU__<row_and_col, row_major , row_and_col> { typedef r_mult t; };
1937  DEFMU__<row_and_col, col_major , row_major > { typedef rcmult t; };
1938  DEFMU__<row_and_col, col_major , col_major > { typedef c_mult t; };
1939  DEFMU__<row_and_col, col_major , col_and_row> { typedef c_mult t; };
1940  DEFMU__<row_and_col, col_major , row_and_col> { typedef c_mult t; };
1941  DEFMU__<row_and_col, col_and_row, row_major > { typedef rcmult t; };
1942  DEFMU__<row_and_col, col_and_row, col_major > { typedef rcmult t; };
1943  DEFMU__<row_and_col, col_and_row, col_and_row> { typedef rcmult t; };
1944  DEFMU__<row_and_col, col_and_row, row_and_col> { typedef rcmult t; };
1945  DEFMU__<row_and_col, row_and_col, row_major > { typedef r_mult t; };
1946  DEFMU__<row_and_col, row_and_col, col_major > { typedef c_mult t; };
1947  DEFMU__<row_and_col, row_and_col, col_and_row> { typedef r_mult t; };
1948  DEFMU__<row_and_col, row_and_col, row_and_col> { typedef r_mult t; };
1949 
1950  template <typename L1, typename L2, typename L3>
1951  void mult_dispatch(const L1& l1, const L2& l2, L3& l3, abstract_matrix) {
1952  typedef typename temporary_matrix<L3>::matrix_type temp_mat_type;
1953  size_type n = mat_ncols(l1);
1954  if (n == 0) { gmm::clear(l3); return; }
1955  GMM_ASSERT2(n == mat_nrows(l2) && mat_nrows(l1) == mat_nrows(l3) &&
1956  mat_ncols(l2) == mat_ncols(l3), "dimensions mismatch");
1957 
1958  if (same_origin(l2, l3) || same_origin(l1, l3)) {
1959  GMM_WARNING2("A temporary is used for mult");
1960  temp_mat_type temp(mat_nrows(l3), mat_ncols(l3));
1961  mult_spec(l1, l2, temp, typename mult_t<
1962  typename linalg_traits<L1>::sub_orientation,
1963  typename linalg_traits<L2>::sub_orientation,
1964  typename linalg_traits<temp_mat_type>::sub_orientation>::t());
1965  copy(temp, l3);
1966  }
1967  else
1968  mult_spec(l1, l2, l3, typename mult_t<
1969  typename linalg_traits<L1>::sub_orientation,
1970  typename linalg_traits<L2>::sub_orientation,
1971  typename linalg_traits<L3>::sub_orientation>::t());
1972  }
1973 
1974  // Completely generic but inefficient
1975 
1976  template <typename L1, typename L2, typename L3>
1977  void mult_spec(const L1& l1, const L2& l2, L3& l3, g_mult) {
1978  typedef typename linalg_traits<L3>::value_type T;
1979  GMM_WARNING2("Inefficient generic matrix-matrix mult is used");
1980  for (size_type i = 0; i < mat_nrows(l3) ; ++i)
1981  for (size_type j = 0; j < mat_ncols(l3) ; ++j) {
1982  T a(0);
1983  for (size_type k = 0; k < mat_nrows(l2) ; ++k) a += l1(i, k)*l2(k, j);
1984  l3(i, j) = a;
1985  }
1986  }
1987 
1988  // row x col matrix-matrix mult
1989 
1990  template <typename L1, typename L2, typename L3>
1991  void mult_row_col_with_temp(const L1& l1, const L2& l2, L3& l3, col_major) {
1992  typedef typename temporary_col_matrix<L1>::matrix_type temp_col_mat;
1993  temp_col_mat temp(mat_nrows(l1), mat_ncols(l1));
1994  copy(l1, temp);
1995  mult(temp, l2, l3);
1996  }
1997 
1998  template <typename L1, typename L2, typename L3>
1999  void mult_row_col_with_temp(const L1& l1, const L2& l2, L3& l3, row_major) {
2000  typedef typename temporary_row_matrix<L2>::matrix_type temp_row_mat;
2001  temp_row_mat temp(mat_nrows(l2), mat_ncols(l2));
2002  copy(l2, temp);
2003  mult(l1, temp, l3);
2004  }
2005 
2006  template <typename L1, typename L2, typename L3>
2007  void mult_spec(const L1& l1, const L2& l2, L3& l3, rcmult) {
2008  if (is_sparse(l1) && is_sparse(l2)) {
2009  GMM_WARNING3("Inefficient row matrix - col matrix mult for "
2010  "sparse matrices, using temporary");
2011  mult_row_col_with_temp(l1, l2, l3,
2012  typename principal_orientation_type<typename
2013  linalg_traits<L3>::sub_orientation>::potype());
2014  }
2015  else {
2016  auto it2b = linalg_traits<L2>::col_begin(l2), it2 = it2b,
2017  ite = linalg_traits<L2>::col_end(l2);
2018  size_type i,j, k = mat_nrows(l1);
2019 
2020  for (i = 0; i < k; ++i) {
2021  typename linalg_traits<L1>::const_sub_row_type r1=mat_const_row(l1, i);
2022  for (it2 = it2b, j = 0; it2 != ite; ++it2, ++j)
2023  l3(i,j) = vect_sp(r1, linalg_traits<L2>::col(it2));
2024  }
2025  }
2026  }
2027 
2028  // row - row matrix-matrix mult
2029 
2030  template <typename L1, typename L2, typename L3> inline
2031  void mult_spec(const L1& l1, const L2& l2, L3& l3, r_mult) {
2032  mult_spec(l1, l2, l3,r_mult(),typename linalg_traits<L1>::storage_type());
2033  }
2034 
2035  template <typename L1, typename L2, typename L3>
2036  void mult_spec(const L1& l1, const L2& l2, L3& l3, r_mult, abstract_dense) {
2037  // optimizable
2038  clear(l3);
2039  size_type nn = mat_nrows(l3), mm = mat_nrows(l2);
2040  for (size_type i = 0; i < nn; ++i) {
2041  for (size_type j = 0; j < mm; ++j) {
2042  add(scaled(mat_const_row(l2, j), l1(i, j)), mat_row(l3, i));
2043  }
2044  }
2045  }
2046 
2047  template <typename L1, typename L2, typename L3>
2048  void mult_spec(const L1& l1, const L2& l2, L3& l3, r_mult, abstract_sparse) {
2049  // optimizable
2050  clear(l3);
2051  size_type nn = mat_nrows(l3);
2052  for (size_type i = 0; i < nn; ++i) {
2053  typename linalg_traits<L1>::const_sub_row_type rl1=mat_const_row(l1, i);
2054  auto it = vect_const_begin(rl1), ite = vect_const_end(rl1);
2055  for (; it != ite; ++it)
2056  add(scaled(mat_const_row(l2, it.index()), *it), mat_row(l3, i));
2057  }
2058  }
2059 
2060  template <typename L1, typename L2, typename L3>
2061  void mult_spec(const L1& l1, const L2& l2, L3& l3, r_mult, abstract_skyline)
2062  { mult_spec(l1, l2, l3, r_mult(), abstract_sparse()); }
2063 
2064  // col - col matrix-matrix mult
2065 
2066  template <typename L1, typename L2, typename L3> inline
2067  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult) {
2068  mult_spec(l1, l2,l3,c_mult(),typename linalg_traits<L2>::storage_type(),
2069  typename linalg_traits<L2>::sub_orientation());
2070  }
2071 
2072 
2073  template <typename L1, typename L2, typename L3, typename ORIEN>
2074  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult,
2075  abstract_dense, ORIEN) {
2076  typedef typename linalg_traits<L2>::value_type T;
2077  size_type nn = mat_ncols(l3), mm = mat_ncols(l1);
2078 
2079  for (size_type i = 0; i < nn; ++i) {
2080  clear(mat_col(l3, i));
2081  for (size_type j = 0; j < mm; ++j) {
2082  T b = l2(j, i);
2083  if (b != T(0)) add(scaled(mat_const_col(l1, j), b), mat_col(l3, i));
2084  }
2085  }
2086  }
2087 
2088  template <typename L1, typename L2, typename L3, typename ORIEN>
2089  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult,
2090  abstract_sparse, ORIEN) {
2091  // optimizable
2092  clear(l3);
2093  size_type nn = mat_ncols(l3);
2094  for (size_type i = 0; i < nn; ++i) {
2095  typename linalg_traits<L2>::const_sub_col_type rc2 = mat_const_col(l2, i);
2096  auto it = vect_const_begin(rc2), ite = vect_const_end(rc2);
2097  for (; it != ite; ++it)
2098  add(scaled(mat_const_col(l1, it.index()), *it), mat_col(l3, i));
2099  }
2100  }
2101 
2102  template <typename L1, typename L2, typename L3>
2103  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult,
2104  abstract_sparse, row_major) {
2105  typedef typename linalg_traits<L2>::value_type T;
2106  GMM_WARNING3("Inefficient matrix-matrix mult for sparse matrices");
2107  clear(l3);
2108  size_type mm = mat_nrows(l2), nn = mat_ncols(l3);
2109  for (size_type i = 0; i < nn; ++i)
2110  for (size_type j = 0; j < mm; ++j) {
2111  T a = l2(i,j);
2112  if (a != T(0)) add(scaled(mat_const_col(l1, j), a), mat_col(l3, i));
2113  }
2114  }
2115 
2116  template <typename L1, typename L2, typename L3, typename ORIEN> inline
2117  void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult,
2118  abstract_skyline, ORIEN)
2119  { mult_spec(l1, l2, l3, c_mult(), abstract_sparse(), ORIEN()); }
2120 
2121 
2122  // col - row matrix-matrix mult
2123 
2124  template <typename L1, typename L2, typename L3> inline
2125  void mult_spec(const L1& l1, const L2& l2, L3& l3, crmult)
2126  { mult_spec(l1,l2,l3,crmult(), typename linalg_traits<L1>::storage_type()); }
2127 
2128 
2129  template <typename L1, typename L2, typename L3>
2130  void mult_spec(const L1& l1, const L2& l2, L3& l3, crmult, abstract_dense) {
2131  // optimizable
2132  clear(l3);
2133  size_type nn = mat_ncols(l1), mm = mat_nrows(l1);
2134  for (size_type i = 0; i < nn; ++i) {
2135  for (size_type j = 0; j < mm; ++j)
2136  add(scaled(mat_const_row(l2, i), l1(j, i)), mat_row(l3, j));
2137  }
2138  }
2139 
2140  template <typename L1, typename L2, typename L3>
2141  void mult_spec(const L1& l1, const L2& l2, L3& l3, crmult, abstract_sparse) {
2142  // optimizable
2143  clear(l3);
2144  size_type nn = mat_ncols(l1);
2145  for (size_type i = 0; i < nn; ++i) {
2146  typename linalg_traits<L1>::const_sub_col_type rc1 = mat_const_col(l1, i);
2147  auto it = vect_const_begin(rc1), ite = vect_const_end(rc1);
2148  for (; it != ite; ++it)
2149  add(scaled(mat_const_row(l2, i), *it), mat_row(l3, it.index()));
2150  }
2151  }
2152 
2153  template <typename L1, typename L2, typename L3> inline
2154  void mult_spec(const L1& l1, const L2& l2, L3& l3, crmult, abstract_skyline)
2155  { mult_spec(l1, l2, l3, crmult(), abstract_sparse()); }
2156 
2157 
2158  /* ******************************************************************** */
2159  /* Symmetry test. */
2160  /* ******************************************************************** */
2161 
2162  ///@endcond
2163  /** test if A is symmetric.
2164  @param A a matrix.
2165  @param tol a threshold.
2166  */
2167  template <typename MAT> inline
2168  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol
2169  = magnitude_of_linalg(MAT)(-1)) {
2170  typedef magnitude_of_linalg(MAT) R;
2171  if (tol < R(0)) tol = default_tol(R()) * mat_maxnorm(A);
2172  if (mat_nrows(A) != mat_ncols(A)) return false;
2173  return is_symmetric(A, tol, typename linalg_traits<MAT>::storage_type());
2174  }
2175  ///@cond DOXY_SHOW_ALL_FUNCTIONS
2176 
2177  template <typename MAT>
2178  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2179  abstract_dense) {
2180  size_type m = mat_nrows(A);
2181  for (size_type i = 1; i < m; ++i)
2182  for (size_type j = 0; j < i; ++j)
2183  if (gmm::abs(A(i, j)-A(j, i)) > tol) return false;
2184  return true;
2185  }
2186 
2187  template <typename MAT>
2188  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2189  abstract_sparse) {
2190  return is_symmetric(A, tol, typename principal_orientation_type<typename
2191  linalg_traits<MAT>::sub_orientation>::potype());
2192  }
2193 
2194  template <typename MAT>
2195  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2196  row_major) {
2197  for (size_type i = 0; i < mat_nrows(A); ++i) {
2198  typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, i);
2199  auto it = vect_const_begin(row), ite = vect_const_end(row);
2200  for (; it != ite; ++it)
2201  if (gmm::abs(*it - A(it.index(), i)) > tol) return false;
2202  }
2203  return true;
2204  }
2205 
2206  template <typename MAT>
2207  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2208  col_major) {
2209  for (size_type i = 0; i < mat_ncols(A); ++i) {
2210  typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, i);
2211  auto it = vect_const_begin(col), ite = vect_const_end(col);
2212  for (; it != ite; ++it)
2213  if (gmm::abs(*it - A(i, it.index())) > tol) return false;
2214  }
2215  return true;
2216  }
2217 
2218  template <typename MAT>
2219  bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol,
2220  abstract_skyline)
2221  { return is_symmetric(A, tol, abstract_sparse()); }
2222 
2223  ///@endcond
2224  /** test if A is Hermitian.
2225  @param A a matrix.
2226  @param tol a threshold.
2227  */
2228  template <typename MAT> inline
2229  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol
2230  = magnitude_of_linalg(MAT)(-1)) {
2231  typedef magnitude_of_linalg(MAT) R;
2232  if (tol < R(0)) tol = default_tol(R()) * mat_maxnorm(A);
2233  if (mat_nrows(A) != mat_ncols(A)) return false;
2234  return is_hermitian(A, tol, typename linalg_traits<MAT>::storage_type());
2235  }
2236  ///@cond DOXY_SHOW_ALL_FUNCTIONS
2237 
2238  template <typename MAT>
2239  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2240  abstract_dense) {
2241  size_type m = mat_nrows(A);
2242  for (size_type i = 1; i < m; ++i)
2243  for (size_type j = 0; j < i; ++j)
2244  if (gmm::abs(A(i, j)-gmm::conj(A(j, i))) > tol) return false;
2245  return true;
2246  }
2247 
2248  template <typename MAT>
2249  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2250  abstract_sparse) {
2251  return is_hermitian(A, tol, typename principal_orientation_type<typename
2252  linalg_traits<MAT>::sub_orientation>::potype());
2253  }
2254 
2255  template <typename MAT>
2256  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2257  row_major) {
2258  for (size_type i = 0; i < mat_nrows(A); ++i) {
2259  typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, i);
2260  auto it = vect_const_begin(row), ite = vect_const_end(row);
2261  for (; it != ite; ++it)
2262  if (gmm::abs(gmm::conj(*it) - A(it.index(), i)) > tol) return false;
2263  }
2264  return true;
2265  }
2266 
2267  template <typename MAT>
2268  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2269  col_major) {
2270  for (size_type i = 0; i < mat_ncols(A); ++i) {
2271  typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, i);
2272  auto it = vect_const_begin(col), ite = vect_const_end(col);
2273  for (; it != ite; ++it)
2274  if (gmm::abs(gmm::conj(*it) - A(i, it.index())) > tol) return false;
2275  }
2276  return true;
2277  }
2278 
2279  template <typename MAT>
2280  bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol,
2281  abstract_skyline)
2282  { return is_hermitian(A, tol, abstract_sparse()); }
2283  ///@endcond
2284 }
2285 
2286 
2287 #endif // GMM_BLAS_H__
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::is_hermitian
bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol=magnitude_of_linalg(MAT)(-1))
*‍/
Definition: gmm_blas.h:2229
gmm::fill
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*‍/
Definition: gmm_blas.h:103
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::vect_hp
strongest_value_type< V1, V2 >::value_type vect_hp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:511
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::mat_norm1
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_norm1(const M &m)
*‍/
Definition: gmm_blas.h:782
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:263
gmm::clean
void clean(L &l, double threshold)
Clean a vector or matrix (replace near-zero entries with zeroes).
gmm::is_symmetric
bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol=magnitude_of_linalg(MAT)(-1))
*‍/
Definition: gmm_blas.h:2168
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::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
gmm::mat_maxnorm
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_maxnorm(const M &m)
*‍/
Definition: gmm_blas.h:870
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_scaled.h
get a scaled view of a vector/matrix.
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::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::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:209
gmm_conjugated.h
handle conjugation of complex matrices/vectors.
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::reshape
void reshape(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:250
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_transposed.h
Generic transposed matrices.
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::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_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
gmm::mat_trace
linalg_traits< M >::value_type mat_trace(const M &m)
Trace of a matrix.
Definition: gmm_blas.h:528
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::mat_norminf
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_norminf(const M &m)
*‍/
Definition: gmm_blas.h:836
gmm::mult_add
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1781
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1268