diff --git a/dune/istl/matrix.hh b/dune/istl/matrix.hh index 946d28ea4c678b52efaeb19bff2ec1c6e80f0bb5..e8f8e156322647610259245c9fe1c565a35bbdfd 100644 --- a/dune/istl/matrix.hh +++ b/dune/istl/matrix.hh @@ -738,9 +738,9 @@ namespace MatrixImp /** \brief Return the transpose of the matrix */ Matrix transpose() const { - Matrix out(N(), M()); - for (size_type i=0; i<M(); i++) - for (size_type j=0; j<N(); j++) + Matrix out(M(), N()); + for (size_type i=0; i<N(); i++) + for (size_type j=0; j<M(); j++) out[j][i] = (*this)[i][j]; return out; diff --git a/dune/istl/test/matrixtest.cc b/dune/istl/test/matrixtest.cc index 21ac0482b21db580d049d2404fa86596cd85ca69..3214063ea5db9ea693d66f418558cff0f4b0f877 100644 --- a/dune/istl/test/matrixtest.cc +++ b/dune/istl/test/matrixtest.cc @@ -280,6 +280,20 @@ void testSolve(const MatrixType& matrix) DUNE_THROW(ISTLError, "Solve() method doesn't appear to produce the solution!"); } +// ////////////////////////////////////////////////////////////// +// Test transposing the matrix +// ////////////////////////////////////////////////////////////// +template <class MatrixType> +void testTranspose(const MatrixType& matrix) +{ + MatrixType transposedMatrix = matrix.transpose(); + + for(size_t i = 0; i < matrix.N(); i++) + for(size_t j = 0; j < matrix.M(); j++) + if(fabs(transposedMatrix[j][i] - matrix[i][j]) > 1e-10) + DUNE_THROW(ISTLError, "transpose() method produces wrong result!"); +} + int main() { @@ -313,6 +327,17 @@ int main() testSuperMatrix(matrix); + Matrix<FieldMatrix<double,1,1> > nonquadraticMatrix(1,2); + { + size_t n = 1; + for (size_t i=0; i<1; i++) + for (size_t j=0; j<2; j++) + nonquadraticMatrix[i][j] = n++; + } + + testTranspose(nonquadraticMatrix); + + // //////////////////////////////////////////////////////////// // Test the BCRSMatrix class -- a sparse dynamic matrix // ////////////////////////////////////////////////////////////