GetFEM  5.4.3
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 <Yves.Renard@insa-lyon.fr>
34  @date October 13, 2002.
35  @brief Basic linear algebra functions.
36 */
37 
38 #ifndef GMM_BLAS_H__
39 #define GMM_BLAS_H__
40 
41 #include "gmm_scaled.h"
42 #include "gmm_transposed.h"
43 #include "gmm_conjugated.h"
44 
45 namespace gmm {
46 
47  /* ******************************************************************** */
48  /* */
49  /* Generic algorithms */
50  /* */
51  /* ******************************************************************** */
52 
53 
54  /* ******************************************************************** */
55  /* Miscellaneous */
56  /* ******************************************************************** */
57 
58  /** clear (fill with zeros) a vector or matrix. */
59  template <typename L> inline void clear(L &l)
60  { linalg_traits<L>::do_clear(l); }
61  /** @cond DOXY_SHOW_ALL_FUNCTIONS
62  skip all these redundant definitions in doxygen documentation..
63  */
64  template <typename L> inline void clear(const L &l)
65  { linalg_traits<L>::do_clear(linalg_const_cast(l)); }
66 
67  ///@endcond
68  /** count the number of non-zero entries of a vector or matrix. */
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  size_type nr = mat_nrows(ps);
299  typename linalg_traits<V2>::const_iterator
300  it = vect_const_begin(v2), ite = vect_const_end(v2);
301  typename strongest_value_type3<V1,V2,MATSP>::value_type res(0);
302  for (; it != ite; ++it)
303  res += vect_sp(mat_const_row(ps, it.index()), v1)* (*it);
304  return res;
305  }
306 
307  template <typename MATSP, typename V1, typename V2> inline
308  typename strongest_value_type3<V1,V2,MATSP>::value_type
309  vect_sp_with_matr(const MATSP &ps, const V1 &v1, const V2 &v2,
310  abstract_skyline)
311  { return vect_sp_with_matr(ps, v1, v2, abstract_sparse()); }
312 
313  template <typename MATSP, typename V1, typename V2> inline
314  typename strongest_value_type3<V1,V2,MATSP>::value_type
315  vect_sp_with_matr(const MATSP &ps, const V1 &v1, const V2 &v2,
316  abstract_dense) {
317  GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) &&
318  vect_size(v2) == mat_nrows(ps), "dimensions mismatch");
319  typename linalg_traits<V2>::const_iterator
320  it = vect_const_begin(v2), ite = vect_const_end(v2);
321  typename strongest_value_type3<V1,V2,MATSP>::value_type res(0);
322  for (size_type i = 0; it != ite; ++i, ++it)
323  res += vect_sp(mat_const_row(ps, i), v1) * (*it);
324  return res;
325  }
326 
327  template <typename MATSP, typename V1, typename V2> inline
328  typename strongest_value_type3<V1,V2,MATSP>::value_type
329  vect_sp_with_mat(const MATSP &ps, const V1 &v1,const V2 &v2,row_and_col)
330  { return vect_sp_with_mat(ps, v1, v2, row_major()); }
331 
332  template <typename MATSP, typename V1, typename V2> inline
333  typename strongest_value_type3<V1,V2,MATSP>::value_type
334  vect_sp_with_mat(const MATSP &ps, const V1 &v1, const V2 &v2,col_major){
335  return vect_sp_with_matc(ps, v1, v2,
336  typename linalg_traits<V1>::storage_type());
337  }
338 
339  template <typename MATSP, typename V1, typename V2> inline
340  typename strongest_value_type3<V1,V2,MATSP>::value_type
341  vect_sp_with_matc(const MATSP &ps, const V1 &v1, const V2 &v2,
342  abstract_sparse) {
343  GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) &&
344  vect_size(v2) == mat_nrows(ps), "dimensions mismatch");
345  typename linalg_traits<V1>::const_iterator
346  it = vect_const_begin(v1), ite = vect_const_end(v1);
347  typename strongest_value_type3<V1,V2,MATSP>::value_type res(0);
348  for (; it != ite; ++it)
349  res += vect_sp(mat_const_col(ps, it.index()), v2) * (*it);
350  return res;
351  }
352 
353  template <typename MATSP, typename V1, typename V2> inline
354  typename strongest_value_type3<V1,V2,MATSP>::value_type
355  vect_sp_with_matc(const MATSP &ps, const V1 &v1, const V2 &v2,
356  abstract_skyline)
357  { return vect_sp_with_matc(ps, v1, v2, abstract_sparse()); }
358 
359  template <typename MATSP, typename V1, typename V2> inline
360  typename strongest_value_type3<V1,V2,MATSP>::value_type
361  vect_sp_with_matc(const MATSP &ps, const V1 &v1, const V2 &v2,
362  abstract_dense) {
363  GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) &&
364  vect_size(v2) == mat_nrows(ps), "dimensions mismatch");
365  typename linalg_traits<V1>::const_iterator
366  it = vect_const_begin(v1), ite = vect_const_end(v1);
367  typename strongest_value_type3<V1,V2,MATSP>::value_type res(0);
368  for (size_type i = 0; it != ite; ++i, ++it)
369  res += vect_sp(mat_const_col(ps, i), v2) * (*it);
370  return res;
371  }
372 
373  template <typename MATSP, typename V1, typename V2> inline
374  typename strongest_value_type3<V1,V2,MATSP>::value_type
375  vect_sp_with_mat(const MATSP &ps, const V1 &v1,const V2 &v2,col_and_row)
376  { return vect_sp_with_mat(ps, v1, v2, col_major()); }
377 
378  template <typename MATSP, typename V1, typename V2> inline
379  typename strongest_value_type3<V1,V2,MATSP>::value_type
380  vect_sp_with_mat(const MATSP &ps, const V1 &v1, const V2 &v2,
381  abstract_null_type) {
382  typename temporary_vector<V1>::vector_type w(mat_nrows(ps));
383  GMM_WARNING2("Warning, a temporary is used in scalar product\n");
384  mult(ps, v1, w);
385  return vect_sp(w, v2);
386  }
387 
388  template <typename IT1, typename IT2> inline
389  typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
390  typename std::iterator_traits<IT2>::value_type>::T
391  vect_sp_dense_(IT1 it, IT1 ite, IT2 it2) {
392  typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
393  typename std::iterator_traits<IT2>::value_type>::T res(0);
394  for (; it != ite; ++it, ++it2) res += (*it) * (*it2);
395  return res;
396  }
397 
398  template <typename IT1, typename V> inline
399  typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
400  typename linalg_traits<V>::value_type>::T
401  vect_sp_sparse_(IT1 it, IT1 ite, const V &v) {
402  typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
403  typename linalg_traits<V>::value_type>::T res(0);
404  for (; it != ite; ++it) res += (*it) * v[it.index()];
405  return res;
406  }
407 
408  template <typename V1, typename V2> inline
409  typename strongest_value_type<V1,V2>::value_type
410  vect_sp(const V1 &v1, const V2 &v2, abstract_dense, abstract_dense) {
411  return vect_sp_dense_(vect_const_begin(v1), vect_const_end(v1),
412  vect_const_begin(v2));
413  }
414 
415  template <typename V1, typename V2> inline
416  typename strongest_value_type<V1,V2>::value_type
417  vect_sp(const V1 &v1, const V2 &v2, abstract_skyline, abstract_dense) {
418  typename linalg_traits<V1>::const_iterator it1 = vect_const_begin(v1),
419  ite = vect_const_end(v1);
420  typename linalg_traits<V2>::const_iterator it2 = vect_const_begin(v2);
421  return vect_sp_dense_(it1, ite, it2 + it1.index());
422  }
423 
424  template <typename V1, typename V2> inline
425  typename strongest_value_type<V1,V2>::value_type
426  vect_sp(const V1 &v1, const V2 &v2, abstract_dense, abstract_skyline) {
427  typename linalg_traits<V2>::const_iterator it1 = vect_const_begin(v2),
428  ite = vect_const_end(v2);
429  typename linalg_traits<V1>::const_iterator it2 = vect_const_begin(v1);
430  return vect_sp_dense_(it1, ite, it2 + it1.index());
431  }
432 
433  template <typename V1, typename V2> inline
434  typename strongest_value_type<V1,V2>::value_type
435  vect_sp(const V1 &v1, const V2 &v2, abstract_skyline, abstract_skyline) {
436  typedef typename strongest_value_type<V1,V2>::value_type T;
437  auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
438  auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
439  size_type n = std::min(ite1.index(), ite2.index());
440  size_type l = std::max(it1.index(), it2.index());
441 
442  if (l < n) {
443  size_type m = l - it1.index(), p = l - it2.index(), q = m + n - l;
444  return vect_sp_dense_(it1+m, it1+q, it2 + p);
445  }
446  return T(0);
447  }
448 
449  template <typename V1, typename V2> inline
450  typename strongest_value_type<V1,V2>::value_type
451  vect_sp(const V1 &v1, const V2 &v2,abstract_sparse,abstract_dense) {
452  return vect_sp_sparse_(vect_const_begin(v1), vect_const_end(v1), v2);
453  }
454 
455  template <typename V1, typename V2> inline
456  typename strongest_value_type<V1,V2>::value_type
457  vect_sp(const V1 &v1, const V2 &v2, abstract_sparse, abstract_skyline) {
458  return vect_sp_sparse_(vect_const_begin(v1), vect_const_end(v1), v2);
459  }
460 
461  template <typename V1, typename V2> inline
462  typename strongest_value_type<V1,V2>::value_type
463  vect_sp(const V1 &v1, const V2 &v2, abstract_skyline, abstract_sparse) {
464  return vect_sp_sparse_(vect_const_begin(v2), vect_const_end(v2), v1);
465  }
466 
467  template <typename V1, typename V2> inline
468  typename strongest_value_type<V1,V2>::value_type
469  vect_sp(const V1 &v1, const V2 &v2, abstract_dense,abstract_sparse) {
470  return vect_sp_sparse_(vect_const_begin(v2), vect_const_end(v2), v1);
471  }
472 
473 
474  template <typename V1, typename V2> inline
475  typename strongest_value_type<V1,V2>::value_type
476  vect_sp_sparse_sparse(const V1 &v1, const V2 &v2, linalg_true) {
477  typename linalg_traits<V1>::const_iterator it1 = vect_const_begin(v1),
478  ite1 = vect_const_end(v1);
479  typename linalg_traits<V2>::const_iterator it2 = vect_const_begin(v2),
480  ite2 = vect_const_end(v2);
481  typename strongest_value_type<V1,V2>::value_type res(0);
482 
483  while (it1 != ite1 && it2 != ite2) {
484  if (it1.index() == it2.index())
485  { res += (*it1) * *it2; ++it1; ++it2; }
486  else if (it1.index() < it2.index()) ++it1; else ++it2;
487  }
488  return res;
489  }
490 
491  template <typename V1, typename V2> inline
492  typename strongest_value_type<V1,V2>::value_type
493  vect_sp_sparse_sparse(const V1 &v1, const V2 &v2, linalg_false) {
494  return vect_sp_sparse_(vect_const_begin(v1), vect_const_end(v1), v2);
495  }
496 
497  template <typename V1, typename V2> inline
498  typename strongest_value_type<V1,V2>::value_type
499  vect_sp(const V1 &v1, const V2 &v2,abstract_sparse,abstract_sparse) {
500  return vect_sp_sparse_sparse(v1, v2,
501  typename linalg_and<typename linalg_traits<V1>::index_sorted,
502  typename linalg_traits<V2>::index_sorted>::bool_type());
503  }
504 
505  /* ******************************************************************** */
506  /* Hermitian product */
507  /* ******************************************************************** */
508  ///@endcond
509  /** Hermitian product. */
510  template <typename V1, typename V2>
511  inline typename strongest_value_type<V1,V2>::value_type
512  vect_hp(const V1 &v1, const V2 &v2)
513  { return vect_sp(v1, conjugated(v2)); }
514 
515  /** Hermitian product with a matrix. */
516  template <typename MATSP, typename V1, typename V2> inline
517  typename strongest_value_type3<V1,V2,MATSP>::value_type
518  vect_hp(const MATSP &ps, const V1 &v1, const V2 &v2) {
519  return vect_sp(ps, v1, gmm::conjugated(v2));
520  }
521 
522  /* ******************************************************************** */
523  /* Trace of a matrix */
524  /* ******************************************************************** */
525 
526  /** Trace of a matrix */
527  template <typename M>
528  typename linalg_traits<M>::value_type
529  mat_trace(const M &m) {
530  typedef typename linalg_traits<M>::value_type T;
531  T res(0);
532  for (size_type i = 0; i < std::min(mat_nrows(m), mat_ncols(m)); ++i)
533  res += m(i,i);
534  return res;
535  }
536 
537  /* ******************************************************************** */
538  /* Euclidean norm */
539  /* ******************************************************************** */
540 
541  /** squared Euclidean norm of a vector. */
542  template <typename V>
543  typename number_traits<typename linalg_traits<V>::value_type>
544  ::magnitude_type
545  vect_norm2_sqr(const V &v) {
546  typedef typename linalg_traits<V>::value_type T;
547  typedef typename number_traits<T>::magnitude_type R;
548  auto it = vect_const_begin(v), ite = vect_const_end(v);
549  R res(0);
550  for (; it != ite; ++it) res += gmm::abs_sqr(*it);
551  return res;
552  }
553 
554  /** Euclidean norm of a vector. */
555  template <typename V> inline
556  typename number_traits<typename linalg_traits<V>::value_type>
557  ::magnitude_type
558  vect_norm2(const V &v)
559  { return sqrt(vect_norm2_sqr(v)); }
560 
561 
562  /** squared Euclidean distance between two vectors */
563  template <typename V1, typename V2> inline
564  typename number_traits<typename linalg_traits<V1>::value_type>
565  ::magnitude_type
566  vect_dist2_sqr(const V1 &v1, const V2 &v2) { // not fully optimized
567  typedef typename linalg_traits<V1>::value_type T;
568  typedef typename number_traits<T>::magnitude_type R;
569  auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
570  auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
571  size_type k1(0), k2(0);
572  R res(0);
573  while (it1 != ite1 && it2 != ite2) {
574  size_type i1 = index_of_it(it1, k1,
575  typename linalg_traits<V1>::storage_type());
576  size_type i2 = index_of_it(it2, k2,
577  typename linalg_traits<V2>::storage_type());
578 
579  if (i1 == i2) {
580  res += gmm::abs_sqr(*it2 - *it1); ++it1; ++k1; ++it2; ++k2;
581  }
582  else if (i1 < i2) {
583  res += gmm::abs_sqr(*it1); ++it1; ++k1;
584  }
585  else {
586  res += gmm::abs_sqr(*it2); ++it2; ++k2;
587  }
588  }
589  while (it1 != ite1) { res += gmm::abs_sqr(*it1); ++it1; }
590  while (it2 != ite2) { res += gmm::abs_sqr(*it2); ++it2; }
591  return res;
592  }
593 
594  /** Euclidean distance between two vectors */
595  template <typename V1, typename V2> inline
596  typename number_traits<typename linalg_traits<V1>::value_type>
597  ::magnitude_type
598  vect_dist2(const V1 &v1, const V2 &v2)
599  { return sqrt(vect_dist2_sqr(v1, v2)); }
600  ///@cond DOXY_SHOW_ALL_FUNCTIONS
601  template <typename M>
602  typename number_traits<typename linalg_traits<M>::value_type>
603  ::magnitude_type
604  mat_euclidean_norm_sqr(const M &m, row_major) {
605  typename number_traits<typename linalg_traits<M>::value_type>
606  ::magnitude_type res(0);
607  for (size_type i = 0; i < mat_nrows(m); ++i)
608  res += vect_norm2_sqr(mat_const_row(m, i));
609  return res;
610  }
611 
612  template <typename M>
613  typename number_traits<typename linalg_traits<M>::value_type>
614  ::magnitude_type
615  mat_euclidean_norm_sqr(const M &m, col_major) {
616  typename number_traits<typename linalg_traits<M>::value_type>
617  ::magnitude_type res(0);
618  for (size_type i = 0; i < mat_ncols(m); ++i)
619  res += vect_norm2_sqr(mat_const_col(m, i));
620  return res;
621  }
622  ///@endcond
623  /** squared Euclidean norm of a matrix. */
624  template <typename M> inline
625  typename number_traits<typename linalg_traits<M>::value_type>
626  ::magnitude_type
627  mat_euclidean_norm_sqr(const M &m) {
628  return mat_euclidean_norm_sqr(m,
629  typename principal_orientation_type<typename
630  linalg_traits<M>::sub_orientation>::potype());
631  }
632 
633  /** Euclidean norm of a matrix. */
634  template <typename M> inline
635  typename number_traits<typename linalg_traits<M>::value_type>
636  ::magnitude_type
637  mat_euclidean_norm(const M &m)
638  { return gmm::sqrt(mat_euclidean_norm_sqr(m)); }
639 
640  /* ******************************************************************** */
641  /* vector norm1 */
642  /* ******************************************************************** */
643  /** 1-norm of a vector */
644  template <typename V>
645  typename number_traits<typename linalg_traits<V>::value_type>
646  ::magnitude_type
647  vect_norm1(const V &v) {
648  auto it = vect_const_begin(v), ite = vect_const_end(v);
649  typename number_traits<typename linalg_traits<V>::value_type>
650  ::magnitude_type res(0);
651  for (; it != ite; ++it) res += gmm::abs(*it);
652  return res;
653  }
654 
655  /** 1-distance between two vectors */
656  template <typename V1, typename V2> inline
657  typename number_traits<typename linalg_traits<V1>::value_type>
658  ::magnitude_type
659  vect_dist1(const V1 &v1, const V2 &v2) { // not fully optimized
660  typedef typename linalg_traits<V1>::value_type T;
661  typedef typename number_traits<T>::magnitude_type R;
662  auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
663  auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
664  size_type k1(0), k2(0);
665  R res(0);
666  while (it1 != ite1 && it2 != ite2) {
667  size_type i1 = index_of_it(it1, k1,
668  typename linalg_traits<V1>::storage_type());
669  size_type i2 = index_of_it(it2, k2,
670  typename linalg_traits<V2>::storage_type());
671 
672  if (i1 == i2) {
673  res += gmm::abs(*it2 - *it1); ++it1; ++k1; ++it2; ++k2;
674  }
675  else if (i1 < i2) {
676  res += gmm::abs(*it1); ++it1; ++k1;
677  }
678  else {
679  res += gmm::abs(*it2); ++it2; ++k2;
680  }
681  }
682  while (it1 != ite1) { res += gmm::abs(*it1); ++it1; }
683  while (it2 != ite2) { res += gmm::abs(*it2); ++it2; }
684  return res;
685  }
686 
687  /* ******************************************************************** */
688  /* vector Infinity norm */
689  /* ******************************************************************** */
690  /** Infinity norm of a vector. */
691  template <typename V>
692  typename number_traits<typename linalg_traits<V>::value_type>
693  ::magnitude_type
694  vect_norminf(const V &v) {
695  auto it = vect_const_begin(v), ite = vect_const_end(v);
696  typename number_traits<typename linalg_traits<V>::value_type>
697  ::magnitude_type res(0);
698  for (; it != ite; ++it) res = std::max(res, gmm::abs(*it));
699  return res;
700  }
701 
702  /** Infinity distance between two vectors */
703  template <typename V1, typename V2> inline
704  typename number_traits<typename linalg_traits<V1>::value_type>
705  ::magnitude_type
706  vect_distinf(const V1 &v1, const V2 &v2) { // not fully optimized
707  typedef typename linalg_traits<V1>::value_type T;
708  typedef typename number_traits<T>::magnitude_type R;
709  auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
710  auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
711  size_type k1(0), k2(0);
712  R res(0);
713  while (it1 != ite1 && it2 != ite2) {
714  size_type i1 = index_of_it(it1, k1,
715  typename linalg_traits<V1>::storage_type());
716  size_type i2 = index_of_it(it2, k2,
717  typename linalg_traits<V2>::storage_type());
718 
719  if (i1 == i2) {
720  res = std::max(res, gmm::abs(*it2 - *it1)); ++it1; ++k1; ++it2; ++k2;
721  }
722  else if (i1 < i2) {
723  res = std::max(res, gmm::abs(*it1)); ++it1; ++k1;
724  }
725  else {
726  res = std::max(res, gmm::abs(*it2)); ++it2; ++k2;
727  }
728  }
729  while (it1 != ite1) { res = std::max(res, gmm::abs(*it1)); ++it1; }
730  while (it2 != ite2) { res = std::max(res, gmm::abs(*it2)); ++it2; }
731  return res;
732  }
733 
734  /* ******************************************************************** */
735  /* matrix norm1 */
736  /* ******************************************************************** */
737  ///@cond DOXY_SHOW_ALL_FUNCTIONS
738  template <typename M>
739  typename number_traits<typename linalg_traits<M>::value_type>
740  ::magnitude_type
741  mat_norm1(const M &m, col_major) {
742  typename number_traits<typename linalg_traits<M>::value_type>
743  ::magnitude_type res(0);
744  for (size_type i = 0; i < mat_ncols(m); ++i)
745  res = std::max(res, vect_norm1(mat_const_col(m,i)));
746  return res;
747  }
748 
749  template <typename M>
750  typename number_traits<typename linalg_traits<M>::value_type>
751  ::magnitude_type
752  mat_norm1(const M &m, row_major) {
753  typedef typename linalg_traits<M>::value_type T;
754  typedef typename number_traits<T>::magnitude_type R;
755  typedef typename linalg_traits<M>::storage_type store_type;
756 
757  std::vector<R> aux(mat_ncols(m));
758  for (size_type i = 0; i < mat_nrows(m); ++i) {
759  typename linalg_traits<M>::const_sub_row_type row = mat_const_row(m, i);
760  auto it = vect_const_begin(row), ite = vect_const_end(row);
761  for (size_type k = 0; it != ite; ++it, ++k)
762  aux[index_of_it(it, k, store_type())] += gmm::abs(*it);
763  }
764  return vect_norminf(aux);
765  }
766 
767  template <typename M>
768  typename number_traits<typename linalg_traits<M>::value_type>
769  ::magnitude_type
770  mat_norm1(const M &m, col_and_row)
771  { return mat_norm1(m, col_major()); }
772 
773  template <typename M>
774  typename number_traits<typename linalg_traits<M>::value_type>
775  ::magnitude_type
776  mat_norm1(const M &m, row_and_col)
777  { return mat_norm1(m, col_major()); }
778  ///@endcond
779  /** 1-norm of a matrix */
780  template <typename M>
781  typename number_traits<typename linalg_traits<M>::value_type>
782  ::magnitude_type
783  mat_norm1(const M &m) {
784  return mat_norm1(m, typename linalg_traits<M>::sub_orientation());
785  }
786 
787 
788  /* ******************************************************************** */
789  /* matrix Infinity norm */
790  /* ******************************************************************** */
791  ///@cond DOXY_SHOW_ALL_FUNCTIONS
792  template <typename M>
793  typename number_traits<typename linalg_traits<M>::value_type>
794  ::magnitude_type
795  mat_norminf(const M &m, row_major) {
796  typename number_traits<typename linalg_traits<M>::value_type>
797  ::magnitude_type res(0);
798  for (size_type i = 0; i < mat_nrows(m); ++i)
799  res = std::max(res, vect_norm1(mat_const_row(m,i)));
800  return res;
801  }
802 
803  template <typename M>
804  typename number_traits<typename linalg_traits<M>::value_type>
805  ::magnitude_type
806  mat_norminf(const M &m, col_major) {
807  typedef typename linalg_traits<M>::value_type T;
808  typedef typename number_traits<T>::magnitude_type R;
809  typedef typename linalg_traits<M>::storage_type store_type;
810 
811  std::vector<R> aux(mat_nrows(m));
812  for (size_type i = 0; i < mat_ncols(m); ++i) {
813  typename linalg_traits<M>::const_sub_col_type col = mat_const_col(m, i);
814  auto it = vect_const_begin(col), ite = vect_const_end(col);
815  for (size_type k = 0; it != ite; ++it, ++k)
816  aux[index_of_it(it, k, store_type())] += gmm::abs(*it);
817  }
818  return vect_norminf(aux);
819  }
820 
821  template <typename M>
822  typename number_traits<typename linalg_traits<M>::value_type>
823  ::magnitude_type
824  mat_norminf(const M &m, col_and_row)
825  { return mat_norminf(m, row_major()); }
826 
827  template <typename M>
828  typename number_traits<typename linalg_traits<M>::value_type>
829  ::magnitude_type
830  mat_norminf(const M &m, row_and_col)
831  { return mat_norminf(m, row_major()); }
832  ///@endcond
833  /** infinity-norm of a matrix.*/
834  template <typename M>
835  typename number_traits<typename linalg_traits<M>::value_type>
836  ::magnitude_type
837  mat_norminf(const M &m) {
838  return mat_norminf(m, typename linalg_traits<M>::sub_orientation());
839  }
840 
841  /* ******************************************************************** */
842  /* Max norm for matrices */
843  /* ******************************************************************** */
844  ///@cond DOXY_SHOW_ALL_FUNCTIONS
845  template <typename M>
846  typename number_traits<typename linalg_traits<M>::value_type>
847  ::magnitude_type
848  mat_maxnorm(const M &m, row_major) {
849  typename number_traits<typename linalg_traits<M>::value_type>
850  ::magnitude_type res(0);
851  for (size_type i = 0; i < mat_nrows(m); ++i)
852  res = std::max(res, vect_norminf(mat_const_row(m,i)));
853  return res;
854  }
855 
856  template <typename M>
857  typename number_traits<typename linalg_traits<M>::value_type>
858  ::magnitude_type
859  mat_maxnorm(const M &m, col_major) {
860  typename number_traits<typename linalg_traits<M>::value_type>
861  ::magnitude_type res(0);
862  for (size_type i = 0; i < mat_ncols(m); ++i)
863  res = std::max(res, vect_norminf(mat_const_col(m,i)));
864  return res;
865  }
866  ///@endcond
867  /** max-norm of a matrix. */
868  template <typename M>
869  typename number_traits<typename linalg_traits<M>::value_type>
870  ::magnitude_type
871  mat_maxnorm(const M &m) {
872  return mat_maxnorm
873  (m, typename principal_orientation_type
874  <typename linalg_traits<M>::sub_orientation>::potype());
875  }
876 
877  /* ******************************************************************** */
878  /* Clean */
879  /* ******************************************************************** */
880  /** Clean a vector or matrix (replace near-zero entries with zeroes). */
881 
882  template <typename L> inline void clean(L &l, double threshold);
883 
884  ///@cond DOXY_SHOW_ALL_FUNCTIONS
885 
886  template <typename L, typename T>
887  void clean(L &l, double threshold, abstract_dense, T) {
888  typedef typename number_traits<T>::magnitude_type R;
889  auto it = vect_begin(l), ite = vect_end(l);
890  for (; it != ite; ++it)
891  if (gmm::abs(*it) < R(threshold)) *it = T(0);
892  }
893 
894  template <typename L, typename T>
895  void clean(L &l, double threshold, abstract_skyline, T)
896  { gmm::clean(l, threshold, abstract_dense(), T()); }
897 
898  template <typename L, typename T>
899  void clean(L &l, double threshold, abstract_sparse, T) {
900  typedef typename number_traits<T>::magnitude_type R;
901  auto it = vect_begin(l), ite = vect_end(l);
902  std::vector<size_type> ind;
903  for (; it != ite; ++it)
904  if (gmm::abs(*it) < R(threshold)) ind.push_back(it.index());
905  for (size_type i = 0; i < ind.size(); ++i) l[ind[i]] = T(0);
906  }
907 
908  template <typename L, typename T>
909  void clean(L &l, double threshold, abstract_dense, std::complex<T>) {
910  auto it = vect_begin(l), ite = vect_end(l);
911  for (; it != ite; ++it){
912  if (gmm::abs((*it).real()) < T(threshold))
913  *it = std::complex<T>(T(0), (*it).imag());
914  if (gmm::abs((*it).imag()) < T(threshold))
915  *it = std::complex<T>((*it).real(), T(0));
916  }
917  }
918 
919  template <typename L, typename T>
920  void clean(L &l, double threshold, abstract_skyline, std::complex<T>)
921  { gmm::clean(l, threshold, abstract_dense(), std::complex<T>()); }
922 
923  template <typename L, typename T>
924  void clean(L &l, double threshold, abstract_sparse, std::complex<T>) {
925  auto it = vect_begin(l), ite = vect_end(l);
926  std::vector<size_type> ind;
927  for (; it != ite; ++it) {
928  bool r = (gmm::abs((*it).real()) < T(threshold));
929  bool i = (gmm::abs((*it).imag()) < T(threshold));
930  if (r && i) ind.push_back(it.index());
931  else if (r) *it = std::complex<T>(T(0), (*it).imag());
932  else if (i) *it = std::complex<T>((*it).real(), T(0));
933  }
934  for (size_type i = 0; i < ind.size(); ++i)
935  l[ind[i]] = std::complex<T>(T(0),T(0));
936  }
937 
938  template <typename L> inline void clean(L &l, double threshold,
939  abstract_vector) {
940  gmm::clean(l, threshold, typename linalg_traits<L>::storage_type(),
941  typename linalg_traits<L>::value_type());
942  }
943 
944  template <typename L> inline void clean(const L &l, double threshold);
945 
946  template <typename L> void clean(L &l, double threshold, row_major) {
947  for (size_type i = 0; i < mat_nrows(l); ++i)
948  gmm::clean(mat_row(l, i), threshold);
949  }
950 
951  template <typename L> void clean(L &l, double threshold, col_major) {
952  for (size_type i = 0; i < mat_ncols(l); ++i)
953  gmm::clean(mat_col(l, i), threshold);
954  }
955 
956  template <typename L> inline void clean(L &l, double threshold,
957  abstract_matrix) {
958  gmm::clean(l, threshold,
959  typename principal_orientation_type<typename
960  linalg_traits<L>::sub_orientation>::potype());
961  }
962 
963  template <typename L> inline void clean(L &l, double threshold)
964  { clean(l, threshold, typename linalg_traits<L>::linalg_type()); }
965 
966  template <typename L> inline void clean(const L &l, double threshold)
967  { gmm::clean(linalg_const_cast(l), threshold); }
968 
969  /* ******************************************************************** */
970  /* Copy */
971  /* ******************************************************************** */
972  ///@endcond
973  /** Copy vectors or matrices.
974  @param l1 source vector or matrix.
975  @param l2 destination.
976  */
977  template <typename L1, typename L2> inline
978  void copy(const L1& l1, L2& l2) {
979  if ((const void *)(&l1) != (const void *)(&l2)) {
980  if (same_origin(l1,l2))
981  GMM_WARNING2("Warning : a conflict is possible in copy\n");
982 
983  copy(l1, l2, typename linalg_traits<L1>::linalg_type(),
984  typename linalg_traits<L2>::linalg_type());
985  }
986  }
987  ///@cond DOXY_SHOW_ALL_FUNCTIONS
988 
989  template <typename L1, typename L2> inline
990  void copy(const L1& l1, const L2& l2) { copy(l1, linalg_const_cast(l2)); }
991 
992  template <typename L1, typename L2> inline
993  void copy(const L1& l1, L2& l2, abstract_vector, abstract_vector) {
994  GMM_ASSERT2(vect_size(l1) == vect_size(l2), "dimensions mismatch, "
995  << vect_size(l1) << " !=" << vect_size(l2));
996  copy_vect(l1, l2, typename linalg_traits<L1>::storage_type(),
997  typename linalg_traits<L2>::storage_type());
998  }
999 
1000  template <typename L1, typename L2> inline
1001  void copy(const L1& l1, L2& l2, abstract_matrix, abstract_matrix) {
1002  size_type m = mat_nrows(l1), n = mat_ncols(l1);
1003  if (!m || !n) return;
1004  GMM_ASSERT2(n==mat_ncols(l2) && m==mat_nrows(l2), "dimensions mismatch");
1005  copy_mat(l1, l2, typename linalg_traits<L1>::sub_orientation(),
1006  typename linalg_traits<L2>::sub_orientation());
1007  }
1008 
1009  template <typename V1, typename V2, typename C1, typename C2> inline
1010  void copy_vect(const V1 &v1, const V2 &v2, C1, C2)
1011  { copy_vect(v1, const_cast<V2 &>(v2), C1(), C2()); }
1012 
1013 
1014  template <typename L1, typename L2>
1015  void copy_mat_by_row(const L1& l1, L2& l2) {
1016  size_type nbr = mat_nrows(l1);
1017  for (size_type i = 0; i < nbr; ++i)
1018  copy(mat_const_row(l1, i), mat_row(l2, i));
1019  }
1020 
1021  template <typename L1, typename L2>
1022  void copy_mat_by_col(const L1 &l1, L2 &l2) {
1023  size_type nbc = mat_ncols(l1);
1024  for (size_type i = 0; i < nbc; ++i) {
1025  copy(mat_const_col(l1, i), mat_col(l2, i));
1026  }
1027  }
1028 
1029  template <typename L1, typename L2> inline
1030  void copy_mat(const L1& l1, L2& l2, row_major, row_major)
1031  { copy_mat_by_row(l1, l2); }
1032 
1033  template <typename L1, typename L2> inline
1034  void copy_mat(const L1& l1, L2& l2, row_major, row_and_col)
1035  { copy_mat_by_row(l1, l2); }
1036 
1037  template <typename L1, typename L2> inline
1038  void copy_mat(const L1& l1, L2& l2, row_and_col, row_and_col)
1039  { copy_mat_by_row(l1, l2); }
1040 
1041  template <typename L1, typename L2> inline
1042  void copy_mat(const L1& l1, L2& l2, row_and_col, row_major)
1043  { copy_mat_by_row(l1, l2); }
1044 
1045  template <typename L1, typename L2> inline
1046  void copy_mat(const L1& l1, L2& l2, col_and_row, row_major)
1047  { copy_mat_by_row(l1, l2); }
1048 
1049  template <typename L1, typename L2> inline
1050  void copy_mat(const L1& l1, L2& l2, row_major, col_and_row)
1051  { copy_mat_by_row(l1, l2); }
1052 
1053  template <typename L1, typename L2> inline
1054  void copy_mat(const L1& l1, L2& l2, col_and_row, row_and_col)
1055  { copy_mat_by_row(l1, l2); }
1056 
1057  template <typename L1, typename L2> inline
1058  void copy_mat(const L1& l1, L2& l2, row_and_col, col_and_row)
1059  { copy_mat_by_row(l1, l2); }
1060 
1061  template <typename L1, typename L2> inline
1062  void copy_mat(const L1& l1, L2& l2, col_major, col_major)
1063  { copy_mat_by_col(l1, l2); }
1064 
1065  template <typename L1, typename L2> inline
1066  void copy_mat(const L1& l1, L2& l2, col_major, col_and_row)
1067  { copy_mat_by_col(l1, l2); }
1068 
1069  template <typename L1, typename L2> inline
1070  void copy_mat(const L1& l1, L2& l2, col_major, row_and_col)
1071  { copy_mat_by_col(l1, l2); }
1072 
1073  template <typename L1, typename L2> inline
1074  void copy_mat(const L1& l1, L2& l2, row_and_col, col_major)
1075  { copy_mat_by_col(l1, l2); }
1076 
1077  template <typename L1, typename L2> inline
1078  void copy_mat(const L1& l1, L2& l2, col_and_row, col_major)
1079  { copy_mat_by_col(l1, l2); }
1080 
1081  template <typename L1, typename L2> inline
1082  void copy_mat(const L1& l1, L2& l2, col_and_row, col_and_row)
1083  { copy_mat_by_col(l1, l2); }
1084 
1085  template <typename L1, typename L2> inline
1086  void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i) {
1087  copy_mat_mixed_rc(l1, l2, i, typename linalg_traits<L1>::storage_type());
1088  }
1089 
1090  template <typename L1, typename L2>
1091  void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_sparse) {
1092  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1093  for (; it != ite; ++it)
1094  l2(i, it.index()) = *it;
1095  }
1096 
1097  template <typename L1, typename L2>
1098  void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_skyline) {
1099  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1100  for (; it != ite; ++it)
1101  l2(i, it.index()) = *it;
1102  }
1103 
1104  template <typename L1, typename L2>
1105  void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_dense) {
1106  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1107  for (size_type j = 0; it != ite; ++it, ++j) l2(i, j) = *it;
1108  }
1109 
1110  template <typename L1, typename L2> inline
1111  void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i) {
1112  copy_mat_mixed_cr(l1, l2, i, typename linalg_traits<L1>::storage_type());
1113  }
1114 
1115  template <typename L1, typename L2>
1116  void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_sparse) {
1117  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1118  for (; it != ite; ++it) l2(it.index(), i) = *it;
1119  }
1120 
1121  template <typename L1, typename L2>
1122  void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_skyline) {
1123  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1124  for (; it != ite; ++it) l2(it.index(), i) = *it;
1125  }
1126 
1127  template <typename L1, typename L2>
1128  void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_dense) {
1129  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1130  for (size_type j = 0; it != ite; ++it, ++j) l2(j, i) = *it;
1131  }
1132 
1133  template <typename L1, typename L2>
1134  void copy_mat(const L1& l1, L2& l2, row_major, col_major) {
1135  clear(l2);
1136  size_type nbr = mat_nrows(l1);
1137  for (size_type i = 0; i < nbr; ++i)
1138  copy_mat_mixed_rc(mat_const_row(l1, i), l2, i);
1139  }
1140 
1141  template <typename L1, typename L2>
1142  void copy_mat(const L1& l1, L2& l2, col_major, row_major) {
1143  clear(l2);
1144  size_type nbc = mat_ncols(l1);
1145  for (size_type i = 0; i < nbc; ++i)
1146  copy_mat_mixed_cr(mat_const_col(l1, i), l2, i);
1147  }
1148 
1149  template <typename L1, typename L2> inline
1150  void copy_vect(const L1 &l1, L2 &l2, abstract_dense, abstract_dense) {
1151  std::copy(vect_const_begin(l1), vect_const_end(l1), vect_begin(l2));
1152  }
1153 
1154  template <typename L1, typename L2> inline // to be optimised ?
1155  void copy_vect(const L1 &l1, L2 &l2, abstract_skyline, abstract_skyline) {
1156  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1157  while (it1 != ite1 && *it1 == typename linalg_traits<L1>::value_type(0))
1158  ++it1;
1159 
1160  if (ite1 - it1 > 0) {
1161  clear(l2);
1162  auto it2 = vect_begin(l2), ite2 = vect_end(l2);
1163  while (*(ite1-1) == typename linalg_traits<L1>::value_type(0))
1164  ite1--;
1165 
1166  if (it2 == ite2) {
1167  l2[it1.index()] = *it1;
1168  ++it1;
1169  l2[ite1.index()-1] = *(ite1-1);
1170  --ite1;
1171  if (it1 < ite1) {
1172  it2 = vect_begin(l2);
1173  ++it2;
1174  std::copy(it1, ite1, it2);
1175  }
1176  } else {
1177  ptrdiff_t m = it1.index() - it2.index();
1178  if (m >= 0 && ite1.index() <= ite2.index())
1179  std::copy(it1, ite1, it2 + m);
1180  else {
1181  if (m < 0)
1182  l2[it1.index()] = *it1;
1183  if (ite1.index() > ite2.index())
1184  l2[ite1.index()-1] = *(ite1-1);
1185  it2 = vect_begin(l2);
1186  ite2 = vect_end(l2);
1187  m = it1.index() - it2.index();
1188  std::copy(it1, ite1, it2 + m);
1189  }
1190  }
1191  }
1192  }
1193 
1194  template <typename L1, typename L2>
1195  void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_dense) {
1196  clear(l2);
1197  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1198  for (; it != ite; ++it) { l2[it.index()] = *it; }
1199  }
1200 
1201  template <typename L1, typename L2>
1202  void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_skyline) {
1203  clear(l2);
1204  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1205  for (; it != ite; ++it) l2[it.index()] = *it;
1206  }
1207 
1208  template <typename L1, typename L2>
1209  void copy_vect(const L1& l1, L2& l2, abstract_skyline, abstract_dense) {
1210  typedef typename linalg_traits<L1>::value_type T;
1211  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1212  if (it == ite)
1213  gmm::clear(l2);
1214  else {
1215  auto it2 = vect_begin(l2), ite2 = vect_end(l2);
1216 
1217  size_type i = it.index(), j;
1218  for (j = 0; j < i; ++j, ++it2) *it2 = T(0);
1219  for (; it != ite; ++it, ++it2) *it2 = *it;
1220  for (; it2 != ite2; ++it2) *it2 = T(0);
1221  }
1222  }
1223 
1224  template <typename L1, typename L2>
1225  void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_sparse) {
1226  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1227  clear(l2);
1228  // cout << "copy " << l1 << " of size " << vect_size(l1) << " nnz = " << nnz(l1) << endl;
1229  for (; it != ite; ++it) {
1230  // cout << "*it = " << *it << endl;
1231  // cout << "it.index() = " << it.index() << endl;
1232  if (*it != (typename linalg_traits<L1>::value_type)(0))
1233  l2[it.index()] = *it;
1234  }
1235  }
1236 
1237  template <typename L1, typename L2>
1238  void copy_vect(const L1& l1, L2& l2, abstract_dense, abstract_sparse) {
1239  clear(l2);
1240  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1241  for (size_type i = 0; it != ite; ++it, ++i)
1242  if (*it != (typename linalg_traits<L1>::value_type)(0))
1243  l2[i] = *it;
1244  }
1245 
1246  template <typename L1, typename L2> // to be optimised ...
1247  void copy_vect(const L1& l1, L2& l2, abstract_dense, abstract_skyline) {
1248  clear(l2);
1249  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1250  for (size_type i = 0; it != ite; ++it, ++i)
1251  if (*it != (typename linalg_traits<L1>::value_type)(0))
1252  l2[i] = *it;
1253  }
1254 
1255 
1256  template <typename L1, typename L2>
1257  void copy_vect(const L1& l1, L2& l2, abstract_skyline, abstract_sparse) {
1258  clear(l2);
1259  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1260  for (; it != ite; ++it)
1261  if (*it != (typename linalg_traits<L1>::value_type)(0))
1262  l2[it.index()] = *it;
1263  }
1264 
1265  /* ******************************************************************** */
1266  /* Matrix and vector addition */
1267  /* algorithms are built in order to avoid some conflicts with */
1268  /* repeated arguments or with overlapping part of a same object. */
1269  /* In the latter case, conflicts are still possible. */
1270  /* ******************************************************************** */
1271  ///@endcond
1272  /** Add two vectors or matrices
1273  @param l1
1274  @param l2 contains on output, l2+l1.
1275  */
1276  template <typename L1, typename L2> inline
1277  void add(const L1& l1, L2& l2) {
1278  add_spec(l1, l2, typename linalg_traits<L2>::linalg_type());
1279  }
1280  ///@cond
1281 
1282  template <typename L1, typename L2> inline
1283  void add(const L1& l1, const L2& l2) { add(l1, linalg_const_cast(l2)); }
1284 
1285  template <typename L1, typename L2> inline
1286  void add_spec(const L1& l1, L2& l2, abstract_vector) {
1287  GMM_ASSERT2(vect_size(l1) == vect_size(l2), "dimensions mismatch, "
1288  << vect_size(l1) << " !=" << vect_size(l2));
1289  add(l1, l2, typename linalg_traits<L1>::storage_type(),
1290  typename linalg_traits<L2>::storage_type());
1291  }
1292 
1293  template <typename L1, typename L2> inline
1294  void add_spec(const L1& l1, L2& l2, abstract_matrix) {
1295  GMM_ASSERT2(mat_nrows(l1)==mat_nrows(l2) && mat_ncols(l1)==mat_ncols(l2),
1296  "dimensions mismatch l1 is " << mat_nrows(l1) << "x"
1297  << mat_ncols(l1) << " and l2 is " << mat_nrows(l2)
1298  << "x" << mat_ncols(l2));
1299  add(l1, l2, typename linalg_traits<L1>::sub_orientation(),
1300  typename linalg_traits<L2>::sub_orientation());
1301  }
1302 
1303  template <typename L1, typename L2>
1304  void add(const L1& l1, L2& l2, row_major, row_major) {
1305  auto it1 = mat_row_begin(l1), ite = mat_row_end(l1);
1306  auto it2 = mat_row_begin(l2);
1307  for ( ; it1 != ite; ++it1, ++it2)
1308  add(linalg_traits<L1>::row(it1), linalg_traits<L2>::row(it2));
1309  }
1310 
1311  template <typename L1, typename L2>
1312  void add(const L1& l1, L2& l2, col_major, col_major) {
1313  auto it1 = mat_col_const_begin(l1), ite = mat_col_const_end(l1);
1314  typename linalg_traits<L2>::col_iterator it2 = mat_col_begin(l2);
1315  for ( ; it1 != ite; ++it1, ++it2)
1316  add(linalg_traits<L1>::col(it1), linalg_traits<L2>::col(it2));
1317  }
1318 
1319  template <typename L1, typename L2> inline
1320  void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i) {
1321  add_mat_mixed_rc(l1, l2, i, typename linalg_traits<L1>::storage_type());
1322  }
1323 
1324  template <typename L1, typename L2>
1325  void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_sparse) {
1326  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1327  for (; it != ite; ++it) l2(i, it.index()) += *it;
1328  }
1329 
1330  template <typename L1, typename L2>
1331  void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_skyline) {
1332  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1333  for (; it != ite; ++it) l2(i, it.index()) += *it;
1334  }
1335 
1336  template <typename L1, typename L2>
1337  void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_dense) {
1338  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1339  for (size_type j = 0; it != ite; ++it, ++j) l2(i, j) += *it;
1340  }
1341 
1342  template <typename L1, typename L2> inline
1343  void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i) {
1344  add_mat_mixed_cr(l1, l2, i, typename linalg_traits<L1>::storage_type());
1345  }
1346 
1347  template <typename L1, typename L2>
1348  void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_sparse) {
1349  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1350  for (; it != ite; ++it) l2(it.index(), i) += *it;
1351  }
1352 
1353  template <typename L1, typename L2>
1354  void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_skyline) {
1355  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1356  for (; it != ite; ++it) l2(it.index(), i) += *it;
1357  }
1358 
1359  template <typename L1, typename L2>
1360  void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_dense) {
1361  auto it = vect_const_begin(l1), ite = vect_const_end(l1);
1362  for (size_type j = 0; it != ite; ++it, ++j) l2(j, i) += *it;
1363  }
1364 
1365  template <typename L1, typename L2>
1366  void add(const L1& l1, L2& l2, row_major, col_major) {
1367  size_type nbr = mat_nrows(l1);
1368  for (size_type i = 0; i < nbr; ++i)
1369  add_mat_mixed_rc(mat_const_row(l1, i), l2, i);
1370  }
1371 
1372  template <typename L1, typename L2>
1373  void add(const L1& l1, L2& l2, col_major, row_major) {
1374  size_type nbc = mat_ncols(l1);
1375  for (size_type i = 0; i < nbc; ++i)
1376  add_mat_mixed_cr(mat_const_col(l1, i), l2, i);
1377  }
1378 
1379  template <typename L1, typename L2> inline
1380  void add(const L1& l1, L2& l2, row_and_col, row_major)
1381  { add(l1, l2, row_major(), row_major()); }
1382 
1383  template <typename L1, typename L2> inline
1384  void add(const L1& l1, L2& l2, row_and_col, row_and_col)
1385  { add(l1, l2, row_major(), row_major()); }
1386 
1387  template <typename L1, typename L2> inline
1388  void add(const L1& l1, L2& l2, row_and_col, col_and_row)
1389  { add(l1, l2, row_major(), row_major()); }
1390 
1391  template <typename L1, typename L2> inline
1392  void add(const L1& l1, L2& l2, col_and_row, row_and_col)
1393  { add(l1, l2, row_major(), row_major()); }
1394 
1395  template <typename L1, typename L2> inline
1396  void add(const L1& l1, L2& l2, row_major, row_and_col)
1397  { add(l1, l2, row_major(), row_major()); }
1398 
1399  template <typename L1, typename L2> inline
1400  void add(const L1& l1, L2& l2, col_and_row, row_major)
1401  { add(l1, l2, row_major(), row_major()); }
1402 
1403  template <typename L1, typename L2> inline
1404  void add(const L1& l1, L2& l2, row_major, col_and_row)
1405  { add(l1, l2, row_major(), row_major()); }
1406 
1407  template <typename L1, typename L2> inline
1408  void add(const L1& l1, L2& l2, row_and_col, col_major)
1409  { add(l1, l2, col_major(), col_major()); }
1410 
1411  template <typename L1, typename L2> inline
1412  void add(const L1& l1, L2& l2, col_major, row_and_col)
1413  { add(l1, l2, col_major(), col_major()); }
1414 
1415  template <typename L1, typename L2> inline
1416  void add(const L1& l1, L2& l2, col_and_row, col_major)
1417  { add(l1, l2, col_major(), col_major()); }
1418 
1419  template <typename L1, typename L2> inline
1420  void add(const L1& l1, L2& l2, col_and_row, col_and_row)
1421  { add(l1, l2, col_major(), col_major()); }
1422 
1423  template <typename L1, typename L2> inline
1424  void add(const L1& l1, L2& l2, col_major, col_and_row)
1425  { add(l1, l2, col_major(), col_major()); }
1426 
1427  ///@endcond
1428  /** Addition of two vectors/matrices
1429  @param l1
1430  @param l2
1431  @param l3 contains l1+l2 on output
1432  */
1433  template <typename L1, typename L2, typename L3> inline
1434  void add(const L1& l1, const L2& l2, L3& l3) {
1435  add_spec(l1, l2, l3, typename linalg_traits<L2>::linalg_type());
1436  }
1437  ///@cond DOXY_SHOW_ALL_FUNCTIONS
1438 
1439  template <typename L1, typename L2, typename L3> inline
1440  void add(const L1& l1, const L2& l2, const L3& l3)
1441  { add(l1, l2, linalg_const_cast(l3)); }
1442 
1443  template <typename L1, typename L2, typename L3> inline
1444  void add_spec(const L1& l1, const L2& l2, L3& l3, abstract_matrix)
1445  { copy(l2, l3); add(l1, l3); }
1446 
1447  template <typename L1, typename L2, typename L3> inline
1448  void add_spec(const L1& l1, const L2& l2, L3& l3, abstract_vector) {
1449  GMM_ASSERT2(vect_size(l1) == vect_size(l2), "dimensions mismatch, "
1450  << vect_size(l1) << " !=" << vect_size(l2));
1451  GMM_ASSERT2(vect_size(l1) == vect_size(l3), "dimensions mismatch, "
1452  << vect_size(l1) << " !=" << vect_size(l3));
1453  if ((const void *)(&l1) == (const void *)(&l3))
1454  add(l2, l3);
1455  else if ((const void *)(&l2) == (const void *)(&l3))
1456  add(l1, l3);
1457  else
1458  add(l1, l2, l3, typename linalg_traits<L1>::storage_type(),
1459  typename linalg_traits<L2>::storage_type(),
1460  typename linalg_traits<L3>::storage_type());
1461  }
1462 
1463  template <typename IT1, typename IT2, typename IT3>
1464  void add_full_(IT1 it1, IT2 it2, IT3 it3, IT3 ite) {
1465  for (; it3 != ite; ++it3, ++it2, ++it1) *it3 = *it1 + *it2;
1466  }
1467 
1468  template <typename IT1, typename IT2, typename IT3>
1469  void add_almost_full_(IT1 it1, IT1 ite1, IT2 it2, IT3 it3, IT3 ite3) {
1470  IT3 it = it3;
1471  for (; it != ite3; ++it, ++it2) *it = *it2;
1472  for (; it1 != ite1; ++it1)
1473  *(it3 + it1.index()) += *it1;
1474  }
1475 
1476  template <typename IT1, typename IT2, typename IT3>
1477  void add_to_full_(IT1 it1, IT1 ite1, IT2 it2, IT2 ite2,
1478  IT3 it3, IT3 ite3) {
1479  typedef typename std::iterator_traits<IT3>::value_type T;
1480  IT3 it = it3;
1481  for (; it != ite3; ++it) *it = T(0);
1482  for (; it1 != ite1; ++it1) *(it3 + it1.index()) = *it1;
1483  for (; it2 != ite2; ++it2) *(it3 + it2.index()) += *it2;
1484  }
1485 
1486  template <typename L1, typename L2, typename L3> inline
1487  void add(const L1& l1, const L2& l2, L3& l3,
1488  abstract_dense, abstract_dense, abstract_dense) {
1489  add_full_(vect_const_begin(l1), vect_const_begin(l2),
1490  vect_begin(l3), vect_end(l3));
1491  }
1492 
1493  // generic function for add(v1, v2, v3).
1494  // Need to be specialized to optimize particular additions.
1495  template <typename L1, typename L2, typename L3,
1496  typename ST1, typename ST2, typename ST3>
1497  inline void add(const L1& l1, const L2& l2, L3& l3, ST1, ST2, ST3)
1498  { copy(l2, l3); add(l1, l3, ST1(), ST3()); }
1499 
1500  template <typename L1, typename L2, typename L3> inline
1501  void add(const L1& l1, const L2& l2, L3& l3,
1502  abstract_sparse, abstract_dense, abstract_dense) {
1503  add_almost_full_(vect_const_begin(l1), vect_const_end(l1),
1504  vect_const_begin(l2), vect_begin(l3), vect_end(l3));
1505  }
1506 
1507  template <typename L1, typename L2, typename L3> inline
1508  void add(const L1& l1, const L2& l2, L3& l3,
1509  abstract_dense, abstract_sparse, abstract_dense)
1510  { add(l2, l1, l3, abstract_sparse(), abstract_dense(), abstract_dense()); }
1511 
1512  template <typename L1, typename L2, typename L3> inline
1513  void add(const L1& l1, const L2& l2, L3& l3,
1514  abstract_sparse, abstract_sparse, abstract_dense) {
1515  add_to_full_(vect_const_begin(l1), vect_const_end(l1),
1516  vect_const_begin(l2), vect_const_end(l2),
1517  vect_begin(l3), vect_end(l3));
1518  }
1519 
1520 
1521  template <typename L1, typename L2, typename L3>
1522  void add_spspsp(const L1& l1, const L2& l2, L3& l3, linalg_true) {
1523  auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1);
1524  auto it2 = vect_const_begin(l2), ite2 = vect_const_end(l2);
1525  clear(l3);
1526  while (it1 != ite1 && it2 != ite2) {
1527  ptrdiff_t d = it1.index() - it2.index();
1528  if (d < 0)
1529  { l3[it1.index()] += *it1; ++it1; }
1530  else if (d > 0)
1531  { l3[it2.index()] += *it2; ++it2; }
1532  else
1533  { l3[it1.index()] = *it1 + *it2; ++it1; ++it2; }
1534  }
1535  for (; it1 != ite1; ++it1) l3[it1.index()] += *it1;
1536  for (; it2 != ite2; ++it2) l3[it2.index()] += *it2;
1537  }
1538 
1539  template <typename L1, typename L2, typename L3>
1540  void add_spspsp(const L1& l1, const L2& l2, L3& l3, linalg_false)
1541  { copy(l2, l3); add(l2, l3); }
1542 
1543  template <typename L1, typename L2, typename L3>
1544  void add(const L1& l1, const L2& l2, L3& l3,
1545  abstract_sparse, abstract_sparse, abstract_sparse) {
1546  add_spspsp(l1, l2, l3, typename linalg_and<typename
1547  linalg_traits<L1>::index_sorted,
1548  typename linalg_traits<L2>::index_sorted>::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:978
gmm::mat_norm1
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_norm1(const M &m)
*‍/
Definition: gmm_blas.h:783
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:558
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:659
gmm::mat_maxnorm
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_maxnorm(const M &m)
*‍/
Definition: gmm_blas.h:871
gmm::vect_hp
strongest_value_type< V1, V2 >::value_type vect_hp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:512
gmm::mat_trace
linalg_traits< M >::value_type mat_trace(const M &m)
Trace of a matrix.
Definition: gmm_blas.h:529
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:637
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:706
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:566
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:545
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:694
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:598
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:647
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:1277
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:627
gmm::mat_norminf
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_norminf(const M &m)
*‍/
Definition: gmm_blas.h:837
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