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