Re: Linux/Intel C++: version 6.0.1 shows up MTL bugs
"Paul C. Leopardi" <[email protected]>
| Newsgroups | gmane.comp.lib.mtl.devel |
|---|---|
| Message-ID | <[email protected]> |
Hi all, With the help of Duraid Madina, I have discovered that there was a bug in the <map> file used by Intel C++ Version 6.0, which was fixed in Version 6.0.1. I now strongly suspect that the bugs I found in MTL 2.1.2-21 with respect to Intel C++ really are bugs in MTL and not the compiler, and the reason these bugs have not been found until now is that GluCat instantiates templates in MTL which are not instantiated by MTL's own test suite. I've attached an 18KB truncated error file and a 57KB patch to MTL which fixes these errors. Could someone please test this with icc and other compilers? Best regards
make-mtl-2.1.2-21-summary.out
(text/plain, 17.7 KB)
make all-recursive
make[1]: Entering directory `/c/leopardi/Download/glucat-0.0.7-icc'
Making all in glucat
make[2]: Entering directory `/c/leopardi/Download/glucat-0.0.7-icc/glucat'
make[3]: Entering directory `/c/leopardi/Download/glucat-0.0.7-icc/glucat'
make[3]: Nothing to be done for `all-am'.
make[3]: Leaving directory `/c/leopardi/Download/glucat-0.0.7-icc/glucat'
make[2]: Leaving directory `/c/leopardi/Download/glucat-0.0.7-icc/glucat'
Making all in test
make[2]: Entering directory `/c/leopardi/Download/glucat-0.0.7-icc/test'
source='driver.cpp' object='driver.o' libtool=no \
depfile='.deps/driver.Po' tmpdepfile='.deps/driver.TPo' \
depmode=none /bin/sh ../admin/depcomp \
icc -D__ICL -DHAVE_CONFIG_H -I. -I. -I.. -D__ICL -mp -Xc -O0 -c -o driver.o `test -f 'driver.cpp' || echo './'`driver.cpp
driver.cpp
./mtl/mtl_algo.h(217): error: invalid destructor name for type "std::iterator_traits<_Iter>::value_type"
sum1.~sum_type();
^
./mtl/mtl_algo.h(240): error: invalid destructor name for type "std::iterator_traits<_Iter>::value_type"
sum1.~sum_type();
^
./mtl/reverse_iter.h(71): error: identifier "current" is undefined
Iter tmp = current;
^
./mtl/reverse_iter.h(75): error: identifier "current" is undefined
Iter tmp = current;
^
./mtl/reverse_iter.h(79): error: identifier "current" is undefined
Iter tmp = current;
^
./mtl/transform_iterator.h(267): error: identifier "iter" is undefined
stride_gen.inc(iter);
^
./mtl/transform_iterator.h(271): error: identifier "iter" is undefined
stride_gen.dec(iter);
^
./mtl/transform_iterator.h(275): error: identifier "iter" is undefined
stride_gen.advance(iter, n);
^
./mtl/transform_iterator.h(279): error: identifier "iter" is undefined
stride_gen.advance(iter, -n);
^
./mtl/transform_iterator.h(284): error: identifier "iter" is undefined
return stride_gen.diff(iter, y.iter, y.stride_gen);
^
./mtl/linalg_vec.h(149): error: identifier "Vec_value_type" is undefined
while (*i == Vec_value_type(0)) ++i;
^
./mtl/linalg_vec.h(153): error: identifier "Vec_value_type" is undefined
++i; while (*i == Vec_value_type(0) && i != end) ++i;
^
./mtl/linalg_vec.h(157): error: identifier "Vec_value_type" is undefined
--i; while (*i == Vec_value_type(0) && i != end) --i;
^
./mtl/linalg_vec.h(176): error: identifier "Vec_value_type" is undefined
while (*i == Vec_value_type(0) && i != end) ++i;
^
./mtl/linalg_vec.h(180): error: identifier "Vec_value_type" is undefined
++i; while (*i == Vec_value_type(0) && i != end) ++i;
^
./mtl/linalg_vec.h(185): error: identifier "Vec_value_type" is undefined
--i; while (*i == Vec_value_type(0)) --i;
^
./mtl/linalg_vec.h(213): error: identifier "Vec_value_type" is undefined
if (*i != Vec_value_type(0)) ++s;
^
./mtl/matrix_implementation.h(351): error: namespace "mtl" has no member "initialize"
mtl::initialize(me, m_in);
^
./mtl/matrix_implementation.h(358): error: namespace "mtl" has no member "initialize"
mtl::initialize(me, m_in);
^
./mtl/matrix_implementation.h(369): error: namespace "mtl" has no member "initialize"
mtl::initialize(me, m_in);
^
./mtl/matrix_implementation.h(380): error: namespace "mtl" has no member "initialize"
mtl::initialize(me, m_in);
^
./mtl/matrix_implementation.h(539): error: identifier "band_type" is undefined
: Base(dim_type(m, n), band_type(sub, super)) { }
^
./mtl/matrix_implementation.h(555): error: identifier "band_type" is undefined
band_type(sub, super)) { }
^
./mtl/matrix_implementation.h(559): error: identifier "band_type" is undefined
band_type(sub, super)) { }
^
./mtl/matrix_implementation.h(572): error: identifier "band_type" is undefined
: Base(x, band_type(sub, super)) { }
^
./mtl/matrix_implementation.h(585): error: identifier "band_type" is undefined
: Base(m_in, band_type(sub, super), *this) { }
^
./mtl/matrix_implementation.h(587): error: identifier "band_type" is undefined
: Base(m_in, band_type(sub, super), *this) { }
^
./mtl/matrix_implementation.h(595): error: identifier "band_type" is undefined
: Base(m_in, band_type(sub, super), s) { }
^
./mtl/matrix_implementation.h(598): error: identifier "band_type" is undefined
: Base(m_in, band_type(sub, super), s) { }
^
./mtl/matrix_implementation.h(627): error: identifier "twod" is undefined
twod.resize(m, n);
^
./mtl/matrix_implementation.h(628): error: identifier "indexer" is undefined
indexer.dim = dim_type(m, n);
^
./mtl/matrix_implementation.h(657): error: identifier "indexer" is undefined
dim_type starts = indexer.at(dim_type(row_start, col_start));
^
./mtl/matrix_implementation.h(662): error: identifier "twod" is undefined
return submatrix_type((value_type*)twod.data()
^
./mtl/matrix_implementation.h(666): error: identifier "twod" is undefined
return submatrix_type((value_type*)twod.data()
^
./mtl/matrix_implementation.h(854): error: identifier "twod" is undefined
twod.resize(n, m);
^
./mtl/matrix_implementation.h(855): error: identifier "indexer" is undefined
indexer.dim = dim_type(n, m);
^
./mtl/matrix_implementation.h(880): error: identifier "indexer" is undefined
dim_type starts = indexer.at(dim_type(row_start, col_start));
^
./mtl/matrix_implementation.h(885): error: identifier "twod" is undefined
return submatrix_type((value_type*)twod.data()
^
./mtl/matrix_implementation.h(889): error: identifier "twod" is undefined
return submatrix_type((value_type*)twod.data()
^
./mtl/matrix_implementation.h(1231): error: identifier "m" is undefined
: Base(m, n, Uplo::bandwidth(uplo_, sub, sub)) , uplo(uplo_) { }
^
./mtl/matrix_implementation.h(1235): error: identifier "m" is undefined
: Base(d, m, n,
^
./mtl/matrix_implementation.h(1311): error: identifier "indexer" is undefined
return MTL_MAX(indexer.super(), indexer.sub());
^
./mtl/matrix_implementation.h(1314): error: identifier "indexer" is undefined
return MTL_MAX(indexer.super(), indexer.sub());
^
./mtl/initialize.h(61): error: identifier "cout" is undefined
cout << " matrix type is symmetric but the matrix in the file is not" << endl;
^
./mtl/initialize.h(61): error: identifier "endl" is undefined
cout << " matrix type is symmetric but the matrix in the file is not" << endl;
^
./mtl/block2D.h(318): error: identifier "s" is undefined
return block_vector(s, offset.oned_length(i * block_dim.first()),
^
./mtl/dense2D.h(1238): error: identifier "offset" is undefined
size_type M = MTL_MIN(m, offset.major());
^
./mtl/dense2D.h(1247): error: identifier "data_" is undefined
data_ = newdata;
^
./mtl/dense2D.h(1248): error: identifier "ld_" is undefined
ld_ = n;
^
./mtl/dense2D.h(1366): error: identifier "data_" is undefined
data_ = &rep;
^
./mtl/utils.h(53): error: namespace "mtl" has no member "print_vector"
mtl::print_vector(*i);
^
./mtl/utils.h(66): error: namespace "mtl" has no member "print_all_matrix"
mtl::print_all_matrix( A(i,j) );
^
./mtl/utils.h(83): error: namespace "mtl" has no member "print_all_matrix"
mtl::print_all_matrix( *A_i );
^
./mtl/utils.h(102): error: namespace "mtl" has no member "print_all_matrix"
mtl::print_all_matrix( *A_i );
^
./mtl/compressed2D.h(615): error: identifier "cout" is undefined
cout << "values ";
^
./mtl/compressed2D.h(669): error: identifier "dim" is undefined
vals.reserve(dim.first() * 5);
^
./mtl/compressed2D.h(686): error: identifier "dim" is undefined
vals.reserve(dim.first() * 5);
^
./mtl/compressed2D.h(705): error: identifier "dim" is undefined
dim = x.dim;
^
./mtl/compressed2D.h(706): error: identifier "values" is undefined
values = &vals; indices = &inds; starts = &ptrs;
^
./mtl/compressed2D.h(706): error: identifier "indices" is undefined
values = &vals; indices = &inds; starts = &ptrs;
^
./mtl/compressed2D.h(706): error: identifier "starts" is undefined
values = &vals; indices = &inds; starts = &ptrs;
^
./mtl/compressed2D.h(804): error: identifier "dim" is undefined
dim = x.dim;
^
./mtl/compressed2D.h(805): error: identifier "values" is undefined
values = &vals; indices = &inds; starts = &ptrs;
^
./mtl/compressed2D.h(805): error: identifier "indices" is undefined
values = &vals; indices = &inds; starts = &ptrs;
^
./mtl/compressed2D.h(805): error: identifier "starts" is undefined
values = &vals; indices = &inds; starts = &ptrs;
^
./mtl/envelope2D.h(299): error: identifier "m" is undefined
: dim(d), val(1, T(0)), ptr(m, size_type(0))
^
./mtl/mtl.h(1542): error: namespace "mtl" has no member "add"
mtl::add(mtl::scaled(A[xi.index()], *xi), y);
^
./mtl/mtl.h(1555): error: namespace "mtl" has no member "dot"
y[Ai.index()] = mtl::dot(*Ai, x); // this is a sparse dot
^
./mtl/mtl.h(1583): error: namespace "mtl" has no member "add"
mtl::add(z, x, z);
^
./mtl/mtl.h(1585): error: namespace "mtl" has no member "add"
mtl::add(x, z, z);
^
./mtl/mtl.h(1675): error: namespace "mtl" has no member "copy"
mtl::copy(y, z);
^
./mtl/mtl.h(2627): error: no instance of function template "mtl::row_matrix<TwoDGen, IndexerGen>::fast_copy [with TwoDGen=mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::storage_t, IndexerGen=mtl::gen_rect_indexer<mtl::row_orien, 0, 0, mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::size_type={mtl::compressed2D<double, mtl::compressed<int, 0, 0>::size_type={int}, 0>::size_type={mtl::compressed<int, 0, 0>::size_type={int}}}>]" matches
the argument list
argument types are: (const mtl::row_matrix<mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::storage_t, mtl::gen_rect_indexer<mtl::row_orien, 0, 0, mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::size_type>>)
B.fast_copy(A);
^
detected during:
instantiation of "void mtl::twod_copy(const MatA &, MatB &, mtl::row_tag, mtl::row_tag) [with MatA=mtl::row_matrix<mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::storage_t, mtl::gen_rect_indexer<mtl::row_orien, 0, 0, mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::size_type={mtl::compressed2D<double, mtl::compressed<int, 0, 0>::size_type={int}, 0>::size_type={mtl::compressed<int, 0, 0>::size_type={int}}}>>,
MatB=glucat::matrix_multi<double, (glucat::index_t={short})(-15), (glucat::index_t={short})15>::matrix_t]" at line 2655
instantiation of "void mtl::twod_copy__(const MatA &, MatB &, mtl::sparse_tag) [with MatA=mtl::row_matrix<mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::storage_t, mtl::gen_rect_indexer<mtl::row_orien, 0, 0, mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::size_type={mtl::compressed2D<double, mtl::compressed<int, 0, 0>::size_type={int}, 0>::size_type={mtl::compressed<int, 0, 0>::size_type={int}}}>>,
MatB=glucat::matrix_multi<double, (glucat::index_t={short})(-15), (glucat::index_t={short})15>::matrix_t]" at line 2666
instantiation of "void mtl::copy__(const MatA &, MatB &, mtl::twod_tag) [with MatA=mtl::row_matrix<mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::storage_t, mtl::gen_rect_indexer<mtl::row_orien, 0, 0, mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::size_type={mtl::compressed2D<double, mtl::compressed<int, 0, 0>::size_type={int}, 0>::size_type={mtl::compressed<int, 0, 0>::size_type={int}}}>>,
MatB=glucat::matrix_multi<double, (glucat::index_t={short})(-15), (glucat::index_t={short})15>::matrix_t]" at line 2689
instantiation of "void mtl::copy(const LinalgA &, const LinalgB &) [with LinalgA=mtl::row_matrix<mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::storage_t, mtl::gen_rect_indexer<mtl::row_orien, 0, 0, mtl::generate_rect<double, mtl::row_major, mtl::compressed<int, 0, 0>, 0, 0>::size_type={mtl::compressed2D<double, mtl::compressed<int, 0, 0>::size_type={int}, 0>::size_type={mtl::compressed<int, 0, 0>::size_type={int}}}>>,
LinalgB=glucat::matrix_multi<double, (glucat::index_t={short})(-15), (glucat::index_t={short})15>::matrix_t]" at line 199 of "../glucat/generation_imp.h"
instantiation of "void glucat::generator_table<Matrix_T>::gen_from_pp4_qm4(const std::vector<Scalar_T, std::allocator<_Elem>> &, glucat::signature_t) [with Matrix_T=glucat::matrix_multi<double, (glucat::index_t={short})(-15), (glucat::index_t={short})15>::matrix_t]" at line 88 of "../glucat/generation_imp.h"
instantiation of "const std::vector<Scalar_T, std::allocator<_Elem>> &glucat::generator_table<Matrix_T>::gen_vector(glucat::index_t={short}, glucat::index_t={short}) [with Matrix_T=glucat::matrix_multi<double, (glucat::index_t={short})(-15), (glucat::index_t={short})15>::matrix_t]" at line 61 of "../glucat/generation_imp.h"
instantiation of "const Matrix_T *glucat::generator_table<Matrix_T>::operator()(glucat::index_t={short}, glucat::index_t={short}) [with Matrix_T=glucat::matrix_multi<double, (glucat::index_t={short})(-15), (glucat::index_t={short})15>::matrix_t]" at line 753 of "../glucat/matrix_multi_imp.h"
instantiation of "void glucat::basis_element(const glucat::index_set<LO, HI> &, const glucat::index_set<LO, HI> &, glucat::matrix_multi<Scalar_T, LO, HI>::matrix_t &) [with Scalar_T=double, LO=(short)(-15), HI=(short)15]" at line 142 of "../glucat/framed_multi_imp.h"
instantiation of "glucat::framed_multi<Scalar_T, LO, HI>::framed_multi(const glucat::framed_multi<Scalar_T, LO, HI>::matrix_multi_t &) [with Scalar_T=double, LO=(glucat::index_t={short})(-15), HI=(glucat::index_t={short})15]" at line 412 of "../glucat/matrix_multi_imp.h"
instantiation of "glucat::matrix_multi<Scalar_T, LO, HI>::multivector_t &glucat::matrix_multi<Scalar_T, LO, HI>::operator/=(const glucat::matrix_multi<Scalar_T, LO, HI>::multivector_t &) [with Scalar_T=double, LO=(glucat::index_t={short})(-15), HI=(glucat::index_t={short})15]" at line 358 of "../glucat/framed_multi_imp.h"
instantiation of "glucat::framed_multi<Scalar_T, LO, HI>::multivector_t &glucat::framed_multi<Scalar_T, LO, HI>::operator/=(const glucat::framed_multi<Scalar_T, LO, HI>::multivector_t &) [with Scalar_T=double, LO=(glucat::index_t={short})(-15), HI=(glucat::index_t={short})15]" at line 40 of "../glucat/framed_multi_imp.h"
instantiation of "glucat::framed_multi<Scalar_T, LO, HI>::framed_multi() [with Scalar_T=double, LO=(glucat::index_t={short})(-15), HI=(glucat::index_t={short})15]" at line 246 of "../glucat/framed_multi_imp.h"
instantiation of "glucat::framed_multi<Scalar_T, LO, HI>::multivector_t &glucat::framed_multi<Scalar_T, LO, HI>::operator*=(const glucat::framed_multi<Scalar_T, LO, HI>::multivector_t &) [with Scalar_T=double, LO=(glucat::index_t={short})(-15), HI=(glucat::index_t={short})15]" at line 176 of "../glucat/framed_multi_imp.h"
instantiation of "bool glucat::framed_multi<Scalar_T, LO, HI>::operator==(const Scalar_T &) const [with Scalar_T=double, LO=(glucat::index_t={short})(-15), HI=(glucat::index_t={short})15]" at line 226 of "../glucat/clifford_algebra_imp.h"
implicit generation of "glucat::framed_multi<Scalar_T, LO, HI>::framed_multi(const glucat::framed_multi<double, (glucat::index_t={short})(-15), (glucat::index_t={short})15> &) [with Scalar_T=double, LO=(glucat::index_t={short})(-15), HI=(glucat::index_t={short})15]" at line 226 of "../glucat/clifford_algebra_imp.h"
instantiation of "const Multivector<Scalar_T, LO, HI> glucat::operator^<Multivector,RHS,Scalar_T,LO,HI>(const Multivector<Scalar_T, LO, HI> &, const RHS<Scalar_T, LO, HI> &) [with Multivector=glucat::framed_multi, RHS=glucat::framed_multi, Scalar_T=double, LO=(glucat::index_t={short})(-15), HI=(glucat::index_t={short})15]" at line 43 of "../test/peg01.h"
instantiation of "void peg01::do_test1<Multivector_T>() [with Multivector_T=glucat::framed_multi<double, (glucat::index_t={short})(-15), (glucat::index_t={short})15>]"
[...]
compilation aborted for driver.cpp (code 2)
make[2]: *** [driver.o] Error 2
make[2]: Leaving directory `/c/leopardi/Download/glucat-0.0.7-icc/test'
make[1]: *** [all-recursive] Error 1
make[1]: Leaving directory `/c/leopardi/Download/glucat-0.0.7-icc'
make: *** [all] Error 2
mtl-2.1.2-21-icc.diff
(text/x-diff, 56.5 KB)
diff -drub mtl-2.1.2-21/mtl/block2D.h mtl-2.1.2-21-icc/mtl/block2D.h
--- mtl-2.1.2-21/mtl/block2D.h Sun Jul 7 01:48:34 2002
+++ mtl-2.1.2-21-icc/mtl/block2D.h Fri Aug 23 20:22:29 2002
@@ -315,7 +315,7 @@
//: Return a reference to the ith 1D container
inline block_vector operator[](size_type i) {
- return block_vector(s, offset.oned_length(i * block_dim.first()),
+ return block_vector(data_, offset.oned_length(i * block_dim.first()),
block_dim, offset.stride(), ld_);
}
diff -drub mtl-2.1.2-21/mtl/compressed2D.h mtl-2.1.2-21-icc/mtl/compressed2D.h
--- mtl-2.1.2-21/mtl/compressed2D.h Sun Jul 7 01:48:34 2002
+++ mtl-2.1.2-21-icc/mtl/compressed2D.h Fri Aug 23 20:32:20 2002
@@ -612,11 +612,11 @@
}
void print() const {
- cout << "values ";
+ std::cout << "values ";
print_vector(*values);
- cout << "indices ";
+ std::cout << "indices ";
print_vector(*indices);
- cout << "starts ";
+ std::cout << "starts ";
print_vector(*starts);
}
@@ -666,8 +666,8 @@
: Base(d, &vals, &inds, &ptrs),
ptrs(d.first() + 1, -IND_OFFSET) /* F to C */
{
- vals.reserve(dim.first() * 5);
- inds.reserve(dim.first() * 5);
+ vals.reserve(this->dim.first() * 5);
+ inds.reserve(this->dim.first() * 5);
}
inline compressed2D(dim_type d, size_type nnz)
@@ -683,8 +683,8 @@
: Base(d, &vals, &inds, &ptrs),
ptrs(d.first() + 1, -IND_OFFSET) /* F to C */
{
- vals.reserve(dim.first() * 5);
- inds.reserve(dim.first() * 5);
+ vals.reserve(this->dim.first() * 5);
+ inds.reserve(this->dim.first() * 5);
}
//: Copy Constructor
@@ -702,8 +702,8 @@
inline self& operator=(const self& x) {
vals = x.vals; ptrs = x.ptrs; inds = x.inds;
/* fill out inhereted part */
- dim = x.dim;
- values = &vals; indices = &inds; starts = &ptrs;
+ this->dim = x.dim;
+ this->values = &vals; this->indices = &inds; this->starts = &ptrs;
return *this;
}
@@ -801,8 +801,8 @@
inline self& operator=(const self& x) {
vals = x.vals; inds = x.inds; ptrs = x.ptrs;
/* fill out inhereted part */
- dim = x.dim;
- values = &vals; indices = &inds; starts = &ptrs;
+ this->dim = x.dim;
+ this->values = &vals; this->indices = &inds; this->starts = &ptrs;
return *this;
}
diff -drub mtl-2.1.2-21/mtl/dense2D.h mtl-2.1.2-21-icc/mtl/dense2D.h
--- mtl-2.1.2-21/mtl/dense2D.h Sun Jul 7 01:48:34 2002
+++ mtl-2.1.2-21-icc/mtl/dense2D.h Fri Aug 23 20:29:47 2002
@@ -1235,8 +1235,8 @@
inline void resize(size_type m, size_type n) {
rep_ptr newdata = new reptype(Offset::size(m, n, 0, 0));
size_type i, j;
- size_type M = MTL_MIN(m, offset.major());
- size_type N = MTL_MIN(n, offset.minor());
+ size_type M = MTL_MIN(m, this->offset.major());
+ size_type N = MTL_MIN(n, this->offset.minor());
for (i = 0; i < M; ++i)
for (j = 0; j < N; ++j)
(*newdata)[i * n + j] = (*this)(i,j);
@@ -1244,10 +1244,10 @@
for (; j < n; ++j)
(*newdata)[i * n + j] = T();
- data_ = newdata;
- ld_ = n;
- offset.dim = dim_type(m, n);
- offset.ld = n;
+ this->data_ = newdata;
+ this->ld_ = n;
+ this->offset.dim = dim_type(m, n);
+ this->offset.ld = n;
}
};
@@ -1363,7 +1363,7 @@
inline external2D& operator=(const external2D& x) {
rep = x.rep;
super::operator=(x);
- data_ = &rep;
+ this->data_ = &rep;
return *this;
}
diff -drub mtl-2.1.2-21/mtl/envelope2D.h mtl-2.1.2-21-icc/mtl/envelope2D.h
--- mtl-2.1.2-21/mtl/envelope2D.h Sun Jul 7 01:48:34 2002
+++ mtl-2.1.2-21-icc/mtl/envelope2D.h Fri Aug 23 20:34:30 2002
@@ -296,9 +296,9 @@
inline envelope2D() : dim(0,0) { }
//: Constructor from dimension pair
inline envelope2D(dim_type d)
- : dim(d), val(1, T(0)), ptr(m, size_type(0))
+ : dim(d), val(1, T(0)), ptr(d.first(), size_type(0))
{
- val.reserve(m * 5);
+ val.reserve(dim.first() * 5);
}
//: Constructor from dimension abd bandwidth pairs
inline envelope2D(dim_type d, dim_type)
Only in mtl-2.1.2-21-icc/mtl: grep.MTL_CONJ.txt
Only in mtl-2.1.2-21-icc/mtl: grep.array2D.txt
Only in mtl-2.1.2-21-icc/mtl: grep.conj.txt
Only in mtl-2.1.2-21-icc/mtl: grep.fast_copy.txt
Only in mtl-2.1.2-21-icc/mtl: grep.initialize.txt
Only in mtl-2.1.2-21-icc/mtl: grep.print_.txt
diff -drub mtl-2.1.2-21/mtl/initialize.h mtl-2.1.2-21-icc/mtl/initialize.h
--- mtl-2.1.2-21/mtl/initialize.h Sun Jul 7 01:48:34 2002
+++ mtl-2.1.2-21-icc/mtl/initialize.h Fri Aug 23 20:18:47 2002
@@ -58,7 +58,7 @@
A(row, col) = e.value;
}
} else {
- cout << " matrix type is symmetric but the matrix in the file is not" << endl;
+ std::cout << " matrix type is symmetric but the matrix in the file is not" << std::endl;
assert(0);
}
}
diff -drub mtl-2.1.2-21/mtl/linalg_vec.h mtl-2.1.2-21-icc/mtl/linalg_vec.h
--- mtl-2.1.2-21/mtl/linalg_vec.h Sun Jul 7 01:48:35 2002
+++ mtl-2.1.2-21-icc/mtl/linalg_vec.h Fri Aug 23 20:01:50 2002
@@ -130,6 +130,7 @@
typedef difference_type Vec_difference_type;
typedef iterator Vec_iterator;
typedef const_iterator Vec_const_iterator;
+ typedef value_type Vec_value_type;
class IndexArray {
public:
diff -drub mtl-2.1.2-21/mtl/matrix.h mtl-2.1.2-21-icc/mtl/matrix.h
--- mtl-2.1.2-21/mtl/matrix.h Sun Jul 7 01:48:35 2002
+++ mtl-2.1.2-21-icc/mtl/matrix.h Sat Aug 24 16:12:33 2002
@@ -44,7 +44,6 @@
#include "mtl/matrix_traits.h"
#include "mtl/compressed1D.h"
-#include "mtl/matrix_implementation.h"
#include "mtl/rect_indexer.h"
#include "mtl/banded_indexer.h"
#include "mtl/diagonal_indexer.h"
@@ -58,6 +57,7 @@
#include "mtl/sparse1D.h"
#include "mtl/entry.h"
#include "mtl/uplo.h"
+#include "mtl/matrix_implementation.h"
namespace mtl {
diff -drub mtl-2.1.2-21/mtl/matrix_implementation.h mtl-2.1.2-21-icc/mtl/matrix_implementation.h
--- mtl-2.1.2-21/mtl/matrix_implementation.h Sun Jul 7 01:48:35 2002
+++ mtl-2.1.2-21-icc/mtl/matrix_implementation.h Sat Aug 24 20:04:21 2002
@@ -492,6 +492,7 @@
typedef typename Base::OneD OneD;
typedef typename Base::OneDRef OneDRef;
typedef typename Base::Indexer Indexer;
+ typedef typename Base::band_type band_type;
typedef OneD Row;
typedef OneDRef RowRef;
enum { M = Indexer::M, N = Indexer::N };
@@ -524,6 +525,9 @@
};
#endif
+ template <class Matrix>
+ inline void fast_copy(const Matrix& x) { Base::fast_copy(x); }
+
inline row_matrix() { }
inline row_matrix& operator=(const row_matrix& x) {
@@ -624,8 +628,8 @@
inline ~row_matrix() { }
inline void resize(size_type m, size_type n) {
- twod.resize(m, n);
- indexer.dim = dim_type(m, n);
+ this->twod.resize(m, n);
+ this->indexer.dim = dim_type(m, n);
}
/* submatrix */
@@ -654,18 +658,18 @@
inline submatrix_type sub_matrix(size_type row_start, size_type row_finish,
size_type col_start, size_type col_finish) const
{
- dim_type starts = indexer.at(dim_type(row_start, col_start));
+ dim_type starts = this->indexer.at(dim_type(row_start, col_start));
size_type m = row_finish - row_start;
size_type n = col_finish - col_start;
typedef typename TwoD::is_strided IsStrided;
if (IsStrided::id) {
- return submatrix_type((value_type*)twod.data()
- + starts.second() * twod.ld() + starts.first(),
- m, n, twod.ld());
+ return submatrix_type((value_type*)this->twod.data()
+ + starts.second() * this->twod.ld() + starts.first(),
+ m, n, this->twod.ld());
} else {
- return submatrix_type((value_type*)twod.data()
- + starts.first() * twod.ld() + starts.second(),
- m, n, twod.ld());
+ return submatrix_type((value_type*)this->twod.data()
+ + starts.first() * this->twod.ld() + starts.second(),
+ m, n, this->twod.ld());
}
}
@@ -759,6 +763,9 @@
: Base(dim_type(m, n),
band_type(sub, super)) { }
+ template <class Matrix>
+ inline void fast_copy(const Matrix& x) { Base::fast_copy(x); }
+
inline column_matrix& operator=(const column_matrix& x) {
Base::operator=(x);
return *this;
@@ -851,8 +858,8 @@
inline ~column_matrix() { }
inline void resize(size_type m, size_type n) {
- twod.resize(n, m);
- indexer.dim = dim_type(n, m);
+ this->twod.resize(n, m);
+ this->indexer.dim = dim_type(n, m);
}
/* submatrix */
@@ -877,18 +884,18 @@
inline submatrix_type sub_matrix(size_type row_start, size_type row_finish,
size_type col_start, size_type col_finish) const
{
- dim_type starts = indexer.at(dim_type(row_start, col_start));
+ dim_type starts = this->indexer.at(dim_type(row_start, col_start));
size_type m = row_finish - row_start;
size_type n = col_finish - col_start;
typedef typename TwoD::is_strided IsStrided;
if (IsStrided::id) {
- return submatrix_type((value_type*)twod.data()
- + starts.second() * twod.ld() + starts.first(),
- m, n, twod.ld());
+ return submatrix_type((value_type*)this->twod.data()
+ + starts.second() * this->twod.ld() + starts.first(),
+ m, n, this->twod.ld());
} else {
- return submatrix_type((value_type*)twod.data()
- + starts.first() * twod.ld() + starts.second(),
- m, n, twod.ld());
+ return submatrix_type((value_type*)this->twod.data()
+ + starts.first() * this->twod.ld() + starts.second(),
+ m, n, this->twod.ld());
}
}
@@ -1228,11 +1235,11 @@
//: dynamic uplo constructor
inline symmetric_matrix(size_type n, int uplo_, int sub)
- : Base(m, n, Uplo::bandwidth(uplo_, sub, sub)) , uplo(uplo_) { }
+ : Base(n, n, Uplo::bandwidth(uplo_, sub, sub)) , uplo(uplo_) { }
//: constructor for external data with dynamic uplo
inline symmetric_matrix(pointer d, size_type n, int uplo_, int sub)
- : Base(d, m, n,
+ : Base(d, n, n,
Uplo::bandwidth(uplo_, sub, sub).first,
Uplo::bandwidth(uplo_, sub, sub).second) { }
@@ -1308,10 +1315,10 @@
/* bandwidth (is symmetric too) */
inline int sub() const {
- return MTL_MAX(indexer.super(), indexer.sub());
+ return MTL_MAX(this->indexer.super(), this->indexer.sub());
}
inline int super() const {
- return MTL_MAX(indexer.super(), indexer.sub());
+ return MTL_MAX(this->indexer.super(), this->indexer.sub());
}
};
diff -drub mtl-2.1.2-21/mtl/mtl.h mtl-2.1.2-21-icc/mtl/mtl.h
--- mtl-2.1.2-21/mtl/mtl.h Sun Jul 7 01:48:35 2002
+++ mtl-2.1.2-21-icc/mtl/mtl.h Sat Aug 24 16:23:56 2002
@@ -1434,6 +1434,739 @@
}
+template <class VecX, class VecY>
+inline void
+copy__(const VecX& x, VecY& y, fast::count<0>)
+{
+ mtl_algo::copy(x.begin(), x.end(), y.begin());
+}
+#if USE_BLAIS
+template <class VecX, class VecY, int N>
+inline void
+copy__(const VecX& x, VecY& y, fast::count<N>)
+{
+ fast::copy(x.begin(), fast::count<N>(), y.begin());
+}
+#endif
+
+
+template <class VecX, class VecY>
+inline void
+oned_copy(const VecX& x, VecY& y, dense_tag, dense_tag) MTL_THROW_ASSERTION
+{
+ MTL_ASSERT(x.size() <= y.size(), "mtl::copy()");
+ copy__(x, y, dim_n<VecX>::RET());
+}
+
+#if 0
+/* perform a scatter */
+template <class VecX, class VecY>
+inline void
+oned_copy(const VecX& x, VecY y, sparse_tag, dense_tag) MTL_THROW_ASSERTION
+{
+ typename VecX::const_iterator xi;
+ for (xi = x.begin(); xi != x.end(); ++xi)
+ y[xi.index()] = *xi;
+}
+
+
+/* perform a gather JGS, does this really make sense? */
+template <class VecX, class VecY>
+inline void
+oned_copy(const VecX& x, VecY y, dense_tag, sparse_tag) MTL_THROW_ASSERTION
+{
+ typedef typename VecX::value_type T;
+ typename VecY::iterator yi;
+ for (yi = y.begin(); yi != y.end(); ++yi)
+ *yi = x[yi.index()];
+}
+#else
+
+
+template <class VecX, class VecY>
+inline void
+oned_copy(const VecX& x, VecY& y, sparse_tag, dense_tag) MTL_THROW_ASSERTION
+{
+ typedef typename linalg_traits<VecY>::value_type T;
+ mtl::set_value(y, T(0));
+ typename VecX::const_iterator xi;
+ for (xi = x.begin(); xi != x.end(); ++xi)
+ y[xi.index()] = *xi;
+}
+
+
+//: Scatter <tt>y <- x</tt>
+//
+// Scatters the elements of the sparse vector x into
+// the dense vector y.
+//
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+//!complexity: O(n) where n is the size of the sparse vector
+template <class VecX, class VecY>
+inline void
+scatter(const VecX& x, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION
+{
+ VecY& y = const_cast<VecY&>(y_);
+ typename VecX::const_iterator xi;
+ for (xi = x.begin(); xi != x.end(); ++xi)
+ y[xi.index()] = *xi;
+}
+
+//: Gather <tt>y <- x</tt>
+//
+// Gathers the elements of the dense vector x into
+// the sparse vector y, based on the non-zero structure of y.
+//
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+//!complexity: O(n) where n is the size of the sparse vector
+template <class VecX, class VecY>
+inline void
+gather(const VecX& x, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION
+{
+ VecY& y = const_cast<VecY&>(y_);
+ typedef typename VecX::value_type T;
+ typename VecY::iterator yi;
+ for (yi = y.begin(); yi != y.end(); ++yi)
+ *yi = x[yi.index()];
+}
+#endif
+
+template <class VecX, class VecY, class Tag>
+inline void
+oned_copy(const VecX& x, VecY& y, Tag, sparse_tag) MTL_THROW_ASSERTION
+{
+ MTL_ASSERT(x.size() <= y.size(), "mtl::copy()");
+ y.clear();
+ typename VecX::const_iterator i = x.begin(), iend = x.end();
+ for (; i != iend; ++i)
+ y.push_back(i.index(), *i);
+}
+
+
+template <class VecX, class VecY>
+inline void
+copy__(const VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION
+{
+ typedef typename linalg_traits<VecX>::sparsity SpX;
+ typedef typename linalg_traits<VecY>::sparsity SpY;
+ oned_copy(x, y, SpX(), SpY());
+}
+
+
+template <class MatA, class MatB>
+inline void
+twod_copy_default(const MatA& A, MatB& B) MTL_THROW_ASSERTION
+{
+ typename MatA::const_iterator i;
+ typename MatA::OneD::const_iterator j, jend;
+
+ for (i = A.begin(); i != A.end(); ++i) {
+ j = (*i).begin(); jend = (*i).end();
+ for (; j != jend; ++j)
+ B(j.row(),j.column()) = *j;
+ }
+}
+
+template <class MatA, class MatB>
+inline void
+twod_copy(const MatA& A, MatB& B, rectangle_tag) MTL_THROW_ASSERTION
+{
+ twod_copy_default(A, B);
+}
+
+template <class MatA, class MatB>
+inline void
+twod_copy(const MatA& A, MatB& B, banded_tag) MTL_THROW_ASSERTION
+{
+ twod_copy_default(A, B);
+}
+
+template <class MatA, class MatB>
+inline void
+twod_copy(const MatA& A, MatB& B, symmetric_tag) MTL_THROW_ASSERTION
+{
+ typename MatA::const_iterator i;
+ typename MatA::OneD::const_iterator j, jend;
+
+ for (i = A.begin(); i != A.end(); ++i) {
+ j = (*i).begin(); jend = (*i).end();
+ for (; j != jend; ++j) {
+ B(j.row(),j.column()) = *j;
+ B(j.column(),j.row()) = *j;
+ }
+ }
+}
+
+template <class MatA, class MatB>
+inline void
+twod_copy(const MatA& A, MatB& B, triangle_tag) MTL_THROW_ASSERTION
+{
+ typedef typename matrix_traits<MatB>::value_type T;
+
+ if (A.is_unit())
+ set_diagonal(B, T(1));
+
+ twod_copy(A, B, rectangle_tag());
+}
+
+template <class MatA, class MatB>
+inline void
+twod_copy__(const MatA& A, MatB& B, dense_tag)
+{
+ typedef typename matrix_traits<MatA>::shape Shape;
+ twod_copy(A, B, Shape());
+}
+
+
+/*
+ Sparse matrices have specialized copy functions since
+ they need to optimize the creation of the non-zero structure.
+
+ only good for same orientation!!!
+ */
+
+
+template <class MatA, class MatB>
+inline void
+twod_copy(const MatA& A, MatB& B, row_tag, row_tag)
+{
+ B.fast_copy(A);
+}
+template <class MatA, class MatB>
+inline void
+twod_copy(const MatA& A, MatB& B, column_tag, column_tag)
+{
+ B.fast_copy(A);
+}
+
+template <class MatA, class MatB>
+inline void
+twod_copy(const MatA& A, MatB& B, row_tag, column_tag)
+{
+ twod_copy__(A, B, dense_tag());
+}
+template <class MatA, class MatB>
+inline void
+twod_copy(const MatA& A, MatB& B, column_tag, row_tag)
+{
+ twod_copy__(A, B, dense_tag());
+}
+
+template <class MatA, class MatB>
+inline void
+twod_copy__(const MatA& A, MatB& B, sparse_tag)
+{
+ typedef typename matrix_traits<MatA>::orientation OrienA;
+ typedef typename matrix_traits<MatB>::orientation OrienB;
+ twod_copy(A, B, OrienA(), OrienB());
+}
+
+template <class MatA, class MatB>
+inline void
+copy__(const MatA& A, MatB& B, twod_tag) MTL_THROW_ASSERTION
+{
+ MTL_ASSERT(A.nrows() <= B.nrows(), "copy(A, B, twod_tag)");
+ MTL_ASSERT(A.ncols() <= B.ncols(), "copy(A, B, twod_tag)");
+
+ typedef typename matrix_traits<MatB>::sparsity Sparsity;
+ twod_copy__(A, B, Sparsity());
+}
+
+//: Copy: <tt>B <- A or y <- x</tt>
+//
+// Copy the elements of matrix A into matrix B, or copy the elements
+// of vector x into vector y. For shaped and sparse matrices, this
+// copies only the elements stored in A to B. If x is a sparse
+// vector and y is dense, a "scatter" is performed. If y is sparse
+// and x is dense, then a "gather" is performed. If both vectors
+// are sparse, but of different structure the result is undefined.
+//
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+//!complexity: O(m*n) for matrices. O(nnz) if either A or B are sparse and of the same orientation (otherwise it can be O(nnz^2). O(n) for vectors.
+//!example: vecvec_copy.cc
+template <class LinalgA, class LinalgB>
+inline void
+copy(const LinalgA& A, MTL_OUT(LinalgB) B_) MTL_THROW_ASSERTION
+{
+ LinalgB& B = const_cast<LinalgB&>(B_);
+ typedef typename linalg_traits<LinalgA>::dimension Dim;
+ copy__(A, B, Dim());
+}
+
+template <class VecX, class VecY> inline
+void
+add__(const VecX& x, VecY& y, fast::count<0>)
+{
+ typedef typename VecX::value_type T;
+ mtl_algo::transform_add(x.begin(), x.end(), y.begin());
+
+}
+#if USE_BLAIS
+template <class VecX, class VecY, int N> inline
+void
+add__(const VecX& x, VecY& y, fast::count<N>)
+{
+ typedef typename VecX::value_type T;
+ fast::transform(x.begin(), fast::count<N>(), y.begin(),
+ y.begin(), std::plus<T>());
+}
+#endif
+template <class VecX, class VecY> inline
+void
+add__(const VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION
+{
+ MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
+
+ add__(x, y, dim_n<VecX>::RET());
+}
+
+
+template <class VecX, class VecY, class VecZ> inline
+void
+oned_add(const VecX& x, const VecY& y, VecZ& z, fast::count<0>)
+{
+ typedef typename VecX::value_type T;
+ mtl_algo::transform(x.begin(), x.end(), y.begin(), z.begin(), std::plus<T>());
+}
+#if USE_BLAIS
+template <class VecX, class VecY, class VecZ, int N> inline
+void
+oned_add(const VecX& x, const VecY& y, VecZ& z, fast::count<N>)
+{
+ typedef typename VecX::value_type T;
+ fast::transform(x.begin(), fast::count<N>(), y.begin(), z.begin(), std::plus<T>());
+}
+#endif
+
+template <class VecX, class VecY, class VecZ>
+inline void
+oned_add(const VecX& x, const VecY& y, VecZ& z_, sparse_tag)
+{
+
+ typedef typename VecZ::value_type T;
+ compressed1D<T> z;
+ typedef typename VecX::const_iterator xiter;
+ typedef typename VecY::const_iterator yiter;
+
+ xiter xi = x.begin();
+ xiter xiend = x.end();
+ yiter yi = y.begin();
+ yiter yiend = y.end();
+
+ while (xi != xiend && yi != yiend) {
+ if (yi.index() < xi.index()) {
+ z.push_back(yi.index(), *yi);
+ ++yi;
+ } else if (xi.index() < yi.index()) {
+ z.push_back(xi.index(), *xi);
+ ++xi;
+ } else {
+ z.push_back(xi.index(), *yi + *xi);
+ ++xi; ++yi;
+ }
+ }
+ while (xi != xiend) {
+ z.push_back(xi.index(), *xi);
+ ++xi;
+ }
+ while (yi != yiend) {
+ z.push_back(yi.index(), *yi);
+ ++yi;
+ }
+ z_.clear();
+ mtl::copy(z, z_);
+}
+
+template <class VecX, class VecY, class VecZ>
+inline void
+oned_add(const VecX& x, const VecY& y, VecZ& z, dense_tag) MTL_THROW_ASSERTION
+{
+ oned_add(x, y, z, dim_n<VecX>::RET());
+}
+
+
+//: Add: <tt>z <- x + y</tt>
+//
+// Add the elements of x and y and assign into z.
+//
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+//!example: y_ax_y.cc, vecvec_add.cc
+//!typereqs: <tt>VecX::value_type</tt>, <tt>VecY::value_type</tt>, and <tt>VecZ::value_type</tt> should be the same type
+//!typereqs: The addition operator must be defined for the value_type.
+//!complexity: linear time
+template <class VecX, class VecY, class VecZ>
+inline void
+add(const VecX& x, const VecY& y, MTL_OUT(VecZ) z_) MTL_THROW_ASSERTION
+{
+ VecZ& z = const_cast<VecZ&>(z_);
+ MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
+ MTL_ASSERT(x.size() <= z.size(), "mtl::add()");
+ typedef typename linalg_traits<VecZ>::sparsity Sparsity;
+ oned_add(x, y, z, Sparsity());
+}
+
+//: Add: <tt>w <- x + y + z</tt>
+//
+// Add the elements of x, y, and z and assign into w.
+// For now just dense vectors.
+//
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+//!example: vecvec_add3.cc
+//!typereqs: <tt>VecX::value_type</tt>, <tt>VecY::value_type</tt>, <tt>VecZ::value_type</tt>, and <tt>VecW::value_type</tt> should be the same type
+//!typereqs: The addition operator must be defined for the value_type.
+//!complexity: linear time
+template <class VecW, class VecX, class VecY, class VecZ>
+inline void
+add(const VecX& x, const VecY& y, const VecZ& z, MTL_OUT(VecW) w_)
+ MTL_THROW_ASSERTION
+{
+ VecW& w = const_cast<VecW&>(w_);
+ MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
+ MTL_ASSERT(x.size() <= z.size(), "mtl::add()");
+ MTL_ASSERT(x.size() <= w.size(), "mtl::add()");
+
+ typename VecX::const_iterator x_i = x.begin();
+ typename VecY::const_iterator y_i = y.begin();
+ typename VecZ::const_iterator z_i = z.begin();
+ typename VecW::iterator w_i = w.begin();
+
+ while (not_at(x_i, x.end())) {
+ *w_i = *x_i + *y_i + *z_i;
+ ++x_i; ++y_i; ++z_i; ++w_i;
+ }
+}
+
+
+template <class MatA, class MatB>
+inline void
+twod_add_default(const MatA& A, MatB& B)
+{
+ typename MatA::const_iterator i;
+ typename MatA::OneD::const_iterator j, jend;
+
+ for (i = A.begin(); i != A.end(); ++i) {
+ j = (*i).begin(); jend = (*i).end();
+ for (; j != jend; ++j)
+ B(j.row(), j.column()) += *j;
+ }
+}
+
+template <class MatA, class MatB>
+inline void
+twod_add(const MatA& A, MatB& B, banded_tag)
+{
+ twod_add_default(A, B);
+}
+
+template <class MatA, class MatB>
+inline void
+twod_add(const MatA& A, MatB& B, rectangle_tag)
+{
+ twod_add_default(A, B);
+}
+
+template <class MatA, class MatB>
+inline void
+twod_add(const MatA& A, MatB& B, triangle_tag)
+{
+ typedef typename matrix_traits<MatA>::size_type Int;
+ typedef typename matrix_traits<MatA>::value_type T;
+ if (A.is_unit())
+ for (Int i = 0; i < MTL_MIN(A.nrows(), A.ncols()); ++i)
+ B(i,i) += T(1);
+
+ twod_add(A, B, banded_tag());
+}
+
+/* perhaps I should add is_row() and is_column()
+ methods to the matrices
+ */
+template <class MatA, class MatB>
+inline void
+twod_symmetric_add(const MatA& A, MatB& B, row_tag)
+{
+ typename MatA::const_iterator i;
+ typename MatA::Row::const_iterator j, jend;
+
+ for (i = A.begin(); i != A.end(); ++i) {
+ j = (*i).begin();
+ jend = (*i).end();
+ if (A.is_upper()) { /* handle the diagonal elements */
+ B(j.column(), j.row()) += *j;
+ ++j;
+ } else
+ --jend;
+ for (; j != jend; ++j) {
+ B(j.row(), j.column()) += *j;
+ B(j.column(), j.row()) += *j;
+ }
+ if (A.is_lower())
+ B(j.column(), j.row()) += *j;
+ }
+}
+
+template <class MatA, class MatB>
+inline void
+twod_symmetric_add(const MatA& A, MatB& B, column_tag)
+{
+ typename MatA::const_iterator i;
+ typename MatA::Column::const_iterator j, jend;
+
+ for (i = A.begin(); i != A.end(); ++i) {
+ j = (*i).begin();
+ jend = (*i).end();
+ if (A.is_lower()) { /* handle the diagonal elements */
+ B(j.column(), j.row()) += *j;
+ ++j;
+ } else
+ --jend;
+ for (; j != jend; ++j) {
+ B(j.row(), j.column()) += *j;
+ B(j.column(), j.row()) += *j;
+ }
+ if (A.is_upper())
+ B(j.column(), j.row()) += *j;
+ }
+}
+
+
+template <class MatA, class MatB>
+inline void
+twod_add(const MatA& A, MatB& B, symmetric_tag)
+{
+ typedef typename matrix_traits<MatA>::orientation Orien;
+ twod_symmetric_add(A, B, Orien());
+}
+
+
+template <class MatA, class MatB>
+inline void
+add__(const MatA& A, MatB& B, twod_tag) MTL_THROW_ASSERTION
+{
+ MTL_ASSERT(A.nrows() <= B.nrows(), "matmat::add()");
+ MTL_ASSERT(A.ncols() <= B.ncols(), "matmat::add()");
+
+ typedef typename matrix_traits<MatA>::shape Shape;
+ twod_add(A, B, Shape());
+}
+
+//: Add: <tt>B <- A + B or y <- x + y</tt>
+// The function adds the element of A to B, or the elements of x to y.
+//
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+//!complexity: O(m*n) for a dense A, O(nnz) for a sparse A. O(n) for a vector.
+
+template <class LinalgA, class LinalgB>
+inline void
+add(const LinalgA& A, MTL_OUT(LinalgB) B_) MTL_THROW_ASSERTION
+{
+ LinalgB& B = const_cast<LinalgB&>(B_);
+ typedef typename linalg_traits<LinalgA>::dimension Dim;
+ add__(A, B, Dim());
+}
+
+
+template <class VecX, class VecY, class T>
+inline T
+dot(const VecX& x, const VecY& y, T s, fast::count<0>)
+{
+ return mtl_algo::inner_product(x.begin(), x.end(), y.begin(), s);
+}
+#if USE_BLAIS
+template <class VecX, class VecY, class T, int N>
+inline T
+dot(const VecX& x, const VecY& y, T s, fast::count<N>)
+{
+ return fast::inner_product(x.begin(), fast::count<N>(), y.begin(), s);
+}
+#endif
+
+
+template <class VecX, class VecY, class T>
+inline T
+dot(const VecX& x, const VecY& y, T s, dense_tag, dense_tag)
+{
+ return dot(x, y, s, dim_n<VecX>::RET());
+}
+
+template <class InputIterator1, class InputIterator2, class T>
+inline T
+sparse_inner_product(InputIterator1 f1, InputIterator1 l1,
+ InputIterator2 f2, InputIterator2 l2, T init)
+{
+ InputIterator1 first1 = f1;
+ InputIterator1 last1 = l1;
+ InputIterator2 first2 = f2;
+ InputIterator2 last2 = l2;
+
+ while (first1 != last1 && first2 != last2) {
+ if (first1.index() == first2.index())
+ init += (*first1++ * *first2++);
+ else if (first1.index() < first2.index())
+ ++first1;
+ else
+ ++first2;
+ }
+ return init;
+}
+
+template <class IndexedIterator, class RandomAccessIterator, class T>
+inline T
+sparse_dense_inner_product(IndexedIterator f1, IndexedIterator l1,
+ RandomAccessIterator f2, T init)
+{
+ IndexedIterator first1 = f1, last1 = l1;
+ RandomAccessIterator first2 = f2;
+
+ while (first1 != last1) {
+ init += (*first1 * first2[first1.index()]);
+ ++first1;
+ }
+ return init;
+}
+
+
+template <class VecX, class VecY, class T>
+inline T
+dot(const VecX& x, const VecY& y, T s, sparse_tag, sparse_tag)
+{
+ if (x.nnz() < y.nnz())
+ return sparse_inner_product(x.begin(), x.end(), y.begin(), y.end(), s);
+ else
+ return sparse_inner_product(y.begin(), y.end(), x.begin(), x.end(), s);
+}
+
+template <class VecX, class VecY, class T>
+inline T
+dot(const VecX& x, const VecY& y, T s, dense_tag, sparse_tag)
+{
+ return sparse_dense_inner_product(y.begin(), y.end(), x.begin(), s);
+}
+
+template <class VecX, class VecY, class T>
+inline T
+dot(const VecX& x, const VecY& y, T s, sparse_tag, dense_tag)
+{
+ return sparse_dense_inner_product(x.begin(), x.end(), y.begin(), s);
+}
+
+
+//: Dot Product: <tt>s <- x . y + s</tt>
+// The type used for argument s determines the
+// type of the resulting product.
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+template <class VecX, class VecY, class T>
+inline T
+dot(const VecX& x, const VecY& y, T s) MTL_THROW_ASSERTION
+{
+ MTL_ASSERT(x.size() == y.size(), "mtl::dot()");
+ typedef typename linalg_traits<VecX>::sparsity SparseX;
+ typedef typename linalg_traits<VecY>::sparsity SparseY;
+ return dot(x, y, s, SparseX(), SparseY());
+}
+
+
+//: Dot Product: <tt>s <- x . y</tt>
+// The type of the resulting product is <TT>VecX::value_type</TT>.
+//!category: algorithms
+//!component: function
+//!example: vecvec_dot.cc, dot_prod.cc
+//!definition: mtl.h
+template <class VecX, class VecY>
+inline typename VecX::value_type
+dot(const VecX& x, const VecY& y) MTL_THROW_ASSERTION
+{
+ typedef typename VecX::value_type T;
+ return mtl::dot(x, y, T(0));
+}
+
+#ifdef USE_DOUBLE_DOUBLE
+//: Dot Product (extended precision): <tt>s <- x . y + s</tt>
+// The type of the resulting product is double_double
+// Extended precision is used internally.
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+template <class VecX, class VecY>
+inline double_double
+dot(const VecX& x, const VecY& y, double_double s) MTL_THROW_ASSERTION
+{
+ typedef typename VecX::value_type x_type;
+ typedef typename VecY::value_type y_type;
+ // x_type and y_type must be either float or double
+ multiply<double_double, x_type, y_type> m;
+ addition<double_double, double_double, double_double> a;
+ return mtl_algo::inner_product(x.begin(), x.end(), y.begin(), s, a, m);
+}
+#endif /* USE_DOUBLE_DOUBLE */
+
+template <class T>
+struct conj_func {
+ typedef T result_type;
+ inline T operator()(const T& x) const { return MTL_CONJ(x); }
+};
+
+template <class VecX, class VecY, class T>
+inline T
+dot_conj(const VecX& x, const VecY& y, T s, fast::count<0>)
+{
+ return mtl_algo::inner_product(x.begin(), x.end(),
+ trans_iter(y.begin(), conj_func<T>()), s);
+}
+#if USE_BLAIS
+template <class VecX, class VecY, class T, int N>
+inline T
+dot_conj(const VecX& x, const VecY& y, T s, fast::count<N>)
+{
+ return fast::inner_product(x.begin(), x.end(),
+ trans_iter(y.begin(), conj_func<T>()), s);
+}
+#endif
+
+//: Dot Conjugate: <tt>s <- x . conj(y) + s</tt>
+// Similar to dot product. The complex conjugate of the elements of y
+// is used. For real numbers, the conjugate is just that real number.
+// Note that the type of parameter s is the return type of this
+// function.
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+template <class VecX, class VecY, class T>
+inline T
+dot_conj(const VecX& x, const VecY& y, T s) MTL_THROW_ASSERTION
+{
+ MTL_ASSERT(x.size() <= y.size(), "mtl::dot_conj()");
+ return dot_conj(x, y, s, dim_n<VecX>::RET());
+}
+
+//: Dot Conjugate: <tt>s <- x . conj(y)</tt>
+// A slightly simpler version of the dot conjugate.
+// The return type is the element type of vector x.
+//!category: algorithms
+//!component: function
+//!definition: mtl.h
+template <class VecX, class VecY>
+inline typename VecX::value_type
+dot_conj(const VecX& x, const VecY& y) MTL_THROW_ASSERTION
+{
+ typedef typename VecX::value_type T;
+ return mtl::dot_conj(x, y, T(0));
+}
+
/*
This version of the algorithm depends on the compiler
hoisting the reference of z[j.row()] out of the inner loop
@@ -2424,551 +3157,6 @@
}
}
-template <class VecX, class VecY>
-inline void
-copy__(const VecX& x, VecY& y, fast::count<0>)
-{
- mtl_algo::copy(x.begin(), x.end(), y.begin());
-}
-#if USE_BLAIS
-template <class VecX, class VecY, int N>
-inline void
-copy__(const VecX& x, VecY& y, fast::count<N>)
-{
- fast::copy(x.begin(), fast::count<N>(), y.begin());
-}
-#endif
-
-
-template <class VecX, class VecY>
-inline void
-oned_copy(const VecX& x, VecY& y, dense_tag, dense_tag) MTL_THROW_ASSERTION
-{
- MTL_ASSERT(x.size() <= y.size(), "mtl::copy()");
- copy__(x, y, dim_n<VecX>::RET());
-}
-
-#if 0
-/* perform a scatter */
-template <class VecX, class VecY>
-inline void
-oned_copy(const VecX& x, VecY y, sparse_tag, dense_tag) MTL_THROW_ASSERTION
-{
- typename VecX::const_iterator xi;
- for (xi = x.begin(); xi != x.end(); ++xi)
- y[xi.index()] = *xi;
-}
-
-
-/* perform a gather JGS, does this really make sense? */
-template <class VecX, class VecY>
-inline void
-oned_copy(const VecX& x, VecY y, dense_tag, sparse_tag) MTL_THROW_ASSERTION
-{
- typedef typename VecX::value_type T;
- typename VecY::iterator yi;
- for (yi = y.begin(); yi != y.end(); ++yi)
- *yi = x[yi.index()];
-}
-#else
-
-
-template <class VecX, class VecY>
-inline void
-oned_copy(const VecX& x, VecY& y, sparse_tag, dense_tag) MTL_THROW_ASSERTION
-{
- typedef typename linalg_traits<VecY>::value_type T;
- mtl::set_value(y, T(0));
- typename VecX::const_iterator xi;
- for (xi = x.begin(); xi != x.end(); ++xi)
- y[xi.index()] = *xi;
-}
-
-
-//: Scatter <tt>y <- x</tt>
-//
-// Scatters the elements of the sparse vector x into
-// the dense vector y.
-//
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-//!complexity: O(n) where n is the size of the sparse vector
-template <class VecX, class VecY>
-inline void
-scatter(const VecX& x, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION
-{
- VecY& y = const_cast<VecY&>(y_);
- typename VecX::const_iterator xi;
- for (xi = x.begin(); xi != x.end(); ++xi)
- y[xi.index()] = *xi;
-}
-
-//: Gather <tt>y <- x</tt>
-//
-// Gathers the elements of the dense vector x into
-// the sparse vector y, based on the non-zero structure of y.
-//
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-//!complexity: O(n) where n is the size of the sparse vector
-template <class VecX, class VecY>
-inline void
-gather(const VecX& x, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION
-{
- VecY& y = const_cast<VecY&>(y_);
- typedef typename VecX::value_type T;
- typename VecY::iterator yi;
- for (yi = y.begin(); yi != y.end(); ++yi)
- *yi = x[yi.index()];
-}
-#endif
-
-template <class VecX, class VecY, class Tag>
-inline void
-oned_copy(const VecX& x, VecY& y, Tag, sparse_tag) MTL_THROW_ASSERTION
-{
- MTL_ASSERT(x.size() <= y.size(), "mtl::copy()");
- y.clear();
- typename VecX::const_iterator i = x.begin(), iend = x.end();
- for (; i != iend; ++i)
- y.push_back(i.index(), *i);
-}
-
-
-template <class VecX, class VecY>
-inline void
-copy__(const VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION
-{
- typedef typename linalg_traits<VecX>::sparsity SpX;
- typedef typename linalg_traits<VecY>::sparsity SpY;
- oned_copy(x, y, SpX(), SpY());
-}
-
-
-template <class MatA, class MatB>
-inline void
-twod_copy_default(const MatA& A, MatB& B) MTL_THROW_ASSERTION
-{
- typename MatA::const_iterator i;
- typename MatA::OneD::const_iterator j, jend;
-
- for (i = A.begin(); i != A.end(); ++i) {
- j = (*i).begin(); jend = (*i).end();
- for (; j != jend; ++j)
- B(j.row(),j.column()) = *j;
- }
-}
-
-template <class MatA, class MatB>
-inline void
-twod_copy(const MatA& A, MatB& B, rectangle_tag) MTL_THROW_ASSERTION
-{
- twod_copy_default(A, B);
-}
-
-template <class MatA, class MatB>
-inline void
-twod_copy(const MatA& A, MatB& B, banded_tag) MTL_THROW_ASSERTION
-{
- twod_copy_default(A, B);
-}
-
-template <class MatA, class MatB>
-inline void
-twod_copy(const MatA& A, MatB& B, symmetric_tag) MTL_THROW_ASSERTION
-{
- typename MatA::const_iterator i;
- typename MatA::OneD::const_iterator j, jend;
-
- for (i = A.begin(); i != A.end(); ++i) {
- j = (*i).begin(); jend = (*i).end();
- for (; j != jend; ++j) {
- B(j.row(),j.column()) = *j;
- B(j.column(),j.row()) = *j;
- }
- }
-}
-
-template <class MatA, class MatB>
-inline void
-twod_copy(const MatA& A, MatB& B, triangle_tag) MTL_THROW_ASSERTION
-{
- typedef typename matrix_traits<MatB>::value_type T;
-
- if (A.is_unit())
- set_diagonal(B, T(1));
-
- twod_copy(A, B, rectangle_tag());
-}
-
-template <class MatA, class MatB>
-inline void
-twod_copy__(const MatA& A, MatB& B, dense_tag)
-{
- typedef typename matrix_traits<MatA>::shape Shape;
- twod_copy(A, B, Shape());
-}
-
-
-/*
- Sparse matrices have specialized copy functions since
- they need to optimize the creation of the non-zero structure.
-
- only good for same orientation!!!
- */
-
-
-template <class MatA, class MatB>
-inline void
-twod_copy(const MatA& A, MatB& B, row_tag, row_tag)
-{
- B.fast_copy(A);
-}
-template <class MatA, class MatB>
-inline void
-twod_copy(const MatA& A, MatB& B, column_tag, column_tag)
-{
- B.fast_copy(A);
-}
-
-template <class MatA, class MatB>
-inline void
-twod_copy(const MatA& A, MatB& B, row_tag, column_tag)
-{
- twod_copy__(A, B, dense_tag());
-}
-template <class MatA, class MatB>
-inline void
-twod_copy(const MatA& A, MatB& B, column_tag, row_tag)
-{
- twod_copy__(A, B, dense_tag());
-}
-
-template <class MatA, class MatB>
-inline void
-twod_copy__(const MatA& A, MatB& B, sparse_tag)
-{
- typedef typename matrix_traits<MatA>::orientation OrienA;
- typedef typename matrix_traits<MatB>::orientation OrienB;
- twod_copy(A, B, OrienA(), OrienB());
-}
-
-template <class MatA, class MatB>
-inline void
-copy__(const MatA& A, MatB& B, twod_tag) MTL_THROW_ASSERTION
-{
- MTL_ASSERT(A.nrows() <= B.nrows(), "copy(A, B, twod_tag)");
- MTL_ASSERT(A.ncols() <= B.ncols(), "copy(A, B, twod_tag)");
-
- typedef typename matrix_traits<MatB>::sparsity Sparsity;
- twod_copy__(A, B, Sparsity());
-}
-
-//: Copy: <tt>B <- A or y <- x</tt>
-//
-// Copy the elements of matrix A into matrix B, or copy the elements
-// of vector x into vector y. For shaped and sparse matrices, this
-// copies only the elements stored in A to B. If x is a sparse
-// vector and y is dense, a "scatter" is performed. If y is sparse
-// and x is dense, then a "gather" is performed. If both vectors
-// are sparse, but of different structure the result is undefined.
-//
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-//!complexity: O(m*n) for matrices. O(nnz) if either A or B are sparse and of the same orientation (otherwise it can be O(nnz^2). O(n) for vectors.
-//!example: vecvec_copy.cc
-template <class LinalgA, class LinalgB>
-inline void
-copy(const LinalgA& A, MTL_OUT(LinalgB) B_) MTL_THROW_ASSERTION
-{
- LinalgB& B = const_cast<LinalgB&>(B_);
- typedef typename linalg_traits<LinalgA>::dimension Dim;
- copy__(A, B, Dim());
-}
-
-template <class VecX, class VecY> inline
-void
-add__(const VecX& x, VecY& y, fast::count<0>)
-{
- typedef typename VecX::value_type T;
- mtl_algo::transform_add(x.begin(), x.end(), y.begin());
-
-}
-#if USE_BLAIS
-template <class VecX, class VecY, int N> inline
-void
-add__(const VecX& x, VecY& y, fast::count<N>)
-{
- typedef typename VecX::value_type T;
- fast::transform(x.begin(), fast::count<N>(), y.begin(),
- y.begin(), std::plus<T>());
-}
-#endif
-template <class VecX, class VecY> inline
-void
-add__(const VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION
-{
- MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
-
- add__(x, y, dim_n<VecX>::RET());
-}
-
-
-template <class VecX, class VecY, class VecZ> inline
-void
-oned_add(const VecX& x, const VecY& y, VecZ& z, fast::count<0>)
-{
- typedef typename VecX::value_type T;
- mtl_algo::transform(x.begin(), x.end(), y.begin(), z.begin(), std::plus<T>());
-}
-#if USE_BLAIS
-template <class VecX, class VecY, class VecZ, int N> inline
-void
-oned_add(const VecX& x, const VecY& y, VecZ& z, fast::count<N>)
-{
- typedef typename VecX::value_type T;
- fast::transform(x.begin(), fast::count<N>(), y.begin(), z.begin(), std::plus<T>());
-}
-#endif
-
-template <class VecX, class VecY, class VecZ>
-inline void
-oned_add(const VecX& x, const VecY& y, VecZ& z_, sparse_tag)
-{
-
- typedef typename VecZ::value_type T;
- compressed1D<T> z;
- typedef typename VecX::const_iterator xiter;
- typedef typename VecY::const_iterator yiter;
-
- xiter xi = x.begin();
- xiter xiend = x.end();
- yiter yi = y.begin();
- yiter yiend = y.end();
-
- while (xi != xiend && yi != yiend) {
- if (yi.index() < xi.index()) {
- z.push_back(yi.index(), *yi);
- ++yi;
- } else if (xi.index() < yi.index()) {
- z.push_back(xi.index(), *xi);
- ++xi;
- } else {
- z.push_back(xi.index(), *yi + *xi);
- ++xi; ++yi;
- }
- }
- while (xi != xiend) {
- z.push_back(xi.index(), *xi);
- ++xi;
- }
- while (yi != yiend) {
- z.push_back(yi.index(), *yi);
- ++yi;
- }
- z_.clear();
- mtl::copy(z, z_);
-}
-
-template <class VecX, class VecY, class VecZ>
-inline void
-oned_add(const VecX& x, const VecY& y, VecZ& z, dense_tag) MTL_THROW_ASSERTION
-{
- oned_add(x, y, z, dim_n<VecX>::RET());
-}
-
-
-//: Add: <tt>z <- x + y</tt>
-//
-// Add the elements of x and y and assign into z.
-//
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-//!example: y_ax_y.cc, vecvec_add.cc
-//!typereqs: <tt>VecX::value_type</tt>, <tt>VecY::value_type</tt>, and <tt>VecZ::value_type</tt> should be the same type
-//!typereqs: The addition operator must be defined for the value_type.
-//!complexity: linear time
-template <class VecX, class VecY, class VecZ>
-inline void
-add(const VecX& x, const VecY& y, MTL_OUT(VecZ) z_) MTL_THROW_ASSERTION
-{
- VecZ& z = const_cast<VecZ&>(z_);
- MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
- MTL_ASSERT(x.size() <= z.size(), "mtl::add()");
- typedef typename linalg_traits<VecZ>::sparsity Sparsity;
- oned_add(x, y, z, Sparsity());
-}
-
-//: Add: <tt>w <- x + y + z</tt>
-//
-// Add the elements of x, y, and z and assign into w.
-// For now just dense vectors.
-//
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-//!example: vecvec_add3.cc
-//!typereqs: <tt>VecX::value_type</tt>, <tt>VecY::value_type</tt>, <tt>VecZ::value_type</tt>, and <tt>VecW::value_type</tt> should be the same type
-//!typereqs: The addition operator must be defined for the value_type.
-//!complexity: linear time
-template <class VecW, class VecX, class VecY, class VecZ>
-inline void
-add(const VecX& x, const VecY& y, const VecZ& z, MTL_OUT(VecW) w_)
- MTL_THROW_ASSERTION
-{
- VecW& w = const_cast<VecW&>(w_);
- MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
- MTL_ASSERT(x.size() <= z.size(), "mtl::add()");
- MTL_ASSERT(x.size() <= w.size(), "mtl::add()");
-
- typename VecX::const_iterator x_i = x.begin();
- typename VecY::const_iterator y_i = y.begin();
- typename VecZ::const_iterator z_i = z.begin();
- typename VecW::iterator w_i = w.begin();
-
- while (not_at(x_i, x.end())) {
- *w_i = *x_i + *y_i + *z_i;
- ++x_i; ++y_i; ++z_i; ++w_i;
- }
-}
-
-
-template <class MatA, class MatB>
-inline void
-twod_add_default(const MatA& A, MatB& B)
-{
- typename MatA::const_iterator i;
- typename MatA::OneD::const_iterator j, jend;
-
- for (i = A.begin(); i != A.end(); ++i) {
- j = (*i).begin(); jend = (*i).end();
- for (; j != jend; ++j)
- B(j.row(), j.column()) += *j;
- }
-}
-
-template <class MatA, class MatB>
-inline void
-twod_add(const MatA& A, MatB& B, banded_tag)
-{
- twod_add_default(A, B);
-}
-
-template <class MatA, class MatB>
-inline void
-twod_add(const MatA& A, MatB& B, rectangle_tag)
-{
- twod_add_default(A, B);
-}
-
-template <class MatA, class MatB>
-inline void
-twod_add(const MatA& A, MatB& B, triangle_tag)
-{
- typedef typename matrix_traits<MatA>::size_type Int;
- typedef typename matrix_traits<MatA>::value_type T;
- if (A.is_unit())
- for (Int i = 0; i < MTL_MIN(A.nrows(), A.ncols()); ++i)
- B(i,i) += T(1);
-
- twod_add(A, B, banded_tag());
-}
-
-/* perhaps I should add is_row() and is_column()
- methods to the matrices
- */
-template <class MatA, class MatB>
-inline void
-twod_symmetric_add(const MatA& A, MatB& B, row_tag)
-{
- typename MatA::const_iterator i;
- typename MatA::Row::const_iterator j, jend;
-
- for (i = A.begin(); i != A.end(); ++i) {
- j = (*i).begin();
- jend = (*i).end();
- if (A.is_upper()) { /* handle the diagonal elements */
- B(j.column(), j.row()) += *j;
- ++j;
- } else
- --jend;
- for (; j != jend; ++j) {
- B(j.row(), j.column()) += *j;
- B(j.column(), j.row()) += *j;
- }
- if (A.is_lower())
- B(j.column(), j.row()) += *j;
- }
-}
-
-template <class MatA, class MatB>
-inline void
-twod_symmetric_add(const MatA& A, MatB& B, column_tag)
-{
- typename MatA::const_iterator i;
- typename MatA::Column::const_iterator j, jend;
-
- for (i = A.begin(); i != A.end(); ++i) {
- j = (*i).begin();
- jend = (*i).end();
- if (A.is_lower()) { /* handle the diagonal elements */
- B(j.column(), j.row()) += *j;
- ++j;
- } else
- --jend;
- for (; j != jend; ++j) {
- B(j.row(), j.column()) += *j;
- B(j.column(), j.row()) += *j;
- }
- if (A.is_upper())
- B(j.column(), j.row()) += *j;
- }
-}
-
-
-template <class MatA, class MatB>
-inline void
-twod_add(const MatA& A, MatB& B, symmetric_tag)
-{
- typedef typename matrix_traits<MatA>::orientation Orien;
- twod_symmetric_add(A, B, Orien());
-}
-
-
-template <class MatA, class MatB>
-inline void
-add__(const MatA& A, MatB& B, twod_tag) MTL_THROW_ASSERTION
-{
- MTL_ASSERT(A.nrows() <= B.nrows(), "matmat::add()");
- MTL_ASSERT(A.ncols() <= B.ncols(), "matmat::add()");
-
- typedef typename matrix_traits<MatA>::shape Shape;
- twod_add(A, B, Shape());
-}
-
-//: Add: <tt>B <- A + B or y <- x + y</tt>
-// The function adds the element of A to B, or the elements of x to y.
-//
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-//!complexity: O(m*n) for a dense A, O(nnz) for a sparse A. O(n) for a vector.
-
-template <class LinalgA, class LinalgB>
-inline void
-add(const LinalgA& A, MTL_OUT(LinalgB) B_) MTL_THROW_ASSERTION
-{
- LinalgB& B = const_cast<LinalgB&>(B_);
- typedef typename linalg_traits<LinalgA>::dimension Dim;
- add__(A, B, Dim());
-}
-
-
-
template <class VecX, class VecY, class VecZ>
inline void
ele_mult(const VecX& x, const VecY& y, MTL_OUT(VecZ) z_, fast::count<0>)
@@ -3131,200 +3319,6 @@
typedef typename linalg_traits<LinalgA>::dimension Dim;
swap(const_cast<LinalgA&>(A), const_cast<LinalgB&>(B), Dim());
}
-
-
-template <class VecX, class VecY, class T>
-inline T
-dot(const VecX& x, const VecY& y, T s, fast::count<0>)
-{
- return mtl_algo::inner_product(x.begin(), x.end(), y.begin(), s);
-}
-#if USE_BLAIS
-template <class VecX, class VecY, class T, int N>
-inline T
-dot(const VecX& x, const VecY& y, T s, fast::count<N>)
-{
- return fast::inner_product(x.begin(), fast::count<N>(), y.begin(), s);
-}
-#endif
-
-
-template <class VecX, class VecY, class T>
-inline T
-dot(const VecX& x, const VecY& y, T s, dense_tag, dense_tag)
-{
- return dot(x, y, s, dim_n<VecX>::RET());
-}
-
-template <class InputIterator1, class InputIterator2, class T>
-inline T
-sparse_inner_product(InputIterator1 f1, InputIterator1 l1,
- InputIterator2 f2, InputIterator2 l2, T init)
-{
- InputIterator1 first1 = f1;
- InputIterator1 last1 = l1;
- InputIterator2 first2 = f2;
- InputIterator2 last2 = l2;
-
- while (first1 != last1 && first2 != last2) {
- if (first1.index() == first2.index())
- init += (*first1++ * *first2++);
- else if (first1.index() < first2.index())
- ++first1;
- else
- ++first2;
- }
- return init;
-}
-
-template <class IndexedIterator, class RandomAccessIterator, class T>
-inline T
-sparse_dense_inner_product(IndexedIterator f1, IndexedIterator l1,
- RandomAccessIterator f2, T init)
-{
- IndexedIterator first1 = f1, last1 = l1;
- RandomAccessIterator first2 = f2;
-
- while (first1 != last1) {
- init += (*first1 * first2[first1.index()]);
- ++first1;
- }
- return init;
-}
-
-
-template <class VecX, class VecY, class T>
-inline T
-dot(const VecX& x, const VecY& y, T s, sparse_tag, sparse_tag)
-{
- if (x.nnz() < y.nnz())
- return sparse_inner_product(x.begin(), x.end(), y.begin(), y.end(), s);
- else
- return sparse_inner_product(y.begin(), y.end(), x.begin(), x.end(), s);
-}
-
-template <class VecX, class VecY, class T>
-inline T
-dot(const VecX& x, const VecY& y, T s, dense_tag, sparse_tag)
-{
- return sparse_dense_inner_product(y.begin(), y.end(), x.begin(), s);
-}
-
-template <class VecX, class VecY, class T>
-inline T
-dot(const VecX& x, const VecY& y, T s, sparse_tag, dense_tag)
-{
- return sparse_dense_inner_product(x.begin(), x.end(), y.begin(), s);
-}
-
-
-//: Dot Product: <tt>s <- x . y + s</tt>
-// The type used for argument s determines the
-// type of the resulting product.
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-template <class VecX, class VecY, class T>
-inline T
-dot(const VecX& x, const VecY& y, T s) MTL_THROW_ASSERTION
-{
- MTL_ASSERT(x.size() == y.size(), "mtl::dot()");
- typedef typename linalg_traits<VecX>::sparsity SparseX;
- typedef typename linalg_traits<VecY>::sparsity SparseY;
- return dot(x, y, s, SparseX(), SparseY());
-}
-
-
-//: Dot Product: <tt>s <- x . y</tt>
-// The type of the resulting product is <TT>VecX::value_type</TT>.
-//!category: algorithms
-//!component: function
-//!example: vecvec_dot.cc, dot_prod.cc
-//!definition: mtl.h
-template <class VecX, class VecY>
-inline typename VecX::value_type
-dot(const VecX& x, const VecY& y) MTL_THROW_ASSERTION
-{
- typedef typename VecX::value_type T;
- return mtl::dot(x, y, T(0));
-}
-
-#ifdef USE_DOUBLE_DOUBLE
-//: Dot Product (extended precision): <tt>s <- x . y + s</tt>
-// The type of the resulting product is double_double
-// Extended precision is used internally.
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-template <class VecX, class VecY>
-inline double_double
-dot(const VecX& x, const VecY& y, double_double s) MTL_THROW_ASSERTION
-{
- typedef typename VecX::value_type x_type;
- typedef typename VecY::value_type y_type;
- // x_type and y_type must be either float or double
- multiply<double_double, x_type, y_type> m;
- addition<double_double, double_double, double_double> a;
- return mtl_algo::inner_product(x.begin(), x.end(), y.begin(), s, a, m);
-}
-#endif /* USE_DOUBLE_DOUBLE */
-
-template <class T>
-struct conj_func {
- typedef T result_type;
- inline T operator()(const T& x) const { return MTL_CONJ(x); }
-};
-
-template <class VecX, class VecY, class T>
-inline T
-dot_conj(const VecX& x, const VecY& y, T s, fast::count<0>)
-{
- return mtl_algo::inner_product(x.begin(), x.end(),
- trans_iter(y.begin(), conj_func<T>()), s);
-}
-#if USE_BLAIS
-template <class VecX, class VecY, class T, int N>
-inline T
-dot_conj(const VecX& x, const VecY& y, T s, fast::count<N>)
-{
- return fast::inner_product(x.begin(), x.end(),
- trans_iter(y.begin(), conj_func<T>()), s);
-}
-#endif
-
-//: Dot Conjugate: <tt>s <- x . conj(y) + s</tt>
-// Similar to dot product. The complex conjugate of the elements of y
-// is used. For real numbers, the conjugate is just that real number.
-// Note that the type of parameter s is the return type of this
-// function.
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-template <class VecX, class VecY, class T>
-inline T
-dot_conj(const VecX& x, const VecY& y, T s) MTL_THROW_ASSERTION
-{
- MTL_ASSERT(x.size() <= y.size(), "mtl::dot_conj()");
- return dot_conj(x, y, s, dim_n<VecX>::RET());
-}
-
-//: Dot Conjugate: <tt>s <- x . conj(y)</tt>
-// A slightly simpler version of the dot conjugate.
-// The return type is the element type of vector x.
-//!category: algorithms
-//!component: function
-//!definition: mtl.h
-template <class VecX, class VecY>
-inline typename VecX::value_type
-dot_conj(const VecX& x, const VecY& y) MTL_THROW_ASSERTION
-{
- typedef typename VecX::value_type T;
- return mtl::dot_conj(x, y, T(0));
-}
-
-
-
-
} /* namespace mtl */
diff -drub mtl-2.1.2-21/mtl/mtl_algo.h mtl-2.1.2-21-icc/mtl/mtl_algo.h
--- mtl-2.1.2-21/mtl/mtl_algo.h Sun Jul 7 01:48:35 2002
+++ mtl-2.1.2-21-icc/mtl/mtl_algo.h Fri Aug 23 19:55:57 2002
@@ -214,7 +214,7 @@
*result = sum2;
++result;
if( ++first==last ) break;
- sum1.~sum_type();
+ sum1.~value_type();
new((void*)&sum1) sum_type( sum2 + *first );
}
}
@@ -237,7 +237,7 @@
*result = sum2;
++result;
if( ++first==last ) break;
- sum1.~sum_type();
+ sum1.~value_type();
new((void*)&sum1) sum_type( op( sum2, *first ) );
}
}
diff -drub mtl-2.1.2-21/mtl/reverse_iter.h mtl-2.1.2-21-icc/mtl/reverse_iter.h
--- mtl-2.1.2-21/mtl/reverse_iter.h Sun Jul 7 01:48:36 2002
+++ mtl-2.1.2-21-icc/mtl/reverse_iter.h Fri Aug 23 19:58:17 2002
@@ -68,15 +68,15 @@
reverse_iter(Iter x) : super(x) {}
inline difference_type index() const {
- Iter tmp = current;
+ Iter tmp = this->current;
return (--tmp).index();
}
inline difference_type row() const {
- Iter tmp = current;
+ Iter tmp = this->current;
return (--tmp).row();
}
inline difference_type column() const {
- Iter tmp = current;
+ Iter tmp = this->current;
return (--tmp).column();
}
diff -drub mtl-2.1.2-21/mtl/transform_iterator.h mtl-2.1.2-21-icc/mtl/transform_iterator.h
--- mtl-2.1.2-21/mtl/transform_iterator.h Sun Jul 7 01:48:36 2002
+++ mtl-2.1.2-21-icc/mtl/transform_iterator.h Fri Aug 23 19:59:50 2002
@@ -264,24 +264,24 @@
inline general_stride_iterator(RandomAccessIterator i, StrideGen s)
: super(*this, i), stride_gen(s) { }
inline self& operator++() {
- stride_gen.inc(iter);
+ stride_gen.inc(this->iter);
return *this;
}
inline self& operator--() {
- stride_gen.dec(iter);
+ stride_gen.dec(this->iter);
return *this;
}
inline self& operator+=(Distance n) {
- stride_gen.advance(iter, n);
+ stride_gen.advance(this->iter, n);
return *this;
}
inline self& operator-=(Distance n) {
- stride_gen.advance(iter, -n);
+ stride_gen.advance(this->iter, -n);
return *this;
}
inline Distance diff(const self& y) const {
- return stride_gen.diff(iter, y.iter, y.stride_gen);
+ return stride_gen.diff(this->iter, y.iter, y.stride_gen);
}
inline friend Distance operator-(const self& x, const self& y) {
diff -drub mtl-2.1.2-21/mtl/utils.h mtl-2.1.2-21-icc/mtl/utils.h
--- mtl-2.1.2-21/mtl/utils.h Sun Jul 7 01:48:36 2002
+++ mtl-2.1.2-21-icc/mtl/utils.h Sat Aug 24 16:17:37 2002
@@ -43,69 +43,6 @@
/* Utility Functions */
-
-template <class Vector>
-inline void
-print_partitioned_vector(Vector x)
-{
- for (typename Vector::iterator i = x.begin();
- i != x.end(); ++i) {
- mtl::print_vector(*i);
- }
-}
-
-template <class Matrix>
-inline void
-print_partitioned_matrix(const Matrix& A)
-{
- std::cout << "Top mat: " << A.nrows() << "x" << A.ncols() << std::endl;
- typedef typename mtl::matrix_traits<Matrix>::size_type Int;
- Int i,j;
- for (i=0; i < A.nrows(); ++i) {
- for (j=0; j < A.ncols(); ++j)
- mtl::print_all_matrix( A(i,j) );
- std::cout << std::endl;
- }
-}
-
-template <class Matrix>
-inline void
-print_partitioned_by_row(const Matrix& A)
-{
- std::cout << "Top mat: " << A.nrows() << "x" << A.ncols() << std::endl;
- typename Matrix::const_iterator A_kk;
- typename Matrix::Row::const_iterator A_i, A_iend;
- A_kk = A.begin();
- while (not_at(A_kk, A.end())) {
- A_i = (*A_kk).begin();
- A_iend = (*A_kk).end();
- while (not_at(A_i, A_iend)) {
- mtl::print_all_matrix( *A_i );
- ++A_i;
- }
- ++A_kk;
- }
-}
-
-template <class Matrix>
-inline void
-print_partitioned_by_column(const Matrix& A)
-{
- std::cout << "Top mat: " << A.nrows() << "x" << A.ncols() << std::endl;
- typename Matrix::const_iterator A_kk;
- typename Matrix::Column::const_iterator A_i, A_iend;
- A_kk = A.begin();
- while (not_at(A_kk, A.end())) {
- A_i = (*A_kk).begin();
- A_iend = (*A_kk).end();
- while (not_at(A_i, A_iend)) {
- mtl::print_all_matrix( *A_i );
- ++A_i;
- }
- ++A_kk;
- }
-}
-
//: utility function for banded
// should use this in a couple more places
//!noindex:
@@ -119,7 +56,6 @@
return dimension<size_type>(start,fini);
}
-
template <class Iterator>
inline void
print_vector(Iterator y, Iterator y_end)
@@ -372,6 +308,69 @@
++A_kk;
}
std::cout << "]" << std::endl;
+}
+
+
+template <class Vector>
+inline void
+print_partitioned_vector(Vector x)
+{
+ for (typename Vector::iterator i = x.begin();
+ i != x.end(); ++i) {
+ mtl::print_vector(*i);
+ }
+}
+
+template <class Matrix>
+inline void
+print_partitioned_matrix(const Matrix& A)
+{
+ std::cout << "Top mat: " << A.nrows() << "x" << A.ncols() << std::endl;
+ typedef typename mtl::matrix_traits<Matrix>::size_type Int;
+ Int i,j;
+ for (i=0; i < A.nrows(); ++i) {
+ for (j=0; j < A.ncols(); ++j)
+ mtl::print_all_matrix( A(i,j) );
+ std::cout << std::endl;
+ }
+}
+
+template <class Matrix>
+inline void
+print_partitioned_by_row(const Matrix& A)
+{
+ std::cout << "Top mat: " << A.nrows() << "x" << A.ncols() << std::endl;
+ typename Matrix::const_iterator A_kk;
+ typename Matrix::Row::const_iterator A_i, A_iend;
+ A_kk = A.begin();
+ while (not_at(A_kk, A.end())) {
+ A_i = (*A_kk).begin();
+ A_iend = (*A_kk).end();
+ while (not_at(A_i, A_iend)) {
+ mtl::print_all_matrix( *A_i );
+ ++A_i;
+ }
+ ++A_kk;
+ }
+}
+
+template <class Matrix>
+inline void
+print_partitioned_by_column(const Matrix& A)
+{
+ std::cout << "Top mat: " << A.nrows() << "x" << A.ncols() << std::endl;
+ typename Matrix::const_iterator A_kk;
+ typename Matrix::Column::const_iterator A_i, A_iend;
+ A_kk = A.begin();
+ while (not_at(A_kk, A.end())) {
+ A_i = (*A_kk).begin();
+ A_iend = (*A_kk).end();
+ while (not_at(A_i, A_iend)) {
+ mtl::print_all_matrix( *A_i );
+ ++A_i;
+ }
+ ++A_kk;
+ }
}
#if 1