Fix Trilinos mat-mat multiplication

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 60152441619a6d212227718846b7ccd905e5bea3..4878966edcb4747665c4cce1d61063b5b2a7c4d1 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -42,7 +42,7 @@ inconvenience this causes.
   <br>
   (Daniel Arndt, 2016/04/18)
   </li>
-  
+
   <li> Changed: FlatManifold takes as argument a periodicity option. This
   used to be a Point<dim>, but it should have been a Tensor<1,dim>. This
   is now changed.
@@ -201,6 +201,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+ <li> Fixed: The methods TrilinosWrappers::SparseMatrix::(T)mmult previously
+ produced invalid matrix sizes if the final matrix was non-square. This has
+ been fixed.
+ <br>
+ (Martin Kronbichler, Daniel Jodlbauer, 2016/04/21)
+ </li>
+
  <li> New: Added an optional string parameter to the ParameterHandler::read_input ()
  and ParameterHandler::read_input_from_string() functions.
  When a line which equals this string is encountered, the parsing of parameters
@@ -209,7 +216,6 @@ inconvenience this causes.
  (Denis Davydov, 2016/04/20)
  </li>
 
-
  <li> New: Added move operations to IndexSet.
  <br>
  (Daniel Shapero, 2016/04/19)

diff --git a/source/lac/trilinos_sparse_matrix.cc b/source/lac/trilinos_sparse_matrix.cc
index ba82179ff5405356d6d68b47b81a8f67ace519e4..56ad9213e02e3c8db170dfb5a657ebf2e27be865 100644 (file)
--- a/source/lac/trilinos_sparse_matrix.cc
+++ b/source/lac/trilinos_sparse_matrix.cc
@@ -2295,7 +2295,10 @@ namespace TrilinosWrappers
       // matrix that we got as a result.
       Epetra_CrsMatrix *C_mat;
       ML_Operator2EpetraCrsMatrix(C_, C_mat);
-      C_mat->FillComplete();
+      C_mat->FillComplete(mod_B->DomainMap(),
+                          transpose_left ?
+                          inputleft.trilinos_matrix().DomainMap() :
+                          inputleft.trilinos_matrix().RangeMap());
       C_mat->OptimizeStorage();
       result.reinit (*C_mat);
 

diff --git a/tests/trilinos/sparse_matrix_mmult_01.cc b/tests/trilinos/sparse_matrix_mmult_01.cc
new file mode 100644 (file)
index 0000000..857f569
--- /dev/null
+++ b/tests/trilinos/sparse_matrix_mmult_01.cc
@@ -0,0 +1,81 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, 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 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check TrilinosWrappers::SparseMatrix::mmult for non-square matrices
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+
+template<typename MATRIX>
+void out_matrix_size(const MATRIX &M,
+                     const std::string &name)
+{
+  deallog << name << " = " << M.m() << " x " << M.n() << std::endl;
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv);
+
+  mpi_initlog();
+
+  TrilinosWrappers::SparsityPattern P_A;
+  TrilinosWrappers::SparseMatrix A;
+  P_A.reinit(4, 3);
+  P_A.compress();
+  A.reinit(P_A);
+
+  TrilinosWrappers::SparsityPattern P_B;
+  TrilinosWrappers::SparseMatrix B;
+  P_B.reinit(3, 5);
+  P_B.compress();
+  B.reinit(P_B);
+
+  out_matrix_size(A, "A");
+  out_matrix_size(B, "B");
+
+  TrilinosWrappers::SparseMatrix AB;
+  A.mmult(AB, B);
+
+  out_matrix_size(AB, "A.B");
+
+  Assert (AB.m() == A.m(), ExcInternalError());
+  Assert (AB.n() == B.n(), ExcInternalError());
+
+  TrilinosWrappers::SparsityPattern P_C;
+  TrilinosWrappers::SparseMatrix C;
+  P_C.reinit(4, 5);
+  P_C.compress();
+  C.reinit(P_C);
+
+  out_matrix_size(C, "C");
+
+  TrilinosWrappers::SparseMatrix AtC;
+  A.Tmmult(AtC, C);
+
+  out_matrix_size(AtC, "At.C");
+
+  Assert (AtC.m() == A.n(), ExcInternalError());
+  Assert (AtC.n() == C.n(), ExcInternalError());
+}

diff --git a/tests/trilinos/sparse_matrix_mmult_01.output b/tests/trilinos/sparse_matrix_mmult_01.output
new file mode 100644 (file)
index 0000000..686d240
--- /dev/null
+++ b/tests/trilinos/sparse_matrix_mmult_01.output
@@ -0,0 +1,6 @@
+
+DEAL::A = 4 x 3
+DEAL::B = 3 x 5
+DEAL::A.B = 4 x 5
+DEAL::C = 4 x 5
+DEAL::At.C = 3 x 5