Hi,
I figured out what was happening, by liberally sprinking "cout" into the MTL
source code. I had actually had a bug in a different routine, but I als found
out that I don't need the assignment statement before the copy. My code now
looks like:
gengen(const int n, vector< Array< Matrix_t, 1 > >& generators)
...
// Use recursion to define generator array in terms of smaller generator array
Array< Matrix_t, 1 >
result(gen(gengen<Matrix_t>(n-1, generators)));
// Save the resulting generator array.
// This is tedious because Array assignment does not copy shape.
generators.resize(n+1,result);
generators[n].resize(result.shape());
generators[n] = result;
for( int k = result.lbound(firstDim); k <= result.ubound(firstDim); ++k)
{
mtl::copy( result(k), generators[n](k));
}
The code which had the bug was:
kron(const Matrix_t& x, const Matrix_t& y, Matrix_t& z)
{
const int rx(x.nrows()); const int cx(x.ncols());
const int ry(y.nrows()); const int cy(y.ncols());
Matrix_t result(rx*ry, cx*cy);
z = result;
mtl::set(z, 0);
for (typename Matrix_t::const_iterator i = x.begin(); i != x.end(); ++i)
for (typename Matrix_t::Row::const_iterator j = (*i).begin();
j != (*i).end(); ++j)
for (typename Matrix_t::const_iterator k = y.begin(); k != y.end(); ++k)
for (typename Matrix_t::Row::const_iterator l = (*k).begin();
l != (*k).end(); ++l)
z(j.row()*ry + l.row(), j.column()*cy + l.column()) = (*j) * (*l);
}
I had put:
const int ry(x.nrows()); const int cy(x.ncols());
and this probably caused random-seeming errors, as memory was overwritten.
The moral is, that if you have a template library, you have the source, and
even if it seems complicated, you can always insert debug print statements to
see what it is doing.
--------------------- Object Oriented Numerics List --------------------------
* To subscribe/unsubscribe: use the handy web form at
http://oonumerics.org/oon/
* If this doesn't work, please send a note to owner-oon-list@oonumerics.org
This archive was generated by hypermail 2b29 : Wed Feb 20 2002 - 03:20:16 EST