@ -1,11 +1,11 @@
/* -*- c++ -*- (enables emacs c++ mode) */
/*===========================================================================
Copyright ( C ) 2002 - 2015 Yves Renard
This file is a part of GET FEM + +
Getfem + + is free software ; you can redistribute it and / or modify it
Copyright ( C ) 2002 - 2017 Yves Renard
This file is a part of Get FEM + +
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,11 +1643,9 @@ 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
*/
/ / Official Intel Hint on blocked_range vs . linear iterators : http : / / software . intel . com / en - us / forums / topic / 289505
template < typename IT1 , typename IT2 >
class forward_range_mult {
IT1 my_begin ;
@ -1754,23 +1705,119 @@ namespace gmm {
my_l2 ( l2 )
{ }
} ;
# endif
template < typename L1 , typename L2 , typename L3 , typename L4 >
class tbbHelper_mult_add_by_row {
L2 const * my_l2 ;
L3 const * my_l3 ;
/ / Typedefs for Iterator Types
typedef typename linalg_traits < L3 > : : iterator frm_IT1 ;
typedef typename linalg_traits < L1 > : : const_row_iterator frm_IT2 ;
typedef typename linalg_traits < L4 > : : const_iterator frm_IT3 ;
public :
void operator ( ) ( const forward_range_mult < frm_IT1 , frm_IT2 > & r ) const {
L2 const & l2 = * my_l2 ;
frm_IT1 it = r . begin ( ) ;
frm_IT1 ite = r . end ( ) ;
frm_IT2 itr = r . begin_row ( ) ;
frm_IT3 addIt = my_l3 - > begin ( ) ;
for ( ; it ! = ite ; + + it , + + itr , + + addIt ) {
* it = vect_sp ( linalg_traits < L1 > : : row ( itr ) , l2 ,
typename linalg_traits < L1 > : : storage_type ( ) ,
typename linalg_traits < L2 > : : storage_type ( ) ) + * addIt ;
}
}
tbbHelper_mult_add_by_row ( L2 const * l2 , L2 const * l3 ) : my_l2 ( l2 ) , my_l3 ( l3 ) {
/ / Intentionally left empty .
}
} ;
template < typename L1 , typename L2 , typename L3 >
class TbbMultFunctor {
public :
TbbMultFunctor ( L1 const & l1 , L2 const & l2 , L3 & l3 ) : l1 ( l1 ) , l2 ( l2 ) , l3 ( l3 ) {
/ / Intentionally left empty .
}
void operator ( ) ( tbb : : blocked_range < unsigned long > const & range ) const {
auto itr = mat_row_const_begin ( l1 ) + range . begin ( ) ;
auto l2it = l2 . begin ( ) + range . begin ( ) ;
auto l3it = l3 . begin ( ) + range . begin ( ) ;
auto l3ite = l3 . begin ( ) + range . end ( ) ;
for ( ; l3it ! = l3ite ; + + l3it , + + l2it , + + itr ) {
* l3it = vect_sp ( linalg_traits < L1 > : : row ( itr ) , l2 , typename linalg_traits < L1 > : : storage_type ( ) , typename linalg_traits < L2 > : : storage_type ( ) ) ;
}
}
private :
L1 const & l1 ;
L2 const & l2 ;
L3 & l3 ;
} ;
template < typename L1 , typename L2 , typename L3 >
void mult_by_row_parallel ( const L1 & l1 , const L2 & l2 , L3 & l3 , abstract_dense ) {
tbb : : parallel_for ( tbb : : blocked_range < unsigned long > ( 0 , vect_size ( l3 ) , 10 ) , TbbMultFunctor < L1 , L2 , L3 > ( l1 , l2 , l3 ) ) ;
}
template < typename L1 , typename L2 , typename L3 , typename L4 >
class TbbMultAddFunctor {
public :
TbbMultAddFunctor ( L1 const & l1 , L2 const & l2 , L3 const & l3 , L4 & l4 ) : l1 ( l1 ) , l2 ( l2 ) , l3 ( l3 ) , l4 ( l4 ) {
/ / Intentionally left empty .
}
void operator ( ) ( tbb : : blocked_range < unsigned long > const & range ) const {
auto itr = mat_row_const_begin ( l1 ) + range . begin ( ) ;
auto l2it = l2 . begin ( ) + range . begin ( ) ;
auto l3it = l3 . begin ( ) + range . begin ( ) ;
auto l4it = l4 . begin ( ) + range . begin ( ) ;
auto l4ite = l4 . begin ( ) + range . end ( ) ;
for ( ; l4it ! = l4ite ; + + l4it , + + l3it , + + l2it , + + itr ) {
* l4it = vect_sp ( linalg_traits < L1 > : : row ( itr ) , l2 , typename linalg_traits < L1 > : : storage_type ( ) , typename linalg_traits < L2 > : : storage_type ( ) ) + * l3it ;
}
}
private :
L1 const & l1 ;
L2 const & l2 ;
L3 const & l3 ;
L4 & l4 ;
} ;
template < typename L1 , typename L2 , typename L3 , typename L4 >
void mult_add_by_row_parallel ( const L1 & l1 , const L2 & l2 , const L3 & l3 , L4 & l4 , abstract_dense ) {
tbb : : parallel_for ( tbb : : blocked_range < unsigned long > ( 0 , vect_size ( l4 ) , 10 ) , TbbMultAddFunctor < L1 , L2 , L3 , L4 > ( l1 , l2 , l3 , l4 ) ) ;
}
# 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 ) ;
# 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 ( ) ) ;
# endif
auto itr = mat_row_const_begin ( l1 ) ;
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 ( ) ) ;
}
template < typename L1 , typename L2 , typename L3 >
void mult_by_row_bwd ( const L1 & l1 , const L2 & l2 , L3 & l3 , abstract_dense ) {
typename linalg_traits < L3 > : : iterator target_it = vect_end ( l3 ) - 1 , target_ite = vect_begin ( l3 ) - 1 ;
typename linalg_traits < L1 > : : const_row_iterator
itr = mat_row_const_end ( l1 ) - 1 ;
for ( ; target_it ! = target_ite ; - - target_it , - - itr )
* target_it = vect_sp ( linalg_traits < L1 > : : row ( itr ) , l2 ) ;
}
template < typename L1 , typename L2 , typename L3 >
void mult_by_col ( const L1 & l1 , const L2 & l2 , L3 & l3 , abstract_dense ) {
clear ( l3 ) ;
@ -1783,8 +1830,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 +1839,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 ) ;
}
@ -1803,6 +1848,12 @@ namespace gmm {
void mult_spec ( const L1 & l1 , const L2 & l2 , L3 & l3 , row_major )
{ mult_by_row ( l1 , l2 , l3 , typename linalg_traits < L3 > : : storage_type ( ) ) ; }
# ifdef STORM_HAVE_INTELTBB
template < typename L1 , typename L2 , typename L3 > inline
void mult_parallel_spec ( const L1 & l1 , const L2 & l2 , L3 & l3 , row_major )
{ mult_by_row_parallel ( l1 , l2 , l3 , typename linalg_traits < L3 > : : storage_type ( ) ) ; }
# endif
template < typename L1 , typename L2 , typename L3 > inline
void mult_spec ( const L1 & l1 , const L2 & l2 , L3 & l3 , col_major )
{ mult_by_col ( l1 , l2 , l3 , typename linalg_traits < L2 > : : storage_type ( ) ) ; }
@ -1862,21 +1913,60 @@ 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 ( ) ) ;
}
}
# ifdef STORM_HAVE_INTELTBB
/** Multiply. l3 = l1*l2; */
template < typename L1 , typename L2 , typename L3 > inline
void mult_parallel ( const L1 & l1 , const L2 & l2 , L3 & l3 ) {
size_type m = mat_nrows ( l1 ) , n = mat_ncols ( l1 ) ;
if ( ! m | | ! n ) return ;
GMM_ASSERT2 ( n = = vect_size ( l2 ) , " dimensions mismatch " ) ;
if ( ! same_origin ( l2 , l3 ) ) {
mult_parallel_spec ( l1 , l2 , l3 , typename principal_orientation_type < typename
linalg_traits < L1 > : : sub_orientation > : : potype ( ) ) ;
} else {
GMM_WARNING2 ( " Warning, temporaries are used for mult \n " ) ;
typename temporary_vector < L2 > : : vector_type l2tmp ( vect_size ( l2 ) ) ;
copy ( l2 , l2tmp ) ;
mult_parallel_spec ( l1 , l2tmp , l3 , typename principal_orientation_type < typename
linalg_traits < L1 > : : sub_orientation > : : potype ( ) ) ;
}
}
/** Multiply-accumulate. l4 = l1*l2 + l3; */
template < typename L1 , typename L2 , typename L3 , typename L4 > inline
void mult_add_parallel ( 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 , l4 ) & & ! same_origin ( l3 , l4 ) & & ! same_origin ( l2 , l3 ) ) {
mult_add_parallel_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_parallel_spec ( l1 , l2tmp , l3tmp , l4 , typename principal_orientation_type < typename linalg_traits < L1 > : : sub_orientation > : : potype ( ) ) ;
}
}
# endif
/ / / @ cond DOXY_SHOW_ALL_FUNCTIONS
template < typename L1 , typename L2 , typename L3 > inline
@ -1905,23 +1995,32 @@ 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 ) ;
for ( ; add_it ! = add_ite ; + + add_it , + + target_it , + + itr )
* target_it = vect_sp ( linalg_traits < L1 > : : row ( itr ) , l2 ) + * add_it ;
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 ;
}
template < typename L1 , typename L2 , typename L3 , typename L4 >
void mult_add_by_row_bwd ( const L1 & l1 , const L2 & l2 , const L3 & l3 , L4 & l4 , abstract_dense ) {
typename linalg_traits < L3 > : : const_iterator add_it = vect_end ( l3 ) - 1 , add_ite = vect_begin ( l3 ) - 1 ;
typename linalg_traits < L4 > : : iterator target_it = vect_end ( l4 ) - 1 , target_ite = vect_begin ( l4 ) - 1 ;
typename linalg_traits < L1 > : : const_row_iterator
itr = mat_row_const_end ( l1 ) - 1 ;
for ( ; add_it ! = add_ite ; - - add_it , - - target_it , - - itr )
* target_it = vect_sp ( linalg_traits < L1 > : : row ( itr ) , l2 ) + * add_it ;
}
template < typename L1 , typename L2 , typename L3 >
void mult_add_by_col ( const L1 & l1 , const L2 & l2 , L3 & l3 , abstract_dense ) {
size_type nc = mat_ncols ( l1 ) ;
@ -1931,8 +2030,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 +2038,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 +2051,13 @@ 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 ( ) ) ; }
# ifdef STORM_HAVE_INTELTBB
template < typename L1 , typename L2 , typename L3 , typename L4 > inline
void mult_add_parallel_spec ( const L1 & l1 , const L2 & l2 , const L3 & l3 , L4 & l4 , row_major )
{ mult_add_by_row_parallel ( l1 , l2 , l3 , l4 , typename linalg_traits < L4 > : : storage_type ( ) ) ; }
# endif
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 +2216,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 +2254,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 +2295,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 +2346,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 +2398,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 +2410,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 +2459,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 +2471,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 ;
}
cdehnert
8 years ago
GitHub
xxxxxxxxxx