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