|
|
|
@ -1,11 +1,11 @@ |
|
|
|
/* -*- c++ -*- (enables emacs c++ mode) */ |
|
|
|
/*=========================================================================== |
|
|
|
|
|
|
|
Copyright (C) 2002-2015 Yves Renard |
|
|
|
|
|
|
|
This file is a part of GETFEM++ |
|
|
|
|
|
|
|
Getfem++ is free software; you can redistribute it and/or modify it |
|
|
|
|
|
|
|
Copyright (C) 2002-2017 Yves Renard |
|
|
|
|
|
|
|
This file is a part of GetFEM++ |
|
|
|
|
|
|
|
GetFEM++ is free software; you can redistribute it and/or modify it |
|
|
|
under the terms of the GNU Lesser General Public License as published |
|
|
|
by the Free Software Foundation; either version 3 of the License, or |
|
|
|
(at your option) any later version along with the GCC Runtime Library |
|
|
|
@ -17,7 +17,7 @@ |
|
|
|
You should have received a copy of the GNU Lesser General Public License |
|
|
|
along with this program; if not, write to the Free Software Foundation, |
|
|
|
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. |
|
|
|
|
|
|
|
|
|
|
|
As a special exception, you may use this file as it is a part of a free |
|
|
|
software library without restriction. Specifically, if other files |
|
|
|
instantiate templates or use macros or inline functions from this file, |
|
|
|
@ -26,7 +26,7 @@ |
|
|
|
to be covered by the GNU Lesser General Public License. This exception |
|
|
|
does not however invalidate any other reasons why the executable file |
|
|
|
might be covered by the GNU Lesser General Public License. |
|
|
|
|
|
|
|
|
|
|
|
===========================================================================*/ |
|
|
|
|
|
|
|
/**@file gmm_blas.h |
|
|
|
@ -38,16 +38,17 @@ |
|
|
|
#ifndef GMM_BLAS_H__ |
|
|
|
#define GMM_BLAS_H__ |
|
|
|
|
|
|
|
// This Version of GMM was modified for StoRM. |
|
|
|
// To detect whether the usage of TBB is possible, this include is neccessary |
|
|
|
// This Version of gmm++ was modified for Storm. |
|
|
|
// To detect whether the usage of TBB is possible, this include is neccessary. |
|
|
|
#include "storm-config.h" |
|
|
|
|
|
|
|
#ifdef STORM_HAVE_INTELTBB |
|
|
|
# include <new> // This fixes a potential dependency ordering problem between GMM and TBB |
|
|
|
# include "tbb/tbb.h" |
|
|
|
# include <iterator> |
|
|
|
#include <new> // This fixes a potential dependency ordering problem between GMM and TBB |
|
|
|
#include "tbb/tbb.h" |
|
|
|
#include <iterator> |
|
|
|
#endif |
|
|
|
|
|
|
|
|
|
|
|
#include "gmm_scaled.h" |
|
|
|
#include "gmm_transposed.h" |
|
|
|
#include "gmm_conjugated.h" |
|
|
|
@ -79,9 +80,8 @@ namespace gmm { |
|
|
|
{ return nnz(l, typename linalg_traits<L>::linalg_type()); } |
|
|
|
|
|
|
|
///@cond DOXY_SHOW_ALL_FUNCTIONS |
|
|
|
template <typename L> inline size_type nnz(const L& l, abstract_vector) { |
|
|
|
typename linalg_traits<L>::const_iterator it = vect_const_begin(l), |
|
|
|
ite = vect_const_end(l); |
|
|
|
template <typename L> inline size_type nnz(const L& l, abstract_vector) { |
|
|
|
auto it = vect_const_begin(l), ite = vect_const_end(l); |
|
|
|
size_type res(0); |
|
|
|
for (; it != ite; ++it) ++res; |
|
|
|
return res; |
|
|
|
@ -444,10 +444,8 @@ namespace gmm { |
|
|
|
typename strongest_value_type<V1,V2>::value_type |
|
|
|
vect_sp(const V1 &v1, const V2 &v2, abstract_skyline, abstract_skyline) { |
|
|
|
typedef typename strongest_value_type<V1,V2>::value_type T; |
|
|
|
typename linalg_traits<V1>::const_iterator it1 = vect_const_begin(v1), |
|
|
|
ite1 = vect_const_end(v1); |
|
|
|
typename linalg_traits<V2>::const_iterator it2 = vect_const_begin(v2), |
|
|
|
ite2 = vect_const_end(v2); |
|
|
|
auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1); |
|
|
|
auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2); |
|
|
|
size_type n = std::min(ite1.index(), ite2.index()); |
|
|
|
size_type l = std::max(it1.index(), it2.index()); |
|
|
|
|
|
|
|
@ -557,8 +555,7 @@ namespace gmm { |
|
|
|
vect_norm2_sqr(const V &v) { |
|
|
|
typedef typename linalg_traits<V>::value_type T; |
|
|
|
typedef typename number_traits<T>::magnitude_type R; |
|
|
|
typename linalg_traits<V>::const_iterator |
|
|
|
it = vect_const_begin(v), ite = vect_const_end(v); |
|
|
|
auto it = vect_const_begin(v), ite = vect_const_end(v); |
|
|
|
R res(0); |
|
|
|
for (; it != ite; ++it) res += gmm::abs_sqr(*it); |
|
|
|
return res; |
|
|
|
@ -579,10 +576,8 @@ namespace gmm { |
|
|
|
vect_dist2_sqr(const V1 &v1, const V2 &v2) { // not fully optimized |
|
|
|
typedef typename linalg_traits<V1>::value_type T; |
|
|
|
typedef typename number_traits<T>::magnitude_type R; |
|
|
|
typename linalg_traits<V1>::const_iterator |
|
|
|
it1 = vect_const_begin(v1), ite1 = vect_const_end(v1); |
|
|
|
typename linalg_traits<V2>::const_iterator |
|
|
|
it2 = vect_const_begin(v2), ite2 = vect_const_end(v2); |
|
|
|
auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1); |
|
|
|
auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2); |
|
|
|
size_type k1(0), k2(0); |
|
|
|
R res(0); |
|
|
|
while (it1 != ite1 && it2 != ite2) { |
|
|
|
@ -660,8 +655,7 @@ namespace gmm { |
|
|
|
typename number_traits<typename linalg_traits<V>::value_type> |
|
|
|
::magnitude_type |
|
|
|
vect_norm1(const V &v) { |
|
|
|
typename linalg_traits<V>::const_iterator |
|
|
|
it = vect_const_begin(v), ite = vect_const_end(v); |
|
|
|
auto it = vect_const_begin(v), ite = vect_const_end(v); |
|
|
|
typename number_traits<typename linalg_traits<V>::value_type> |
|
|
|
::magnitude_type res(0); |
|
|
|
for (; it != ite; ++it) res += gmm::abs(*it); |
|
|
|
@ -676,10 +670,9 @@ namespace gmm { |
|
|
|
typename number_traits<typename linalg_traits<V>::value_type> |
|
|
|
::magnitude_type |
|
|
|
vect_norminf(const V &v) { |
|
|
|
typename linalg_traits<V>::const_iterator |
|
|
|
it = vect_const_begin(v), ite = vect_const_end(v); |
|
|
|
typename number_traits<typename linalg_traits<V>::value_type> |
|
|
|
::magnitude_type res(0); |
|
|
|
auto it = vect_const_begin(v), ite = vect_const_end(v); |
|
|
|
typename number_traits<typename linalg_traits<V>::value_type> |
|
|
|
::magnitude_type res(0); |
|
|
|
for (; it != ite; ++it) res = std::max(res, gmm::abs(*it)); |
|
|
|
return res; |
|
|
|
} |
|
|
|
@ -709,10 +702,8 @@ namespace gmm { |
|
|
|
|
|
|
|
std::vector<R> aux(mat_ncols(m)); |
|
|
|
for (size_type i = 0; i < mat_nrows(m); ++i) { |
|
|
|
typedef typename linalg_traits<M>::const_sub_row_type row_type; |
|
|
|
row_type row = mat_const_row(m, i); |
|
|
|
typename linalg_traits<row_type>::const_iterator |
|
|
|
it = vect_const_begin(row), ite = vect_const_end(row); |
|
|
|
typename linalg_traits<M>::const_sub_row_type row = mat_const_row(m, i); |
|
|
|
auto it = vect_const_begin(row), ite = vect_const_end(row); |
|
|
|
for (size_type k = 0; it != ite; ++it, ++k) |
|
|
|
aux[index_of_it(it, k, store_type())] += gmm::abs(*it); |
|
|
|
} |
|
|
|
@ -765,10 +756,8 @@ namespace gmm { |
|
|
|
|
|
|
|
std::vector<R> aux(mat_nrows(m)); |
|
|
|
for (size_type i = 0; i < mat_ncols(m); ++i) { |
|
|
|
typedef typename linalg_traits<M>::const_sub_col_type col_type; |
|
|
|
col_type col = mat_const_col(m, i); |
|
|
|
typename linalg_traits<col_type>::const_iterator |
|
|
|
it = vect_const_begin(col), ite = vect_const_end(col); |
|
|
|
typename linalg_traits<M>::const_sub_col_type col = mat_const_col(m, i); |
|
|
|
auto it = vect_const_begin(col), ite = vect_const_end(col); |
|
|
|
for (size_type k = 0; it != ite; ++it, ++k) |
|
|
|
aux[index_of_it(it, k, store_type())] += gmm::abs(*it); |
|
|
|
} |
|
|
|
@ -843,7 +832,7 @@ namespace gmm { |
|
|
|
template <typename L, typename T> |
|
|
|
void clean(L &l, double threshold, abstract_dense, T) { |
|
|
|
typedef typename number_traits<T>::magnitude_type R; |
|
|
|
typename linalg_traits<L>::iterator it = vect_begin(l), ite = vect_end(l); |
|
|
|
auto it = vect_begin(l), ite = vect_end(l); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (gmm::abs(*it) < R(threshold)) *it = T(0); |
|
|
|
} |
|
|
|
@ -855,7 +844,7 @@ namespace gmm { |
|
|
|
template <typename L, typename T> |
|
|
|
void clean(L &l, double threshold, abstract_sparse, T) { |
|
|
|
typedef typename number_traits<T>::magnitude_type R; |
|
|
|
typename linalg_traits<L>::iterator it = vect_begin(l), ite = vect_end(l); |
|
|
|
auto it = vect_begin(l), ite = vect_end(l); |
|
|
|
std::vector<size_type> ind; |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (gmm::abs(*it) < R(threshold)) ind.push_back(it.index()); |
|
|
|
@ -864,7 +853,7 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L, typename T> |
|
|
|
void clean(L &l, double threshold, abstract_dense, std::complex<T>) { |
|
|
|
typename linalg_traits<L>::iterator it = vect_begin(l), ite = vect_end(l); |
|
|
|
auto it = vect_begin(l), ite = vect_end(l); |
|
|
|
for (; it != ite; ++it){ |
|
|
|
if (gmm::abs((*it).real()) < T(threshold)) |
|
|
|
*it = std::complex<T>(T(0), (*it).imag()); |
|
|
|
@ -879,7 +868,7 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L, typename T> |
|
|
|
void clean(L &l, double threshold, abstract_sparse, std::complex<T>) { |
|
|
|
typename linalg_traits<L>::iterator it = vect_begin(l), ite = vect_end(l); |
|
|
|
auto it = vect_begin(l), ite = vect_end(l); |
|
|
|
std::vector<size_type> ind; |
|
|
|
for (; it != ite; ++it) { |
|
|
|
bool r = (gmm::abs((*it).real()) < T(threshold)); |
|
|
|
@ -972,18 +961,14 @@ namespace gmm { |
|
|
|
void copy_mat_by_row(const L1& l1, L2& l2) { |
|
|
|
size_type nbr = mat_nrows(l1); |
|
|
|
for (size_type i = 0; i < nbr; ++i) |
|
|
|
copy_vect(mat_const_row(l1, i), mat_row(l2, i), |
|
|
|
typename linalg_traits<L1>::storage_type(), |
|
|
|
typename linalg_traits<L2>::storage_type()); |
|
|
|
copy(mat_const_row(l1, i), mat_row(l2, i)); |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_mat_by_col(const L1 &l1, L2 &l2) { |
|
|
|
size_type nbc = mat_ncols(l1); |
|
|
|
for (size_type i = 0; i < nbc; ++i) { |
|
|
|
copy_vect(mat_const_col(l1, i), mat_col(l2, i), |
|
|
|
typename linalg_traits<L1>::storage_type(), |
|
|
|
typename linalg_traits<L2>::storage_type()); |
|
|
|
copy(mat_const_col(l1, i), mat_col(l2, i)); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
@ -1050,24 +1035,21 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_sparse) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) |
|
|
|
l2(i, it.index()) = *it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_skyline) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) |
|
|
|
l2(i, it.index()) = *it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_dense) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (size_type j = 0; it != ite; ++it, ++j) l2(i, j) = *it; |
|
|
|
} |
|
|
|
|
|
|
|
@ -1078,22 +1060,19 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_sparse) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) l2(it.index(), i) = *it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_skyline) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) l2(it.index(), i) = *it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_dense) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (size_type j = 0; it != ite; ++it, ++j) l2(j, i) = *it; |
|
|
|
} |
|
|
|
|
|
|
|
@ -1120,15 +1099,13 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> inline // to be optimised ? |
|
|
|
void copy_vect(const L1 &l1, L2 &l2, abstract_skyline, abstract_skyline) { |
|
|
|
typename linalg_traits<L1>::const_iterator it1 = vect_const_begin(l1), |
|
|
|
ite1 = vect_const_end(l1); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
while (it1 != ite1 && *it1 == typename linalg_traits<L1>::value_type(0)) |
|
|
|
++it1; |
|
|
|
|
|
|
|
if (ite1 - it1 > 0) { |
|
|
|
clear(l2); |
|
|
|
typename linalg_traits<L2>::iterator it2 = vect_begin(l2), |
|
|
|
ite2 = vect_end(l2); |
|
|
|
auto it2 = vect_begin(l2), ite2 = vect_end(l2); |
|
|
|
while (*(ite1-1) == typename linalg_traits<L1>::value_type(0)) ite1--; |
|
|
|
|
|
|
|
if (it2 == ite2) { |
|
|
|
@ -1155,29 +1132,25 @@ namespace gmm { |
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_dense) { |
|
|
|
clear(l2); |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) { l2[it.index()] = *it; } |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_skyline) { |
|
|
|
clear(l2); |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) l2[it.index()] = *it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_vect(const L1& l1, L2& l2, abstract_skyline, abstract_dense) { |
|
|
|
typedef typename linalg_traits<L1>::value_type T; |
|
|
|
typedef typename linalg_traits<L1>::const_iterator l1_const_iterator; |
|
|
|
typedef typename linalg_traits<L2>::iterator l2_iterator; |
|
|
|
l1_const_iterator it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
if (it == ite) |
|
|
|
gmm::clear(l2); |
|
|
|
else { |
|
|
|
l2_iterator it2 = vect_begin(l2), ite2 = vect_end(l2); |
|
|
|
auto it2 = vect_begin(l2), ite2 = vect_end(l2); |
|
|
|
|
|
|
|
size_type i = it.index(), j; |
|
|
|
for (j = 0; j < i; ++j, ++it2) *it2 = T(0); |
|
|
|
@ -1188,19 +1161,21 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_sparse) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
clear(l2); |
|
|
|
for (; it != ite; ++it) |
|
|
|
// cout << "copy " << l1 << " of size " << vect_size(l1) << " nnz = " << nnz(l1) << endl; |
|
|
|
for (; it != ite; ++it) { |
|
|
|
// cout << "*it = " << *it << endl; |
|
|
|
// cout << "it.index() = " << it.index() << endl; |
|
|
|
if (*it != (typename linalg_traits<L1>::value_type)(0)) |
|
|
|
l2[it.index()] = *it; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_vect(const L1& l1, L2& l2, abstract_dense, abstract_sparse) { |
|
|
|
clear(l2); |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (size_type i = 0; it != ite; ++it, ++i) |
|
|
|
if (*it != (typename linalg_traits<L1>::value_type)(0)) |
|
|
|
l2[i] = *it; |
|
|
|
@ -1209,8 +1184,7 @@ namespace gmm { |
|
|
|
template <typename L1, typename L2> // to be optimised ... |
|
|
|
void copy_vect(const L1& l1, L2& l2, abstract_dense, abstract_skyline) { |
|
|
|
clear(l2); |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (size_type i = 0; it != ite; ++it, ++i) |
|
|
|
if (*it != (typename linalg_traits<L1>::value_type)(0)) |
|
|
|
l2[i] = *it; |
|
|
|
@ -1220,8 +1194,7 @@ namespace gmm { |
|
|
|
template <typename L1, typename L2> |
|
|
|
void copy_vect(const L1& l1, L2& l2, abstract_skyline, abstract_sparse) { |
|
|
|
clear(l2); |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (*it != (typename linalg_traits<L1>::value_type)(0)) |
|
|
|
l2[it.index()] = *it; |
|
|
|
@ -1258,25 +1231,24 @@ namespace gmm { |
|
|
|
template <typename L1, typename L2> inline |
|
|
|
void add_spec(const L1& l1, L2& l2, abstract_matrix) { |
|
|
|
GMM_ASSERT2(mat_nrows(l1)==mat_nrows(l2) && mat_ncols(l1)==mat_ncols(l2), |
|
|
|
"dimensions mismatch"); |
|
|
|
"dimensions mismatch l1 is " << mat_nrows(l1) << "x" |
|
|
|
<< mat_ncols(l1) << " and l2 is " << mat_nrows(l2) |
|
|
|
<< "x" << mat_ncols(l2)); |
|
|
|
add(l1, l2, typename linalg_traits<L1>::sub_orientation(), |
|
|
|
typename linalg_traits<L2>::sub_orientation()); |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, row_major, row_major) { |
|
|
|
typename linalg_traits<L1>::const_row_iterator it1 = mat_row_begin(l1), |
|
|
|
ite = mat_row_end(l1); |
|
|
|
typename linalg_traits<L2>::row_iterator it2 = mat_row_begin(l2); |
|
|
|
auto it1 = mat_row_begin(l1), ite = mat_row_end(l1); |
|
|
|
auto it2 = mat_row_begin(l2); |
|
|
|
for ( ; it1 != ite; ++it1, ++it2) |
|
|
|
add(linalg_traits<L1>::row(it1), linalg_traits<L2>::row(it2)); |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, col_major, col_major) { |
|
|
|
typename linalg_traits<L1>::const_col_iterator |
|
|
|
it1 = mat_col_const_begin(l1), |
|
|
|
ite = mat_col_const_end(l1); |
|
|
|
auto it1 = mat_col_const_begin(l1), ite = mat_col_const_end(l1); |
|
|
|
typename linalg_traits<L2>::col_iterator it2 = mat_col_begin(l2); |
|
|
|
for ( ; it1 != ite; ++it1, ++it2) |
|
|
|
add(linalg_traits<L1>::col(it1), linalg_traits<L2>::col(it2)); |
|
|
|
@ -1289,22 +1261,19 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_sparse) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) l2(i, it.index()) += *it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_skyline) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) l2(i, it.index()) += *it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add_mat_mixed_rc(const L1& l1, L2& l2, size_type i, abstract_dense) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (size_type j = 0; it != ite; ++it, ++j) l2(i, j) += *it; |
|
|
|
} |
|
|
|
|
|
|
|
@ -1315,22 +1284,19 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_sparse) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) l2(it.index(), i) += *it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_skyline) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (; it != ite; ++it) l2(it.index(), i) += *it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add_mat_mixed_cr(const L1& l1, L2& l2, size_type i, abstract_dense) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
auto it = vect_const_begin(l1), ite = vect_const_end(l1); |
|
|
|
for (size_type j = 0; it != ite; ++it, ++j) l2(j, i) += *it; |
|
|
|
} |
|
|
|
|
|
|
|
@ -1490,10 +1456,8 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2, typename L3> |
|
|
|
void add_spspsp(const L1& l1, const L2& l2, L3& l3, linalg_true) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
typename linalg_traits<L2>::const_iterator |
|
|
|
it2 = vect_const_begin(l2), ite2 = vect_const_end(l2); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
auto it2 = vect_const_begin(l2), ite2 = vect_const_end(l2); |
|
|
|
clear(l3); |
|
|
|
while (it1 != ite1 && it2 != ite2) { |
|
|
|
ptrdiff_t d = it1.index() - it2.index(); |
|
|
|
@ -1522,23 +1486,20 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, abstract_dense, abstract_dense) { |
|
|
|
typename linalg_traits<L1>::const_iterator it1 = vect_const_begin(l1); |
|
|
|
typename linalg_traits<L2>::iterator |
|
|
|
it2 = vect_begin(l2), ite = vect_end(l2); |
|
|
|
auto it1 = vect_const_begin(l1); |
|
|
|
auto it2 = vect_begin(l2), ite = vect_end(l2); |
|
|
|
for (; it2 != ite; ++it2, ++it1) *it2 += *it1; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, abstract_dense, abstract_skyline) { |
|
|
|
typedef typename linalg_traits<L1>::const_iterator const_l1_iterator; |
|
|
|
typedef typename linalg_traits<L2>::iterator l2_iterator; |
|
|
|
typedef typename linalg_traits<L1>::value_type T; |
|
|
|
|
|
|
|
const_l1_iterator it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
size_type i1 = 0, ie1 = vect_size(l1); |
|
|
|
while (it1 != ite1 && *it1 == T(0)) { ++it1; ++i1; } |
|
|
|
if (it1 != ite1) { |
|
|
|
l2_iterator it2 = vect_begin(l2), ite2 = vect_end(l2); |
|
|
|
auto it2 = vect_begin(l2), ite2 = vect_end(l2); |
|
|
|
while (ie1 && *(ite1-1) == T(0)) { ite1--; --ie1; } |
|
|
|
|
|
|
|
if (it2 == ite2 || i1 < it2.index()) { |
|
|
|
@ -1558,10 +1519,9 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, abstract_skyline, abstract_dense) { |
|
|
|
typename linalg_traits<L1>::const_iterator it1 = vect_const_begin(l1), |
|
|
|
ite1 = vect_const_end(l1); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
if (it1 != ite1) { |
|
|
|
typename linalg_traits<L2>::iterator it2 = vect_begin(l2); |
|
|
|
auto it2 = vect_begin(l2); |
|
|
|
it2 += it1.index(); |
|
|
|
for (; it1 != ite1; ++it2, ++it1) *it2 += *it1; |
|
|
|
} |
|
|
|
@ -1570,30 +1530,26 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, abstract_sparse, abstract_dense) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
for (; it1 != ite1; ++it1) l2[it1.index()] += *it1; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, abstract_sparse, abstract_sparse) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
for (; it1 != ite1; ++it1) l2[it1.index()] += *it1; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, abstract_sparse, abstract_skyline) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
for (; it1 != ite1; ++it1) l2[it1.index()] += *it1; |
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, abstract_skyline, abstract_sparse) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
for (; it1 != ite1; ++it1) |
|
|
|
if (*it1 != typename linalg_traits<L1>::value_type(0)) |
|
|
|
l2[it1.index()] += *it1; |
|
|
|
@ -1601,16 +1557,14 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, abstract_skyline, abstract_skyline) { |
|
|
|
typedef typename linalg_traits<L1>::const_iterator const_l1_iterator; |
|
|
|
typedef typename linalg_traits<L2>::iterator l2_iterator; |
|
|
|
typedef typename linalg_traits<L1>::value_type T1; |
|
|
|
typedef typename linalg_traits<L2>::value_type T2; |
|
|
|
|
|
|
|
const_l1_iterator it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
|
|
|
|
while (it1 != ite1 && *it1 == T1(0)) ++it1; |
|
|
|
if (ite1 != it1) { |
|
|
|
l2_iterator it2 = vect_begin(l2), ite2 = vect_end(l2); |
|
|
|
auto it2 = vect_begin(l2), ite2 = vect_end(l2); |
|
|
|
while (*(ite1-1) == T1(0)) ite1--; |
|
|
|
if (it2 == ite2 || it1.index() < it2.index()) { |
|
|
|
l2[it1.index()] = T2(0); |
|
|
|
@ -1627,8 +1581,7 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2> |
|
|
|
void add(const L1& l1, L2& l2, abstract_dense, abstract_sparse) { |
|
|
|
typename linalg_traits<L1>::const_iterator |
|
|
|
it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); |
|
|
|
for (size_type i = 0; it1 != ite1; ++it1, ++i) |
|
|
|
if (*it1 != typename linalg_traits<L1>::value_type(0)) l2[i] += *it1; |
|
|
|
} |
|
|
|
@ -1690,7 +1643,7 @@ namespace gmm { |
|
|
|
if (aux != T(0)) l3[i] = aux; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
#ifdef STORM_HAVE_INTELTBB |
|
|
|
/* Official Intel Hint on blocked_range vs. linear iterators: http://software.intel.com/en-us/forums/topic/289505 |
|
|
|
|
|
|
|
@ -1755,19 +1708,18 @@ namespace gmm { |
|
|
|
{} |
|
|
|
}; |
|
|
|
#endif |
|
|
|
|
|
|
|
|
|
|
|
template <typename L1, typename L2, typename L3> |
|
|
|
void mult_by_row(const L1& l1, const L2& l2, L3& l3, abstract_dense) { |
|
|
|
typename linalg_traits<L3>::iterator it=vect_begin(l3), ite=vect_end(l3); |
|
|
|
typename linalg_traits<L1>::const_row_iterator |
|
|
|
itr = mat_row_const_begin(l1); |
|
|
|
auto itr = mat_row_const_begin(l1); |
|
|
|
#ifdef STORM_HAVE_INTELTBB |
|
|
|
tbb::parallel_for(forward_range_mult<typename linalg_traits<L3>::iterator, typename linalg_traits<L1>::const_row_iterator>(it, ite, itr), tbbHelper_mult_by_row<L1, L2, L3>(&l2)); |
|
|
|
#else |
|
|
|
for (; it != ite; ++it, ++itr) |
|
|
|
*it = vect_sp(linalg_traits<L1>::row(itr), l2, |
|
|
|
typename linalg_traits<L1>::storage_type(), |
|
|
|
typename linalg_traits<L2>::storage_type()); |
|
|
|
typename linalg_traits<L1>::storage_type(), |
|
|
|
typename linalg_traits<L2>::storage_type()); |
|
|
|
#endif |
|
|
|
} |
|
|
|
|
|
|
|
@ -1783,8 +1735,7 @@ namespace gmm { |
|
|
|
void mult_by_col(const L1& l1, const L2& l2, L3& l3, abstract_sparse) { |
|
|
|
typedef typename linalg_traits<L2>::value_type T; |
|
|
|
clear(l3); |
|
|
|
typename linalg_traits<L2>::const_iterator it = vect_const_begin(l2), |
|
|
|
ite = vect_const_end(l2); |
|
|
|
auto it = vect_const_begin(l2), ite = vect_const_end(l2); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (*it != T(0)) add(scaled(mat_const_col(l1, it.index()), *it), l3); |
|
|
|
} |
|
|
|
@ -1793,8 +1744,7 @@ namespace gmm { |
|
|
|
void mult_by_col(const L1& l1, const L2& l2, L3& l3, abstract_skyline) { |
|
|
|
typedef typename linalg_traits<L2>::value_type T; |
|
|
|
clear(l3); |
|
|
|
typename linalg_traits<L2>::const_iterator it = vect_const_begin(l2), |
|
|
|
ite = vect_const_end(l2); |
|
|
|
auto it = vect_const_begin(l2), ite = vect_const_end(l2); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (*it != T(0)) add(scaled(mat_const_col(l1, it.index()), *it), l3); |
|
|
|
} |
|
|
|
@ -1862,20 +1812,21 @@ namespace gmm { |
|
|
|
/** Multiply-accumulate. l4 = l1*l2 + l3; */ |
|
|
|
template <typename L1, typename L2, typename L3, typename L4> inline |
|
|
|
void mult_add(const L1& l1, const L2& l2, const L3& l3, L4& l4) { |
|
|
|
size_type m = mat_nrows(l1), n = mat_ncols(l1); |
|
|
|
if (!m || !n) return; |
|
|
|
GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l3) && vect_size(l3) == vect_size(l4), "dimensions mismatch"); |
|
|
|
if (!same_origin(l2, l3)) { |
|
|
|
mult_add_spec(l1, l2, l3, l4, typename principal_orientation_type<typename |
|
|
|
linalg_traits<L1>::sub_orientation>::potype()); |
|
|
|
} |
|
|
|
else { |
|
|
|
GMM_WARNING2("Warning, A temporary is used for mult\n"); |
|
|
|
typename temporary_vector<L3>::vector_type temp(vect_size(l2)); |
|
|
|
copy(l2, temp); |
|
|
|
mult_add_spec(l1, temp, l3, l4, typename principal_orientation_type<typename |
|
|
|
linalg_traits<L1>::sub_orientation>::potype()); |
|
|
|
} |
|
|
|
size_type m = mat_nrows(l1), n = mat_ncols(l1); |
|
|
|
if (!m || !n) return; |
|
|
|
GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l3) && vect_size(l3) == vect_size(l4), "dimensions mismatch"); |
|
|
|
if (!same_origin(l2, l4) && !same_origin(l3, l4) && !same_origin(l2, l3)) { |
|
|
|
mult_add_spec(l1, l2, l3, l4, typename principal_orientation_type<typename |
|
|
|
linalg_traits<L1>::sub_orientation>::potype()); |
|
|
|
} else { |
|
|
|
GMM_WARNING2("Warning, Multiple temporaries are used for mult\n"); |
|
|
|
typename temporary_vector<L2>::vector_type l2tmp(vect_size(l2)); |
|
|
|
copy(l2, l2tmp); |
|
|
|
typename temporary_vector<L3>::vector_type l3tmp(vect_size(l3)); |
|
|
|
copy(l3, l3tmp); |
|
|
|
mult_add_spec(l1, l2tmp, l3tmp, l4, typename principal_orientation_type<typename |
|
|
|
linalg_traits<L1>::sub_orientation>::potype()); |
|
|
|
} |
|
|
|
} |
|
|
|
///@cond DOXY_SHOW_ALL_FUNCTIONS |
|
|
|
|
|
|
|
@ -1905,21 +1856,20 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2, typename L3> |
|
|
|
void mult_add_by_row(const L1& l1, const L2& l2, L3& l3, abstract_dense) { |
|
|
|
typename linalg_traits<L3>::iterator it=vect_begin(l3), ite=vect_end(l3); |
|
|
|
typename linalg_traits<L1>::const_row_iterator |
|
|
|
itr = mat_row_const_begin(l1); |
|
|
|
auto it=vect_begin(l3), ite=vect_end(l3); |
|
|
|
auto itr = mat_row_const_begin(l1); |
|
|
|
for (; it != ite; ++it, ++itr) |
|
|
|
*it += vect_sp(linalg_traits<L1>::row(itr), l2); |
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
template <typename L1, typename L2, typename L3, typename L4> |
|
|
|
void mult_add_by_row(const L1& l1, const L2& l2, const L3& l3, L4& l4, abstract_dense) { |
|
|
|
typename linalg_traits<L3>::const_iterator add_it=vect_begin(l3), add_ite=vect_end(l3); |
|
|
|
typename linalg_traits<L4>::iterator target_it=vect_begin(l4), target_ite=vect_end(l4); |
|
|
|
typename linalg_traits<L1>::const_row_iterator |
|
|
|
itr = mat_row_const_begin(l1); |
|
|
|
typename linalg_traits<L3>::const_iterator add_it=vect_begin(l3), add_ite=vect_end(l3); |
|
|
|
typename linalg_traits<L4>::iterator target_it=vect_begin(l4), target_ite=vect_end(l4); |
|
|
|
typename linalg_traits<L1>::const_row_iterator |
|
|
|
itr = mat_row_const_begin(l1); |
|
|
|
for (; add_it != add_ite; ++add_it, ++target_it, ++itr) |
|
|
|
*target_it = vect_sp(linalg_traits<L1>::row(itr), l2) + *add_it; |
|
|
|
*target_it = vect_sp(linalg_traits<L1>::row(itr), l2) + *add_it; |
|
|
|
} |
|
|
|
|
|
|
|
template <typename L1, typename L2, typename L3> |
|
|
|
@ -1931,8 +1881,7 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2, typename L3> |
|
|
|
void mult_add_by_col(const L1& l1, const L2& l2, L3& l3, abstract_sparse) { |
|
|
|
typename linalg_traits<L2>::const_iterator it = vect_const_begin(l2), |
|
|
|
ite = vect_const_end(l2); |
|
|
|
auto it = vect_const_begin(l2), ite = vect_const_end(l2); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (*it != typename linalg_traits<L2>::value_type(0)) |
|
|
|
add(scaled(mat_const_col(l1, it.index()), *it), l3); |
|
|
|
@ -1940,8 +1889,7 @@ namespace gmm { |
|
|
|
|
|
|
|
template <typename L1, typename L2, typename L3> |
|
|
|
void mult_add_by_col(const L1& l1, const L2& l2, L3& l3, abstract_skyline) { |
|
|
|
typename linalg_traits<L2>::const_iterator it = vect_const_begin(l2), |
|
|
|
ite = vect_const_end(l2); |
|
|
|
auto it = vect_const_begin(l2), ite = vect_const_end(l2); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (*it != typename linalg_traits<L2>::value_type(0)) |
|
|
|
add(scaled(mat_const_col(l1, it.index()), *it), l3); |
|
|
|
@ -1954,7 +1902,7 @@ namespace gmm { |
|
|
|
template <typename L1, typename L2, typename L3, typename L4> inline |
|
|
|
void mult_add_spec(const L1& l1, const L2& l2, const L3& l3, L4& l4, row_major) |
|
|
|
{ mult_add_by_row(l1, l2, l3, l4, typename linalg_traits<L3>::storage_type()); } |
|
|
|
|
|
|
|
|
|
|
|
template <typename L1, typename L2, typename L3> inline |
|
|
|
void mult_add_spec(const L1& l1, const L2& l2, L3& l3, col_major) |
|
|
|
{ mult_add_by_col(l1, l2, l3, typename linalg_traits<L2>::storage_type()); } |
|
|
|
@ -2113,8 +2061,7 @@ namespace gmm { |
|
|
|
linalg_traits<L3>::sub_orientation>::potype()); |
|
|
|
} |
|
|
|
else { |
|
|
|
typename linalg_traits<L2>::const_col_iterator |
|
|
|
it2b = linalg_traits<L2>::col_begin(l2), it2, |
|
|
|
auto it2b = linalg_traits<L2>::col_begin(l2), it2 = it2b, |
|
|
|
ite = linalg_traits<L2>::col_end(l2); |
|
|
|
size_type i,j, k = mat_nrows(l1); |
|
|
|
|
|
|
|
@ -2152,8 +2099,7 @@ namespace gmm { |
|
|
|
size_type nn = mat_nrows(l3); |
|
|
|
for (size_type i = 0; i < nn; ++i) { |
|
|
|
typename linalg_traits<L1>::const_sub_row_type rl1=mat_const_row(l1, i); |
|
|
|
typename linalg_traits<typename linalg_traits<L1>::const_sub_row_type>:: |
|
|
|
const_iterator it = vect_const_begin(rl1), ite = vect_const_end(rl1); |
|
|
|
auto it = vect_const_begin(rl1), ite = vect_const_end(rl1); |
|
|
|
for (; it != ite; ++it) |
|
|
|
add(scaled(mat_const_row(l2, it.index()), *it), mat_row(l3, i)); |
|
|
|
} |
|
|
|
@ -2194,9 +2140,8 @@ namespace gmm { |
|
|
|
clear(l3); |
|
|
|
size_type nn = mat_ncols(l3); |
|
|
|
for (size_type i = 0; i < nn; ++i) { |
|
|
|
typename linalg_traits<L2>::const_sub_col_type rc2=mat_const_col(l2, i); |
|
|
|
typename linalg_traits<typename linalg_traits<L2>::const_sub_col_type>:: |
|
|
|
const_iterator it = vect_const_begin(rc2), ite = vect_const_end(rc2); |
|
|
|
typename linalg_traits<L2>::const_sub_col_type rc2 = mat_const_col(l2, i); |
|
|
|
auto it = vect_const_begin(rc2), ite = vect_const_end(rc2); |
|
|
|
for (; it != ite; ++it) |
|
|
|
add(scaled(mat_const_col(l1, it.index()), *it), mat_col(l3, i)); |
|
|
|
} |
|
|
|
@ -2246,9 +2191,8 @@ namespace gmm { |
|
|
|
clear(l3); |
|
|
|
size_type nn = mat_ncols(l1); |
|
|
|
for (size_type i = 0; i < nn; ++i) { |
|
|
|
typename linalg_traits<L1>::const_sub_col_type rc1=mat_const_col(l1, i); |
|
|
|
typename linalg_traits<typename linalg_traits<L1>::const_sub_col_type>:: |
|
|
|
const_iterator it = vect_const_begin(rc1), ite = vect_const_end(rc1); |
|
|
|
typename linalg_traits<L1>::const_sub_col_type rc1 = mat_const_col(l1, i); |
|
|
|
auto it = vect_const_begin(rc1), ite = vect_const_end(rc1); |
|
|
|
for (; it != ite; ++it) |
|
|
|
add(scaled(mat_const_row(l2, i), *it), mat_row(l3, it.index())); |
|
|
|
} |
|
|
|
@ -2299,10 +2243,8 @@ namespace gmm { |
|
|
|
bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol, |
|
|
|
row_major) { |
|
|
|
for (size_type i = 0; i < mat_nrows(A); ++i) { |
|
|
|
typedef typename linalg_traits<MAT>::const_sub_row_type row_type; |
|
|
|
row_type row = mat_const_row(A, i); |
|
|
|
typename linalg_traits<row_type>::const_iterator |
|
|
|
it = vect_const_begin(row), ite = vect_const_end(row); |
|
|
|
typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, i); |
|
|
|
auto it = vect_const_begin(row), ite = vect_const_end(row); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (gmm::abs(*it - A(it.index(), i)) > tol) return false; |
|
|
|
} |
|
|
|
@ -2313,10 +2255,8 @@ namespace gmm { |
|
|
|
bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol, |
|
|
|
col_major) { |
|
|
|
for (size_type i = 0; i < mat_ncols(A); ++i) { |
|
|
|
typedef typename linalg_traits<MAT>::const_sub_col_type col_type; |
|
|
|
col_type col = mat_const_col(A, i); |
|
|
|
typename linalg_traits<col_type>::const_iterator |
|
|
|
it = vect_const_begin(col), ite = vect_const_end(col); |
|
|
|
typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, i); |
|
|
|
auto it = vect_const_begin(col), ite = vect_const_end(col); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (gmm::abs(*it - A(i, it.index())) > tol) return false; |
|
|
|
} |
|
|
|
@ -2364,10 +2304,8 @@ namespace gmm { |
|
|
|
bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol, |
|
|
|
row_major) { |
|
|
|
for (size_type i = 0; i < mat_nrows(A); ++i) { |
|
|
|
typedef typename linalg_traits<MAT>::const_sub_row_type row_type; |
|
|
|
row_type row = mat_const_row(A, i); |
|
|
|
typename linalg_traits<row_type>::const_iterator |
|
|
|
it = vect_const_begin(row), ite = vect_const_end(row); |
|
|
|
typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, i); |
|
|
|
auto it = vect_const_begin(row), ite = vect_const_end(row); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (gmm::abs(gmm::conj(*it) - A(it.index(), i)) > tol) return false; |
|
|
|
} |
|
|
|
@ -2378,10 +2316,8 @@ namespace gmm { |
|
|
|
bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol, |
|
|
|
col_major) { |
|
|
|
for (size_type i = 0; i < mat_ncols(A); ++i) { |
|
|
|
typedef typename linalg_traits<MAT>::const_sub_col_type col_type; |
|
|
|
col_type col = mat_const_col(A, i); |
|
|
|
typename linalg_traits<col_type>::const_iterator |
|
|
|
it = vect_const_begin(col), ite = vect_const_end(col); |
|
|
|
typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, i); |
|
|
|
auto it = vect_const_begin(col), ite = vect_const_end(col); |
|
|
|
for (; it != ite; ++it) |
|
|
|
if (gmm::abs(gmm::conj(*it) - A(i, it.index())) > tol) return false; |
|
|
|
} |
dehnert
8 years ago
xxxxxxxxxx