problem with add( const VecX& x, VecY& y) in case x or y is sparse
Walter Daems <[email protected]> Tue, 20 Apr 2004 16:04:42 +0200
| Newsgroups | gmane.comp.lib.mtl.devel |
|---|---|
| Message-ID | <[email protected]> |
Hi all,
The problem I'm seeing concerns a function in mtl.h:
template <class VecX, class VecY> inline
void add__(const VecX& x, VecY& y, oned_tag);
does the job of adding two 1D vectors.
Actually, it gets called by
template <class LinalgA, class LinalgB> inline
void add(const LinalgA& A, MTL_OUT(LinalgB) B_);
In case both vectors are dense, everything goes well.
In case x or y is sparse (or both are sparse), bogus results
occur (due to the fact that mtl_algo::transform_add() is not
suited to deal with sparse vectors).
The patch attached solves the problem.
The behavior (and implementation) in the case of two dense
vectors is not changed.
bye,
Walter.
--
Walter Daems http://www.kimotion.com
mailto:[email protected]
Kimotion Technologies tel: +32 16 298306
Kapeldreef 60, B-3001 Leuven-Heverlee, Belgium fax: +32 16 298319
_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/
patch_add
(text/plain, 1.7 KB)
*** orig/mtl.h 2002-07-06 17:48:35.000000000 +0200
--- adpt/mtl.h 2004-04-20 15:28:17.000000000 +0200
***************
*** 2711,2723 ****
--- 2711,2762 ----
void
add__(const VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
+ typedef typename linalg_traits<VecX>::sparsity XSparsity;
+ typedef typename linalg_traits<VecY>::sparsity YSparsity;
+
+ add__(x, y, XSparsity(), YSparsity() );
+ }
+
+ template <class VecX, class VecY> inline
+ void
+ add__(const VecX& x, VecY& y, dense_tag, dense_tag)
+ {
add__(x, y, dim_n<VecX>::RET());
}
+ template <class VecX, class VecY> inline
+ void
+ add__(const VecX& x, VecY& y, sparse_tag, dense_tag)
+ {
+ typename VecX::const_iterator it = x.begin();
+ const typename VecX::const_iterator endit = x.end();
+ for ( ; not_at( it, endit); ++it )
+ y[ it.index() ] += *it;
+ }
+
+ template <class VecX, class VecY> inline
+ void
+ add__(const VecX& x, VecY& y, sparse_tag, sparse_tag)
+ {
+ typename VecX::const_iterator it = x.begin();
+ const typename VecX::const_iterator endit = x.end();
+ for ( ; not_at( it, endit); ++it )
+ y[ it.index() ] += *it;
+ }
+
+ template <class VecX, class VecY> inline
+ void
+ add__(const VecX& x, VecY& y, dense_tag, sparse_tag)
+ {
+ for ( typename VecX::size_type i = 0; i < x.size(); ++i )
+ if ( x[i] != 0.0 ) /* line may be omitted for speed reasons, but *
+ * in that case there's n point in using a sparse y vector */
+ y[i] += x[i];
+ }
template <class VecX, class VecY, class VecZ> inline
void
oned_add(const VecX& x, const VecY& y, VecZ& z, fast::count<0>)
{