From 6e985c6dc28ef3a695765071e1699ab404ea6959 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 25 Feb 2009 16:16:12 +0000 Subject: [PATCH] Make the Epetra objects in the TrilinosWrapper classes private. We can still handle the situation with the help of some access functions. Hopefully, I didn't forget to set any plain usage of the Trilinos objects... git-svn-id: https://svn.dealii.org/trunk@18430 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/source/grid/grid_refinement.cc | 4 +- deal.II/examples/step-32/step-32.cc | 6 +- deal.II/examples/step-33/step-33.cc | 14 +- .../lac/include/lac/trilinos_sparse_matrix.h | 120 +++++++++++++++- .../include/lac/trilinos_sparsity_pattern.h | 129 +++++++++++++++--- .../lac/include/lac/trilinos_vector_base.h | 51 ++++++- deal.II/lac/include/lac/vector.templates.h | 6 +- deal.II/lac/source/trilinos_precondition.cc | 58 +++++--- .../lac/source/trilinos_precondition_block.cc | 9 +- deal.II/lac/source/trilinos_solver.cc | 30 ++-- deal.II/lac/source/trilinos_solver_block.cc | 27 ++-- deal.II/lac/source/trilinos_sparse_matrix.cc | 41 +++--- deal.II/lac/source/trilinos_vector.cc | 6 +- 13 files changed, 391 insertions(+), 110 deletions(-) diff --git a/deal.II/deal.II/source/grid/grid_refinement.cc b/deal.II/deal.II/source/grid/grid_refinement.cc index 8633c1ea41..ba33ff56a1 100644 --- a/deal.II/deal.II/source/grid/grid_refinement.cc +++ b/deal.II/deal.II/source/grid/grid_refinement.cc @@ -83,7 +83,7 @@ namespace max_element (const TrilinosWrappers::Vector &criteria) { TrilinosScalar m = 0; - criteria.vector->MaxValue(&m); + criteria.trilinos_vector().MaxValue(&m); return m; } @@ -92,7 +92,7 @@ namespace min_element (const TrilinosWrappers::Vector &criteria) { TrilinosScalar m = 0; - criteria.vector->MinValue(&m); + criteria.trilinos_vector().MinValue(&m); return m; } #endif diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 867b43f086..3687414cca 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -263,7 +263,7 @@ namespace LinearSolvers stokes_matrix (&S), m_inverse (&Mpinv), a_preconditioner (Apreconditioner), - tmp (stokes_matrix->block(1,1).matrix->RowMap()) + tmp (stokes_matrix->block(1,1).row_partitioner()) {} @@ -1298,8 +1298,8 @@ void BoussinesqFlowProblem::project_temperature_field () std::vector rhs_values(n_q_points); - TrilinosWrappers::MPI::Vector rhs (temperature_mass_matrix.matrix->RowMap()), - sol (temperature_mass_matrix.matrix->RowMap()); + TrilinosWrappers::MPI::Vector rhs (temperature_mass_matrix.row_partitioner()), + sol (temperature_mass_matrix.row_partitioner()); for (; cell!=endc; ++cell) if (cell->subdomain_id() == diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index b6ced66928..46bfe39dd1 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -2671,13 +2671,16 @@ ConservationLaw::solve (Vector &newton_update) // deal.II wrapper class itself. So // we access to the actual Trilinos // matrix in the Trilinos wrapper - // class, and create a plain pointer - // when passing it in. + // class by the command + // trilinos_matrix(). Trilinos wants + // the matrix to be non-constant, so + // we have to manually remove the + // constantness using a const_cast. case Parameters::Solver::gmres: { - Epetra_Vector x(View, system_matrix.matrix->RowMap(), + Epetra_Vector x(View, system_matrix.domain_partitioner(), newton_update.begin()); - Epetra_Vector b(View, system_matrix.matrix->RowMap(), + Epetra_Vector b(View, system_matrix.range_partitioner(), right_hand_side.begin()); AztecOO solver; @@ -2702,7 +2705,8 @@ ConservationLaw::solve (Vector &newton_update) solver.SetAztecParam(AZ_athresh, parameters.ilut_atol); solver.SetAztecParam(AZ_rthresh, parameters.ilut_rtol); - solver.SetUserMatrix(&*system_matrix.matrix); + solver.SetUserMatrix(const_cast + (&system_matrix.trilinos_matrix())); solver.Iterate(parameters.max_iterations, parameters.linear_residual); diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 922aa01c29..4e4285e835 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -1583,6 +1583,69 @@ namespace TrilinosWrappers TrilinosScalar frobenius_norm () const; //@} +/** + * @name Access to underlying Trilinos data + */ +//@{ + + /** + * Return a const reference to the + * underlying Trilinos + * Epetra_CrsMatrix data. + */ + const Epetra_CrsMatrix & trilinos_matrix () const; + + /** + * Return a const reference to the + * underlying Trilinos + * Epetra_CrsGraph data that stores + * the sparsity pattern of the + * matrix. + */ + const Epetra_CrsGraph & trilinos_sparsity_pattern () const; + + /** + * Return a const reference to the + * underlying Trilinos Epetra_Map + * that sets the partitioning of the + * domain space of this matrix, i.e., + * the partitioning of the vectors + * this matrix has to be multiplied + * with. + */ + const Epetra_Map & domain_partitioner () const; + + /** + * Return a const reference to the + * underlying Trilinos Epetra_Map + * that sets the partitioning of the + * range space of this matrix, i.e., + * the partitioning of the vectors + * that are result from matrix-vector + * products. + */ + const Epetra_Map & range_partitioner () const; + + /** + * Return a const reference to the + * underlying Trilinos Epetra_Map + * that sets the partitioning of the + * matrix rows. Equal to the + * partitioning of the range. + */ + const Epetra_Map & row_partitioner () const; + + /** + * Return a const reference to the + * underlying Trilinos Epetra_Map + * that sets the partitioning of the + * matrix columns. This is in general + * not equal to the partitioner + * Epetra_Map for the domain because + * of overlap in the matrix. + */ + const Epetra_Map & col_partitioner () const; +//@} /** * @name Iterators */ @@ -1793,7 +1856,6 @@ namespace TrilinosWrappers */ std::vector column_values; - public: /** * A sparse matrix object in * Trilinos to be used for @@ -1803,9 +1865,6 @@ namespace TrilinosWrappers * elements. The actual type, * a sparse matrix, is set in * the constructor. - * - * TODO: This object should - * finally become private. */ std::auto_ptr matrix; @@ -2455,6 +2514,59 @@ namespace TrilinosWrappers } + + inline + const Epetra_CrsMatrix & + SparseMatrix::trilinos_matrix () const + { + return static_cast(*matrix); + } + + + + inline + const Epetra_CrsGraph & + SparseMatrix::trilinos_sparsity_pattern () const + { + return matrix->Graph(); + } + + + + inline + const Epetra_Map & + SparseMatrix::domain_partitioner () const + { + return matrix->DomainMap(); + } + + + + inline + const Epetra_Map & + SparseMatrix::range_partitioner () const + { + return matrix->RangeMap(); + } + + + + inline + const Epetra_Map & + SparseMatrix::row_partitioner () const + { + return matrix->RowMap(); + } + + + + inline + const Epetra_Map & + SparseMatrix::col_partitioner () const + { + return matrix->ColMap(); + } + #endif // DOXYGEN } diff --git a/deal.II/lac/include/lac/trilinos_sparsity_pattern.h b/deal.II/lac/include/lac/trilinos_sparsity_pattern.h index 461550176b..bf4e2ca5f5 100755 --- a/deal.II/lac/include/lac/trilinos_sparsity_pattern.h +++ b/deal.II/lac/include/lac/trilinos_sparsity_pattern.h @@ -811,6 +811,63 @@ namespace TrilinosWrappers const unsigned int n_cols, const unsigned int *col_indices); //@} +/** + * @name Access of underlying Trilinos data + */ +//@{ + + /** + * Return a const reference to the + * underlying Trilinos + * Epetra_CrsGraph data that stores + * the sparsity pattern. + */ + const Epetra_CrsGraph & trilinos_sparsity_pattern () const; + + /** + * Return a const reference to the + * underlying Trilinos Epetra_Map + * that sets the parallel + * partitioning of the domain space + * of this sparsity pattern, i.e., + * the partitioning of the vectors + * matrices based on this sparsity + * pattern are multiplied with. + */ + const Epetra_Map & domain_partitioner () const; + + /** + * Return a const reference to the + * underlying Trilinos Epetra_Map + * that sets the partitioning of the + * range space of this sparsity + * pattern, i.e., the partitioning of + * the vectors that are result from + * matrix-vector products. + */ + const Epetra_Map & range_partitioner () const; + + /** + * Return a const reference to the + * underlying Trilinos Epetra_Map + * that sets the partitioning of the + * sparsity pattern rows. Equal to + * the partitioning of the range. + */ + const Epetra_Map & row_partitioner () const; + + /** + * Return a const reference to the + * underlying Trilinos Epetra_Map + * that sets the partitioning of the + * sparsity pattern columns. This is + * in general not equal to the + * partitioner Epetra_Map for the + * domain because of overlap in the + * matrix. + */ + const Epetra_Map & col_partitioner () const; +//@} /** * @name Iterators */ @@ -1021,7 +1078,6 @@ namespace TrilinosWrappers */ unsigned int row_in_cache; - friend class SparseMatrix; friend class SparsityPatternIterators::const_iterator; friend class SparsityPatternIterators::const_iterator::Accessor; }; @@ -1309,7 +1365,7 @@ namespace TrilinosWrappers if (n_cols == 0) return; - int * col_index_ptr = (int*)col_indices; + int * col_index_ptr = (int*)(col_indices); compressed = false; int ierr; @@ -1320,25 +1376,17 @@ namespace TrilinosWrappers // Epetra_CrsGraph input function, which // is much faster than the // Epetra_FECrsGraph function. - if (row_map.MyGID(row) == true) - ierr = graph->Epetra_CrsGraph::InsertGlobalIndices(row, - n_cols, - col_index_ptr); - else + //if (row_map.MyGID(row) == true) + //ierr = graph->Epetra_CrsGraph::InsertGlobalIndices(row, + // n_cols, + // col_index_ptr); + //else { // When we're at off-processor data, we // have to stick with the standard - // SumIntoGlobalValues - // function. Nevertheless, the way we - // call it is the fastest one (any other - // will lead to repeated allocation and - // deallocation of memory in order to - // call the function we already use, - // which is very unefficient if writing - // one element at a time). - - ierr = graph->InsertGlobalIndices (1, (int*)&row, n_cols, - col_index_ptr); + // SumIntoGlobalValues method. + + ierr = graph->InsertGlobalIndices (1, (int*)&row, n_cols, col_index_ptr); } //Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0])); @@ -1346,6 +1394,51 @@ namespace TrilinosWrappers } + + inline + const Epetra_CrsGraph & + SparsityPattern::trilinos_sparsity_pattern () const + { + return static_cast(*graph); + } + + + + inline + const Epetra_Map & + SparsityPattern::domain_partitioner () const + { + return (Epetra_Map &) graph->DomainMap(); + } + + + + inline + const Epetra_Map & + SparsityPattern::range_partitioner () const + { + return (Epetra_Map &) graph->RangeMap(); + } + + + + inline + const Epetra_Map & + SparsityPattern::row_partitioner () const + { + return (Epetra_Map &) graph->RowMap(); + } + + + + inline + const Epetra_Map & + SparsityPattern::col_partitioner () const + { + return (Epetra_Map &) graph->ColMap(); + } + + #endif // DOXYGEN } diff --git a/deal.II/lac/include/lac/trilinos_vector_base.h b/deal.II/lac/include/lac/trilinos_vector_base.h index 268f6e2cc8..a9890a017f 100644 --- a/deal.II/lac/include/lac/trilinos_vector_base.h +++ b/deal.II/lac/include/lac/trilinos_vector_base.h @@ -798,6 +798,28 @@ namespace TrilinosWrappers */ //@{ + /** + * Return a const reference to the + * underlying Trilinos + * Epetra_MultiVector class. + */ + const Epetra_MultiVector & trilinos_vector () const; + + /** + * Return a (modifyable) reference to + * the underlying Trilinos + * Epetra_FEVector class. + */ + Epetra_FEVector & trilinos_vector (); + + /** + * Return a const reference to the + * underlying Trilinos Epetra_Map + * that sets the parallel + * partitioning of the vector. + */ + const Epetra_Map & vector_partitioner () const; + /** * Output of vector in * user-defined format in analogy @@ -914,13 +936,10 @@ namespace TrilinosWrappers */ bool compressed; - public: /** * An Epetra distibuted vector * type. Requires an existing * Epetra_Map for storing data. - * TODO: Should become private - * at some point. */ std::auto_ptr vector; @@ -1191,6 +1210,32 @@ namespace TrilinosWrappers } + + inline + const Epetra_MultiVector & + VectorBase::trilinos_vector () const + { + return static_cast(*vector); + } + + + + inline + Epetra_FEVector & + VectorBase::trilinos_vector () + { + return *vector; + } + + + + inline + const Epetra_Map & + VectorBase::vector_partitioner () const + { + return (Epetra_Map &)vector->Map(); + } + #endif // DOXYGEN } diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index b3f6c2515e..c6aad63827 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -196,7 +196,7 @@ Vector::Vector (const TrilinosWrappers::MPI::Vector &v) // and copy it TrilinosScalar **start_ptr; - int ierr = localized_vector.vector->ExtractView (&start_ptr); + int ierr = localized_vector.trilinos_vector().ExtractView (&start_ptr); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); std::copy (start_ptr[0], start_ptr[0]+vec_size, begin()); @@ -222,7 +222,7 @@ Vector::Vector (const TrilinosWrappers::Vector &v) // and copy it TrilinosScalar **start_ptr; - int ierr = v.vector->ExtractView (&start_ptr); + int ierr = v.trilinos_vector().ExtractView (&start_ptr); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); std::copy (start_ptr[0], start_ptr[0]+vec_size, begin()); @@ -723,7 +723,7 @@ Vector::operator = (const TrilinosWrappers::Vector &v) // get a representation of the vector // and copy it TrilinosScalar **start_ptr; - int ierr = v.vector->ExtractView (&start_ptr); + int ierr = v.trilinos_vector().ExtractView (&start_ptr); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); std::copy (start_ptr[0], start_ptr[0]+vec_size, begin()); diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index 6198f20bcc..5019b3fc8c 100755 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -65,12 +65,13 @@ namespace TrilinosWrappers PreconditionBase::vmult (VectorBase &dst, const VectorBase &src) const { - Assert (dst.vector->Map().SameAs(preconditioner->OperatorRangeMap()), + Assert (dst.vector_partitioner().SameAs(preconditioner->OperatorRangeMap()), ExcNonMatchingMaps("dst")); - Assert (src.vector->Map().SameAs(preconditioner->OperatorDomainMap()), + Assert (src.vector_partitioner().SameAs(preconditioner->OperatorDomainMap()), ExcNonMatchingMaps("src")); - const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector); + const int ierr = preconditioner->ApplyInverse (src.trilinos_vector(), + dst.trilinos_vector()); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -133,8 +134,10 @@ namespace TrilinosWrappers preconditioner.release(); ifpack.release(); - ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation", - &*matrix.matrix, 0)); + ifpack = Teuchos::rcp (Ifpack().Create + ("point relaxation", + const_cast(&matrix.trilinos_matrix()), + 0)); Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this " "preconditioner")); @@ -183,9 +186,10 @@ namespace TrilinosWrappers preconditioner.release(); ifpack.release(); - ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation", - &*matrix.matrix, - additional_data.overlap)); + ifpack = Teuchos::rcp (Ifpack().Create + ("point relaxation", + const_cast(&matrix.trilinos_matrix()), + additional_data.overlap)); Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this " "preconditioner")); @@ -235,9 +239,10 @@ namespace TrilinosWrappers preconditioner.release(); ifpack.release(); - ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation", - &*matrix.matrix, - additional_data.overlap)); + ifpack = Teuchos::rcp (Ifpack().Create + ("point relaxation", + const_cast(&matrix.trilinos_matrix()), + additional_data.overlap)); Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this " "preconditioner")); @@ -289,8 +294,10 @@ namespace TrilinosWrappers preconditioner.release(); ifpack.release(); - ifpack = Teuchos::rcp (Ifpack().Create ("IC", &*matrix.matrix, - additional_data.overlap)); + ifpack = Teuchos::rcp (Ifpack().Create + ("IC", + const_cast(&matrix.trilinos_matrix()), + additional_data.overlap)); Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this " "preconditioner")); @@ -340,8 +347,10 @@ namespace TrilinosWrappers preconditioner.release(); ifpack.release(); - ifpack = Teuchos::rcp (Ifpack().Create ("ILU", &*matrix.matrix, - additional_data.overlap)); + ifpack = Teuchos::rcp (Ifpack().Create + ("ILU", + const_cast(&matrix.trilinos_matrix()), + additional_data.overlap)); Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this " "preconditioner")); @@ -393,8 +402,10 @@ namespace TrilinosWrappers preconditioner.release(); ifpack.release(); - ifpack = Teuchos::rcp (Ifpack().Create ("ILUT", &*matrix.matrix, - additional_data.overlap)); + ifpack = Teuchos::rcp (Ifpack().Create + ("ILUT", + const_cast(&matrix.trilinos_matrix()), + additional_data.overlap)); Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this " "preconditioner")); @@ -439,8 +450,10 @@ namespace TrilinosWrappers preconditioner.release(); ifpack.release(); - ifpack = Teuchos::rcp (Ifpack().Create ("Amesos", &*matrix.matrix, - additional_data.overlap)); + ifpack = Teuchos::rcp (Ifpack().Create + ("Amesos", + const_cast(&matrix.trilinos_matrix()), + additional_data.overlap)); Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this " "preconditioner")); @@ -490,7 +503,7 @@ namespace TrilinosWrappers preconditioner.release(); ifpack.release(); - ifpack = Teuchos::rcp (new Ifpack_Chebyshev(&*matrix.matrix)); + ifpack = Teuchos::rcp (new Ifpack_Chebyshev (&matrix.trilinos_matrix())); Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this " "preconditioner")); @@ -610,7 +623,7 @@ namespace TrilinosWrappers else parameter_list.set("ML output", 0); - const Epetra_Map & domain_map = matrix.matrix->DomainMap(); + const Epetra_Map & domain_map = matrix.domain_partitioner(); Epetra_MultiVector distributed_constant_modes (domain_map, constant_modes_dimension); @@ -649,7 +662,8 @@ namespace TrilinosWrappers } multilevel_operator = Teuchos::rcp (new ML_Epetra::MultiLevelPreconditioner( - *matrix.matrix, parameter_list, true)); + matrix.trilinos_matrix(), parameter_list, + true)); if (additional_data.output_details) multilevel_operator->PrintUnused(0); diff --git a/deal.II/lac/source/trilinos_precondition_block.cc b/deal.II/lac/source/trilinos_precondition_block.cc index b2e6b0cec2..e313aa660d 100755 --- a/deal.II/lac/source/trilinos_precondition_block.cc +++ b/deal.II/lac/source/trilinos_precondition_block.cc @@ -61,7 +61,7 @@ namespace TrilinosWrappers // (Thyra wrapper around the // matrix). Teuchos::RCP > tmpM = - Thyra::epetraLinearOp(Teuchos::rcp(&*input_M.matrix, false)); + Thyra::epetraLinearOp(Teuchos::rcp(&input_M.trilinos_matrix(), false)); Teuchos::RCP > M = Teuchos::rcp(new Thyra::DefaultLinearOpSource(tmpM), true); @@ -333,7 +333,8 @@ namespace TrilinosWrappers else { Teuchos::RCP > - tmpFp = Thyra::epetraLinearOp(Teuchos::rcp(&*Fp_matrix.matrix,false)); + tmpFp = Thyra::epetraLinearOp(Teuchos::rcp(&Fp_matrix.trilinos_matrix(), + false)); Fp_op = tmpFp; } @@ -366,12 +367,12 @@ namespace TrilinosWrappers // gradient and divergence // operators in Thyra format. Teuchos::RCP > tmpS01 = - Thyra::epetraLinearOp(Teuchos::rcp(&*system_matrix.block(0,1).matrix, + Thyra::epetraLinearOp(Teuchos::rcp(&system_matrix.block(0,1).trilinos_matrix(), false)); Thyra::ConstLinearOperator S01 = tmpS01; Teuchos::RCP > tmpS10 = - Thyra::epetraLinearOp(Teuchos::rcp(&*system_matrix.block(1,0).matrix, + Thyra::epetraLinearOp(Teuchos::rcp(&system_matrix.block(1,0).trilinos_matrix(), false)); Thyra::ConstLinearOperator S10 = tmpS10; diff --git a/deal.II/lac/source/trilinos_solver.cc b/deal.II/lac/source/trilinos_solver.cc index d22462b4b7..81097f8f10 100644 --- a/deal.II/lac/source/trilinos_solver.cc +++ b/deal.II/lac/source/trilinos_solver.cc @@ -81,9 +81,10 @@ namespace TrilinosWrappers // to let the AztecOO solver // know about the matrix and // vectors. - linear_problem = std::auto_ptr ( - new Epetra_LinearProblem(&*(A.matrix), &*x.vector, - &*b.vector)); + linear_problem = std::auto_ptr + (new Epetra_LinearProblem(const_cast(&A.trilinos_matrix()), + &x.trilinos_vector(), + const_cast(&b.trilinos_vector()))); // Next we can allocate the // AztecOO solver... @@ -164,11 +165,11 @@ namespace TrilinosWrappers ExcDimensionMismatch(b.size(), A.m())); Assert (A.local_range ().second == A.m(), ExcMessage ("Can only work in serial when using deal.II vectors.")); - Assert (A.matrix->Filled(), + Assert (A.trilinos_matrix().Filled(), ExcMessage ("Matrix is not compressed. Call compress() method.")); - Epetra_Vector ep_x (View, A.matrix->DomainMap(), x.begin()); - Epetra_Vector ep_b (View, A.matrix->RangeMap(), const_cast(b.begin())); + Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin()); + Epetra_Vector ep_b (View, A.range_partitioner(), const_cast(b.begin())); // We need an // Epetra_LinearProblem object @@ -176,7 +177,8 @@ namespace TrilinosWrappers // know about the matrix and // vectors. linear_problem = std::auto_ptr - (new Epetra_LinearProblem(&*(A.matrix), &ep_x, &ep_b)); + (new Epetra_LinearProblem + (const_cast(&A.trilinos_matrix()), &ep_x, &ep_b)); // Next we can allocate the // AztecOO solver... @@ -398,9 +400,10 @@ namespace TrilinosWrappers // to let the AztecOO solver // know about the matrix and // vectors. - linear_problem = std::auto_ptr ( - new Epetra_LinearProblem(&*(A.matrix), &*x.vector, - &*b.vector)); + linear_problem = std::auto_ptr + (new Epetra_LinearProblem(const_cast(&A.trilinos_matrix()), + &x.trilinos_vector(), + const_cast(&b.trilinos_vector()))); // Next we can allocate the // AztecOO solver... @@ -456,8 +459,8 @@ namespace TrilinosWrappers ExcDimensionMismatch(b.size(), A.m())); Assert (A.local_range ().second == A.m(), ExcMessage ("Can only work in serial when using deal.II vectors.")); - Epetra_Vector ep_x (View, A.matrix->DomainMap(), x.begin()); - Epetra_Vector ep_b (View, A.matrix->RangeMap(), const_cast(b.begin())); + Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin()); + Epetra_Vector ep_b (View, A.range_partitioner(), const_cast(b.begin())); // We need an // Epetra_LinearProblem object @@ -465,7 +468,8 @@ namespace TrilinosWrappers // know about the matrix and // vectors. linear_problem = std::auto_ptr - (new Epetra_LinearProblem(&*(A.matrix), &ep_x, &ep_b)); + (new Epetra_LinearProblem + (const_cast(&A.trilinos_matrix()), &ep_x, &ep_b)); // Next we can allocate the // AztecOO solver... diff --git a/deal.II/lac/source/trilinos_solver_block.cc b/deal.II/lac/source/trilinos_solver_block.cc index 63b374fa74..a9f1ed8faa 100644 --- a/deal.II/lac/source/trilinos_solver_block.cc +++ b/deal.II/lac/source/trilinos_solver_block.cc @@ -104,9 +104,9 @@ namespace TrilinosWrappers MPI::BlockVector tmpb (b.n_blocks()); for (unsigned int i=0; iDomainMap(), false); + tmpx.block(i).reinit (A.block(i,i).domain_partitioner(), false); tmpx.block(i) = x.block(i); - tmpb.block(i).reinit (A.block(i,i).matrix->RangeMap(), false); + tmpb.block(i).reinit (A.block(i,i).range_partitioner(), false); tmpb.block(i) = b.block(i); } tmpx.collect_sizes(); @@ -153,9 +153,9 @@ namespace TrilinosWrappers for (unsigned int i=0; iMap().UniqueGIDs() == true, + AssertThrow (input_x.block(i).vector_partitioner().UniqueGIDs() == true, ExcOverlappingMaps("vector", "x")); - AssertThrow (input_b.block(i).vector->Map().UniqueGIDs() == true, + AssertThrow (input_b.block(i).vector_partitioner().UniqueGIDs() == true, ExcOverlappingMaps("vector", "b")); for (unsigned int j=0; jMap().SameAs( - input_A.block(i,j).matrix->DomainMap()) == true, + AssertThrow (input_x.block(j).vector_partitioner().SameAs( + input_A.block(i,j).domain_partitioner()) == true, ExcNonMatchingMaps (error_component.str())); } { @@ -173,8 +173,8 @@ namespace TrilinosWrappers i_str << j; std::ostringstream error_component; error_component << "b.block(" << i_str.str() << ")"; - AssertThrow (input_b.block(i).vector->Map().SameAs( - input_A.block(i,j).matrix->RangeMap()) == true, + AssertThrow (input_b.block(i).vector_partitioner().SameAs( + input_A.block(i,j).range_partitioner()) == true, ExcNonMatchingMaps (error_component.str())); } } @@ -196,8 +196,8 @@ namespace TrilinosWrappers { Teuchos::RCP > tmp_space = Thyra::create_VectorSpace( - Teuchos::rcp(&input_A.block(i,i).matrix->DomainMap(), - false)); + Teuchos::rcp(&input_A.block(i,i).domain_partitioner(), + false)); epetra_vector_spaces.push_back(tmp_space); } @@ -212,7 +212,7 @@ namespace TrilinosWrappers for (unsigned int j=0; j > - A_ij = Thyra::epetraLinearOp(Teuchos::rcp(&*input_A.block(i,j).matrix, + A_ij = Thyra::epetraLinearOp(Teuchos::rcp(&input_A.block(i,j).trilinos_matrix(), false)); tmpA->setBlock(i, j, A_ij); } @@ -229,14 +229,15 @@ namespace TrilinosWrappers Thyra::Vector sol = A.domain().createMember(); for (unsigned int i=0; i + ((input_b.block(i).trilinos_vector())(0)); Teuchos::RCP > tmp_rhs_i = Thyra::create_Vector(Teuchos::rcp(block_ptr, false), epetra_vector_spaces[i]); const Thyra::Vector rhs_i = tmp_rhs_i; rhs.setBlock(i, rhs_i); - block_ptr = (*input_x.block(i).vector)(0); + block_ptr = (input_x.block(i).trilinos_vector())(0); Teuchos::RCP > tmp_sol_i = Thyra::create_Vector(Teuchos::rcp(block_ptr, false), epetra_vector_spaces[i]); diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index c22d34781d..bdd7c77a2b 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -53,8 +53,9 @@ namespace TrilinosWrappers TrilinosScalar *values = new TrilinosScalar(colnums); int ierr; - ierr = matrix->matrix->ExtractGlobalRowCopy((int)this->a_row, colnums, - ncols, &(values[0])); + ierr = matrix->trilinos_matrix().ExtractGlobalRowCopy((int)this->a_row, + colnums, + ncols, &(values[0])); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); // copy it into our caches if the @@ -213,14 +214,16 @@ namespace TrilinosWrappers SparseMatrix::SparseMatrix (const SparsityPattern &InputSP) : Subscriptor(), - row_map (InputSP.row_map), - col_map (InputSP.col_map), + row_map (InputSP.range_partitioner()), + col_map (InputSP.domain_partitioner()), last_action (Zero), compressed (true), matrix (std::auto_ptr - (new Epetra_FECrsMatrix(Copy, *InputSP.graph, false))) + (new Epetra_FECrsMatrix(Copy, + InputSP.trilinos_sparsity_pattern(), + false))) { - Assert(InputSP.graph->Filled() == true, + Assert(InputSP.trilinos_sparsity_pattern().Filled() == true, ExcMessage("The Trilinos sparsity pattern has not been compressed.")); compress(); } @@ -543,14 +546,15 @@ namespace TrilinosWrappers { matrix.reset(); - row_map = sparsity_pattern.row_map; - col_map = sparsity_pattern.col_map; + row_map = sparsity_pattern.range_partitioner(); + col_map = sparsity_pattern.domain_partitioner(); - Assert (sparsity_pattern.graph->Filled() == true, + Assert (sparsity_pattern.trilinos_sparsity_pattern().Filled() == true, ExcMessage("The Trilinos sparsity pattern has not been compressed")); matrix = std::auto_ptr - (new Epetra_FECrsMatrix(Copy, *sparsity_pattern.graph, false)); + (new Epetra_FECrsMatrix(Copy, sparsity_pattern.trilinos_sparsity_pattern(), + false)); compress(); } @@ -565,7 +569,8 @@ namespace TrilinosWrappers row_map = sparse_matrix.row_map; col_map = sparse_matrix.col_map; matrix = std::auto_ptr - (new Epetra_FECrsMatrix(Copy, sparse_matrix.matrix->Graph(), false)); + (new Epetra_FECrsMatrix(Copy, sparse_matrix.trilinos_sparsity_pattern(), + false)); compress(); } @@ -1135,12 +1140,13 @@ namespace TrilinosWrappers if (matrix->Filled() == false) matrix->FillComplete(col_map, row_map, true); - Assert (src.vector->Map().SameAs(matrix->DomainMap()) == true, + Assert (src.vector_partitioner().SameAs(matrix->DomainMap()) == true, ExcMessage ("Column map of matrix does not fit with vector map!")); - Assert (dst.vector->Map().SameAs(matrix->RangeMap()) == true, + Assert (dst.vector_partitioner().SameAs(matrix->RangeMap()) == true, ExcMessage ("Row map of matrix does not fit with vector map!")); - const int ierr = matrix->Multiply (false, *(src.vector), *(dst.vector)); + const int ierr = matrix->Multiply (false, src.trilinos_vector(), + dst.trilinos_vector()); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -1155,12 +1161,13 @@ namespace TrilinosWrappers if (matrix->Filled() == false) matrix->FillComplete(col_map, row_map, true); - Assert (src.vector->Map().SameAs(matrix->RangeMap()) == true, + Assert (src.vector_partitioner().SameAs(matrix->RangeMap()) == true, ExcMessage ("Column map of matrix does not fit with vector map!")); - Assert (dst.vector->Map().SameAs(matrix->DomainMap()) == true, + Assert (dst.vector_partitioner().SameAs(matrix->DomainMap()) == true, ExcMessage ("Row map of matrix does not fit with vector map!")); - const int ierr = matrix->Multiply (true, *(src.vector), *(dst.vector)); + const int ierr = matrix->Multiply (true, src.trilinos_vector(), + dst.trilinos_vector()); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc index 0bab6e5c05..ebeed0896f 100755 --- a/deal.II/lac/source/trilinos_vector.cc +++ b/deal.II/lac/source/trilinos_vector.cc @@ -235,16 +235,16 @@ namespace TrilinosWrappers Vector::import_nonlocal_data_for_fe (const TrilinosWrappers::SparseMatrix &m, const Vector &v) { - Assert (m.matrix->Filled() == true, + Assert (m.trilinos_matrix().Filled() == true, ExcMessage ("Matrix is not compressed. " "Cannot find exchange information!")); Assert (v.vector->Map().UniqueGIDs() == true, ExcMessage ("The input vector has overlapping data, " "which is not allowed.")); - if (vector->Map().SameAs(m.matrix->ColMap()) == false) + if (vector->Map().SameAs(m.col_partitioner()) == false) { - map = m.matrix->ColMap(); + map = m.col_partitioner(); vector = std::auto_ptr (new Epetra_FEVector(map)); } -- 2.39.5