From 09537d22aa55d8d1ef1810a7607240d3c09744c5 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 25 Nov 2008 09:59:11 +0000 Subject: [PATCH] When initializing the Trilinos sparse matrix, use an Epetra_CrsGraph as in intermediate object. That saves some computing time in setup. git-svn-id: https://svn.dealii.org/trunk@17727 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/source/trilinos_sparse_matrix.cc | 127 +++++++++++-------- deal.II/lac/source/trilinos_vector_base.cc | 53 +++++++- 2 files changed, 125 insertions(+), 55 deletions(-) diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index d45ebca374..d8fed0969d 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -256,7 +256,7 @@ namespace TrilinosWrappers { matrix.reset(); - unsigned int n_rows = sparsity_pattern.n_rows(); + const unsigned int n_rows = sparsity_pattern.n_rows(); if (input_row_map.Comm().MyPID() == 0) { @@ -276,21 +276,28 @@ namespace TrilinosWrappers for (unsigned int row=0; row - (new Epetra_FECrsMatrix(Copy, row_map, - &n_entries_per_row[0], false)); + Epetra_CrsGraph graph (Copy, row_map, + &n_entries_per_row[input_row_map.MinMyGID()], true); // TODO: As of now, assume that the @@ -299,10 +306,10 @@ namespace TrilinosWrappers // each processor set its rows. Since // this is wasteful, a better solution // should be found in the future. - Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), - ExcDimensionMismatch (matrix->NumGlobalRows(), - sparsity_pattern.n_rows())); - + Assert (graph.NumGlobalRows() == (int)sparsity_pattern.n_rows(), + ExcDimensionMismatch (graph.NumGlobalRows(), + sparsity_pattern.n_rows())); + // Trilinos has a bug for rectangular // matrices at this point, so do not // check for consistent column numbers @@ -311,26 +318,32 @@ namespace TrilinosWrappers // this bug is filed in the Sandia // bugzilla under #4123 and should be // fixed for version 9.0 -// Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(), -// ExcDimensionMismatch (matrix->NumGlobalCols(), +// Assert (graph.NumGlobalCols() == (int)sparsity_pattern.n_cols(), +// ExcDimensionMismatch (graph.NumGlobalCols(), // sparsity_pattern.n_cols())); - std::vector values; - std::vector row_indices; + std::vector row_indices; for (unsigned int row=0; row + (new Epetra_FECrsMatrix(Copy, graph, false)); + last_action = Zero; // In the end, the matrix needs to @@ -354,7 +367,7 @@ namespace TrilinosWrappers { matrix.reset(); - unsigned int n_rows = sparsity_pattern.n_rows(); + const unsigned int n_rows = sparsity_pattern.n_rows(); if (input_row_map.Comm().MyPID() == 0) { @@ -374,21 +387,29 @@ namespace TrilinosWrappers for (unsigned int row=0; row - (new Epetra_FECrsMatrix(Copy, row_map, - &n_entries_per_row[0], false)); + Epetra_CrsGraph graph (Copy, row_map, + &n_entries_per_row[input_row_map.MinMyGID()], true); // TODO: As of now, assume that the @@ -397,10 +418,10 @@ namespace TrilinosWrappers // each processor set its rows. Since // this is wasteful, a better solution // should be found in the future. - Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), - ExcDimensionMismatch (matrix->NumGlobalRows(), - sparsity_pattern.n_rows())); - + Assert (graph.NumGlobalRows() == (int)sparsity_pattern.n_rows(), + ExcDimensionMismatch (graph.NumGlobalRows(), + sparsity_pattern.n_rows())); + // Trilinos has a bug for rectangular // matrices at this point, so do not // check for consistent column numbers @@ -409,19 +430,17 @@ namespace TrilinosWrappers // this bug is filed in the Sandia // bugzilla under #4123 and should be // fixed for version 9.0 -// Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(), -// ExcDimensionMismatch (matrix->NumGlobalCols(), +// Assert (graph.NumGlobalCols() == (int)sparsity_pattern.n_cols(), +// ExcDimensionMismatch (graph.NumGlobalCols(), // sparsity_pattern.n_cols())); - std::vector values; - std::vector row_indices; + std::vector row_indices; for (unsigned int row=0; row + (new Epetra_FECrsMatrix(Copy, graph, false)); last_action = Zero; diff --git a/deal.II/lac/source/trilinos_vector_base.cc b/deal.II/lac/source/trilinos_vector_base.cc index 37ebf97a76..6147ec10ed 100644 --- a/deal.II/lac/source/trilinos_vector_base.cc +++ b/deal.II/lac/source/trilinos_vector_base.cc @@ -315,7 +315,7 @@ namespace TrilinosWrappers TrilinosScalar VectorBase::operator * (const VectorBase &vec) const { - Assert (size() == vec.size(), + Assert (local_range() == vec.local_range(), ExcDimensionMismatch(size(), vec.size())); TrilinosScalar result; @@ -766,9 +766,31 @@ namespace TrilinosWrappers ExcMessage("The given value is not finite but " "either infinite or Not A Number (NaN)")); - *vector = *v.vector; + // If we don't have the same map, copy. + if (local_range() != v.local_range()) + { + vector.reset(); + last_action = Zero; + *vector = *v.vector; + *this *= a; + } + else + { + // Otherwise, just update + Assert (vector->Map().SameAs(v.vector->Map()) == true, + ExcMessage ("The Epetra maps in the assignment operator =" + " do not match, even though the local_range " + " seems to be the same. Check vector setup!")); + int ierr; + ierr = vector->GlobalAssemble(last_action); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = vector->Update(a, *v.vector, 0.0); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + last_action = Zero; + } - *this *= a; } @@ -790,9 +812,30 @@ namespace TrilinosWrappers ExcMessage("The given value is not finite but " "either infinite or Not A Number (NaN)")); - *vector = *v.vector; + // If we don't have the same map, copy. + if (local_range() != v.local_range()) + { + vector.reset(); + last_action = Zero; + *vector = *v.vector; + sadd (a, b, w); + } + else + { + // Otherwise, just update + Assert (vector->Map().SameAs(v.vector->Map()) == true, + ExcMessage ("The Epetra maps in the assignment operator =" + " do not match, even though the local_range " + " seems to be the same. Check vector setup!")); + int ierr; + ierr = vector->GlobalAssemble(last_action); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - sadd (a, b, w); + ierr = vector->Update(a, *v.vector, b, *w.vector, 0.0); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + last_action = Zero; + } } -- 2.39.5