From 6cb85b744a29204edfdcf272d39180c8a3417a66 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 21 Sep 2009 03:47:43 +0000 Subject: [PATCH] Make sure that we can compile lac/ with PETSc when the latter is using complex scalars. git-svn-id: https://svn.dealii.org/trunk@19482 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/lac/constraint_matrix.templates.h | 183 ++++++------- deal.II/lac/include/lac/petsc_matrix_base.h | 2 +- .../lac/petsc_parallel_sparse_matrix.h | 26 +- deal.II/lac/include/lac/petsc_solver.h | 50 ++-- deal.II/lac/include/lac/petsc_vector_base.h | 94 +++---- deal.II/lac/include/lac/vector.templates.h | 110 +++++--- deal.II/lac/source/petsc_matrix_base.cc | 110 ++++---- .../source/petsc_parallel_sparse_matrix.cc | 58 ++-- deal.II/lac/source/petsc_solver.cc | 98 +++---- deal.II/lac/source/petsc_vector_base.cc | 249 ++++++++++-------- 10 files changed, 518 insertions(+), 462 deletions(-) diff --git a/deal.II/lac/include/lac/constraint_matrix.templates.h b/deal.II/lac/include/lac/constraint_matrix.templates.h index d022682537..3fd91ec336 100644 --- a/deal.II/lac/include/lac/constraint_matrix.templates.h +++ b/deal.II/lac/include/lac/constraint_matrix.templates.h @@ -72,7 +72,7 @@ ConstraintMatrix::condense (const VectorType &uncondensed, Assert (condensed.size()+n_constraints() == uncondensed.size(), ExcDimensionMismatch(condensed.size()+n_constraints(), uncondensed.size())); - + // store for each line of the // vector its new line number after // compression. If the shift is -1, @@ -85,7 +85,7 @@ ConstraintMatrix::condense (const VectorType &uncondensed, unsigned int shift = 0; unsigned int n_rows = uncondensed.size(); - if (next_constraint == lines.end()) + if (next_constraint == lines.end()) // if no constraint is to be handled for (unsigned int row=0; row!=n_rows; ++row) new_line.push_back (row); @@ -127,7 +127,7 @@ ConstraintMatrix::condense (const VectorType &uncondensed, else // line must be distributed { - for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) + for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) condensed(new_line[next_constraint->entries[q].first]) += uncondensed(row) * next_constraint->entries[q].second; @@ -148,9 +148,10 @@ ConstraintMatrix::condense (VectorType &vec) const std::vector::const_iterator constraint_line = lines.begin(); for (; constraint_line!=lines.end(); ++constraint_line) { - for (unsigned int q=0; q!=constraint_line->entries.size(); ++q) + for (unsigned int q=0; q!=constraint_line->entries.size(); ++q) vec(constraint_line->entries[q].first) - += (vec(constraint_line->line) * + += (static_cast + (vec(constraint_line->line)) * constraint_line->entries[q].second); vec(constraint_line->line) = 0.; @@ -175,11 +176,11 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, // check whether we work on real vectors // or we just used a dummy when calling // the other function above. - const bool use_vectors = (uncondensed_vector.size() == 0 && + const bool use_vectors = (uncondensed_vector.size() == 0 && condensed_vector.size() == 0) ? false : true; const SparsityPattern &uncondensed_struct = uncondensed.get_sparsity_pattern (); - + Assert (sorted == true, ExcMatrixNotClosed()); Assert (uncondensed_struct.is_compressed() == true, ExcMatrixNotClosed()); Assert (condensed.get_sparsity_pattern().is_compressed() == true, ExcMatrixNotClosed()); @@ -210,7 +211,7 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, unsigned int shift = 0; const unsigned int n_rows = uncondensed_struct.n_rows(); - if (next_constraint == lines.end()) + if (next_constraint == lines.end()) // if no constraint is to be handled for (unsigned int row=0; row!=n_rows; ++row) new_line.push_back (row); @@ -256,7 +257,7 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1) condensed.add (new_line[row], new_line[uncondensed_struct.get_column_numbers()[j]], uncondensed.global_entry(j)); - else + else { // let c point to the // constraint of this column @@ -277,12 +278,12 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, // row of the inhomogeneous constraint in // the matrix with Gauss elimination if (use_vectors == true) - condensed_vector(new_line[row]) -= uncondensed.global_entry(j) * + condensed_vector(new_line[row]) -= uncondensed.global_entry(j) * c->inhomogeneity; } if (use_vectors == true) - condensed_vector(new_line[row]) += uncondensed_vector(row); + condensed_vector(new_line[row]) += uncondensed_vector(row); } else // line must be distributed @@ -292,12 +293,12 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, // for each column: distribute if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1) // column is not constrained - for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) + for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) condensed.add (new_line[next_constraint->entries[q].first], new_line[uncondensed_struct.get_column_numbers()[j]], uncondensed.global_entry(j) * next_constraint->entries[q].second); - + else // not only this line but // also this col is constrained @@ -307,7 +308,7 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, std::vector::const_iterator c = lines.begin(); while (c->line != uncondensed_struct.get_column_numbers()[j]) ++c; - + for (unsigned int p=0; p!=c->entries.size(); ++p) for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) condensed.add (new_line[next_constraint->entries[q].first], @@ -319,14 +320,14 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, if (use_vectors == true) for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) condensed_vector (new_line[next_constraint->entries[q].first]) - -= uncondensed.global_entry(j) * + -= uncondensed.global_entry(j) * next_constraint->entries[q].second * c->inhomogeneity; }; // condense the vector if (use_vectors == true) - for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) + for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) condensed_vector(new_line[next_constraint->entries[q].first]) += uncondensed_vector(row) * next_constraint->entries[q].second; @@ -355,7 +356,7 @@ ConstraintMatrix::condense (SparseMatrix &uncondensed, ExcNotQuadratic()); if (use_vectors == true) { - Assert (vec.size() == sparsity.n_rows(), + Assert (vec.size() == sparsity.n_rows(), ExcDimensionMismatch(vec.size(), sparsity.n_rows())); } @@ -363,7 +364,7 @@ ConstraintMatrix::condense (SparseMatrix &uncondensed, for (unsigned int i=0; i &uncondensed, // handles this index std::vector distribute (sparsity.n_rows(), numbers::invalid_unsigned_int); - + for (unsigned int c=0; c &uncondensed, entry != uncondensed.end(row); ++entry) { const unsigned int column = entry->column(); - + // end of row reached? // this should not // happen, since we only @@ -395,7 +396,7 @@ ConstraintMatrix::condense (SparseMatrix &uncondensed, // matrices! Assert (column != SparsityPattern::invalid_entry, ExcMatrixNotClosed()); - + if (distribute[column] != numbers::invalid_unsigned_int) // distribute entry at // regular row @p row @@ -417,7 +418,7 @@ ConstraintMatrix::condense (SparseMatrix &uncondensed, // row of the inhomogeneous constraint in // the matrix with Gauss elimination if (use_vectors == true) - vec(row) -= + vec(row) -= entry->value() * lines[distribute[column]].inhomogeneity; // set old value to zero @@ -451,7 +452,7 @@ ConstraintMatrix::condense (SparseMatrix &uncondensed, // old entry to zero { for (unsigned int q=0; - q!=lines[distribute[row]].entries.size(); ++q) + q!=lines[distribute[row]].entries.size(); ++q) uncondensed.add (lines[distribute[row]].entries[q].first, column, entry->value() * @@ -479,7 +480,7 @@ ConstraintMatrix::condense (SparseMatrix &uncondensed, lines[distribute[column]].entries[q].second); if (use_vectors == true) - vec(lines[distribute[row]].entries[p].first) -= + vec(lines[distribute[row]].entries[p].first) -= entry->value() * lines[distribute[row]].entries[p].second * lines[distribute[column]].inhomogeneity; } @@ -493,7 +494,7 @@ ConstraintMatrix::condense (SparseMatrix &uncondensed, // take care of vector if (use_vectors == true) { - for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q) + for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q) vec(lines[distribute[row]].entries[q].first) += (vec(row) * lines[distribute[row]].entries[q].second); @@ -516,7 +517,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, const bool use_vectors = vec.n_blocks() == 0 ? false : true; const unsigned int blocks = uncondensed.n_block_rows(); - + const BlockSparsityPattern & sparsity = uncondensed.get_sparsity_pattern (); @@ -533,7 +534,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, if (use_vectors == true) { - Assert (vec.size() == sparsity.n_rows(), + Assert (vec.size() == sparsity.n_rows(), ExcDimensionMismatch(vec.size(), sparsity.n_rows())); Assert (vec.n_blocks() == sparsity.n_block_rows(), ExcDimensionMismatch(vec.n_blocks(), sparsity.n_block_rows())); @@ -547,7 +548,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, const BlockIndices & index_mapping = sparsity.get_column_indices(); - + // store for each index whether it must be // distributed or not. If entry is // numbers::invalid_unsigned_int, @@ -557,7 +558,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, // index std::vector distribute (sparsity.n_rows(), numbers::invalid_unsigned_int); - + for (unsigned int c=0; c &uncondensed, const std::pair block_index = index_mapping.global_to_local(row); const unsigned int block_row = block_index.first; - + if (distribute[row] == numbers::invalid_unsigned_int) // regular line. loop over // all columns and see @@ -592,7 +593,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, { const unsigned int global_col = index_mapping.local_to_global(block_col,entry->column()); - + if (distribute[global_col] != numbers::invalid_unsigned_int) // distribute entry at // regular row @p row @@ -601,7 +602,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, // entry to zero { const double old_value = entry->value (); - + for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q) uncondensed.add (row, @@ -615,7 +616,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, // row of the inhomogeneous constraint in // the matrix with Gauss elimination if (use_vectors == true) - vec(row) -= entry->value() * + vec(row) -= entry->value() * lines[distribute[global_col]].inhomogeneity; entry->value() = 0.; @@ -639,7 +640,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, { const unsigned int global_col = index_mapping.local_to_global (block_col, entry->column()); - + if (distribute[global_col] == numbers::invalid_unsigned_int) // distribute @@ -653,9 +654,9 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, // zero { const double old_value = entry->value(); - + for (unsigned int q=0; - q!=lines[distribute[row]].entries.size(); ++q) + q!=lines[distribute[row]].entries.size(); ++q) uncondensed.add (lines[distribute[row]].entries[q].first, global_col, old_value * @@ -673,7 +674,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, // otherwise { const double old_value = entry->value (); - + for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p) { for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q) @@ -684,7 +685,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, lines[distribute[global_col]].entries[q].second); if (use_vectors == true) - vec(lines[distribute[row]].entries[p].first) -= + vec(lines[distribute[row]].entries[p].first) -= old_value * lines[distribute[row]].entries[p].second * lines[distribute[global_col]].inhomogeneity; } @@ -697,7 +698,7 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, // take care of vector if (use_vectors == true) { - for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q) + for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q) vec(lines[distribute[row]].entries[q].first) += (vec(row) * lines[distribute[row]].entries[q].second); @@ -749,7 +750,7 @@ distribute_local_to_global (const Vector &local_vector, // vector that tells about whether a // certain constraint exists. then, we // simply copy over the data. - const std::vector::const_iterator position = + const std::vector::const_iterator position = find_constraint(local_dof_indices[i]); if (position==lines.end()) @@ -764,7 +765,7 @@ distribute_local_to_global (const Vector &local_vector, if (val != 0) for (unsigned int j=0; j::const_iterator + const std::vector::const_iterator position_j = find_constraint(local_dof_indices[j]); if (position_j == lines.end()) @@ -821,7 +822,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, // feature further down. Vector dummy(0); distribute_local_to_global (local_matrix, dummy, local_dof_indices, - global_matrix, dummy, + global_matrix, dummy, internal::bool2type::value>()); } @@ -841,7 +842,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, // actual implementation follows further // down. distribute_local_to_global (local_matrix, local_vector, local_dof_indices, - global_matrix, global_vector, + global_matrix, global_vector, internal::bool2type::value>()); } @@ -888,7 +889,7 @@ ConstraintMatrix::distribute (const VectorType &condensed, unsigned int shift = 0; unsigned int n_rows = uncondensed.size(); - if (next_constraint == lines.end()) + if (next_constraint == lines.end()) // if no constraint is to be handled for (unsigned int row=0; row!=n_rows; ++row) old_line.push_back (row); @@ -921,7 +922,7 @@ ConstraintMatrix::distribute (const VectorType &condensed, // only evaluated so often as there are // entries in new_line[*] which tells us // which constraints exist - for (unsigned int line=0; line::const_iterator next_constraint = lines.begin(); - for (; next_constraint != lines.end(); ++next_constraint) + for (; next_constraint != lines.end(); ++next_constraint) { // fill entry in line // next_constraint.line by adding the // different contributions - double new_value = next_constraint->inhomogeneity; + typename VectorType::value_type + new_value = next_constraint->inhomogeneity; for (unsigned int i=0; ientries.size(); ++i) - new_value += (vec(next_constraint->entries[i].first) * + new_value += (static_cast + (vec(next_constraint->entries[i].first)) * next_constraint->entries[i].second); vec(next_constraint->line) = new_value; } @@ -972,7 +975,7 @@ namespace internals // which global entries (global_row) are // given rise by local entries // (local_row) or some constraints. - struct distributing + struct distributing { distributing (const unsigned int global_row = deal_II_numbers::invalid_unsigned_int, const unsigned int local_row = deal_II_numbers::invalid_unsigned_int); @@ -1070,7 +1073,7 @@ namespace internals } if (&*pos1->constraints == 0) - pos1->constraints = + pos1->constraints = new std::vector > (1,constraint); else pos1->constraints->push_back (constraint); @@ -1113,7 +1116,7 @@ namespace internals j2 = j-istep; temp = my_indices[i].global_row; templ = my_indices[i].local_row; - if (my_indices[j2].global_row > temp) + if (my_indices[j2].global_row > temp) { while ((j >= istep) && (my_indices[j2].global_row > temp)) { @@ -1154,7 +1157,7 @@ namespace internals block_starts[0] = 0; for (unsigned int i=1;i &local_matrix, // check whether we work on real vectors // or we just used a dummy when calling // the other function above. - const bool use_vectors = (local_vector.size() == 0 && + const bool use_vectors = (local_vector.size() == 0 && global_vector.size() == 0) ? false : true; Assert (local_matrix.n() == local_dof_indices.size(), @@ -1284,7 +1287,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, { have_indirect_rows = true; internals::insert_index(my_indices, position->entries[q].first, - std::make_pair + std::make_pair (local_row, position->entries[q].second)); } @@ -1313,7 +1316,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, // hanging nodes in 3d). however, in the // line below, we do actually do // something with this dof - const typename MatrixType::value_type new_diagonal + const typename MatrixType::value_type new_diagonal = (std::fabs(local_matrix(local_row,local_row)) != 0 ? std::fabs(local_matrix(local_row,local_row)) : average_diagonal); global_matrix.add(global_row, global_row, new_diagonal); @@ -1329,7 +1332,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, typedef std::vector > constraint_format; - // now do the actual job. + // now do the actual job. for (unsigned int i=0; i &local_matrix, // all on this set of dofs (saves a lot // of checks). the only check we actually // need to perform is whether the matrix - // element is zero. + // element is zero. if (have_indirect_rows == false) { Assert(loc_row < n_local_dofs, ExcInternalError()); @@ -1416,7 +1419,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, { constraint_format &constraint_j = *my_indices[j].constraints; for (unsigned int p=0; p &local_matrix, { double add_this = loc_col != deal_II_numbers::invalid_unsigned_int ? local_matrix(constraint_i[q].first, loc_col) : 0; - + if (my_indices[j].constraints != 0) { constraint_format &constraint_j = *my_indices[j].constraints; for (unsigned int p=0; p &local_matrix, } } - + // if we got some nontrivial value, // append it to the array of values. @@ -1488,7 +1491,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, if (my_indices[i].constraints != 0) { - std::vector > &constraint_i = + std::vector > &constraint_i = *my_indices[i].constraints; for (unsigned int q=0; q &local_matrix, // the other function for additional // comments. - const bool use_vectors = (local_vector.size() == 0 && + const bool use_vectors = (local_vector.size() == 0 && global_vector.size() == 0) ? false : true; Assert (local_matrix.n() == local_dof_indices.size(), @@ -1600,11 +1603,11 @@ distribute_local_to_global (const FullMatrix &local_matrix, { have_indirect_rows = true; internals::insert_index(my_indices, position->entries[q].first, - std::make_pair + std::make_pair (local_row, position->entries[q].second)); } - const typename MatrixType::value_type new_diagonal + const typename MatrixType::value_type new_diagonal = (std::fabs(local_matrix(local_row,local_row)) != 0 ? std::fabs(local_matrix(local_row,local_row)) : average_diagonal); global_matrix.add(global_row, global_row, new_diagonal); @@ -1684,12 +1687,12 @@ distribute_local_to_global (const FullMatrix &local_matrix, // constraints if (my_indices[j].constraints != 0) { - constraint_format &constraint_j = + constraint_format &constraint_j = *my_indices[j].constraints; for (unsigned int p=0; p &local_matrix, for (unsigned int q=0; q &local_matrix, if (col_val != 0) { *col_ptr++ = localized_indices[j]; - *val_ptr++ = + *val_ptr++ = static_cast(col_val); } } @@ -1762,7 +1765,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, if (my_indices[i].constraints != 0) { - std::vector > &constraint_i = + std::vector > &constraint_i = *my_indices[i].constraints; for (unsigned int q=0; q &local_matrix, } } if (val != 0) - global_vector(my_indices[i].global_row) += + global_vector(my_indices[i].global_row) += static_cast(val); } } @@ -1854,7 +1857,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, } else { - std::vector::iterator it = + std::vector::iterator it = std::lower_bound(actual_dof_indices.begin(), actual_dof_indices.end(), new_index); @@ -1943,7 +1946,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, { have_indirect_rows = true; internals::insert_index(my_indices, position->entries[q].first, - std::make_pair + std::make_pair (local_row, position->entries[q].second)); } @@ -2032,7 +2035,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, // constraints if (my_indices[j].constraints != 0) { - std::vector > &constraint_j = + std::vector > &constraint_j = *my_indices[j].constraints; for (unsigned int p=0; p &local_dof_indices, // given column. if (my_indices[i].constraints != 0) { - std::vector > &constraint_i = + std::vector > &constraint_i = *my_indices[i].constraints; for (unsigned int q=0; q &local_dof_indices, if (dof_mask[constraint_i[q].first][loc_col] == true) goto add_this_index; } - + if (my_indices[j].constraints != 0) { - std::vector > &constraint_j = + std::vector > &constraint_j = *my_indices[j].constraints; for (unsigned int p=0; p &local_dof_indices, } else { - std::vector::iterator it = + std::vector::iterator it = std::lower_bound(actual_dof_indices.begin(), actual_dof_indices.end(), new_index); @@ -2220,8 +2223,8 @@ add_entries_local_to_global (const std::vector &local_dof_indices, const unsigned int next_block_col = block_starts[block_col+1]; sparsity_pattern.block(block,block_col). add_entries(row, - index_it, - actual_dof_indices.begin() + next_block_col, + index_it, + actual_dof_indices.begin() + next_block_col, true); index_it = actual_dof_indices.begin() + next_block_col; } @@ -2288,7 +2291,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, { have_indirect_rows = true; internals::insert_index(my_indices, position->entries[q].first, - std::make_pair + std::make_pair (local_row, position->entries[q].second)); } @@ -2384,7 +2387,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, // constraints if (my_indices[j].constraints != 0) { - std::vector > + std::vector > &constraint_j = *my_indices[j].constraints; for (unsigned int p=0; p &local_dof_indices, // given column. if (my_indices[i].constraints != 0) { - std::vector > + std::vector > &constraint_i = *my_indices[i].constraints; for (unsigned int q=0; q &local_dof_indices, if (dof_mask[constraint_i[q].first][loc_col] == true) goto add_this_index; } - + if (my_indices[j].constraints != 0) { - std::vector > + std::vector > &constraint_j = *my_indices[j].constraints; for (unsigned int p=0; p &local_dof_indices, // that accumulated under the given // process into the global matrix row and // into the vector - sparsity_pattern.block(block, block_col).add_entries(row, - cols.begin(), + sparsity_pattern.block(block, block_col).add_entries(row, + cols.begin(), col_ptr, true); } diff --git a/deal.II/lac/include/lac/petsc_matrix_base.h b/deal.II/lac/include/lac/petsc_matrix_base.h index 592b81777c..17121547a0 100644 --- a/deal.II/lac/include/lac/petsc_matrix_base.h +++ b/deal.II/lac/include/lac/petsc_matrix_base.h @@ -323,7 +323,7 @@ namespace PETScWrappers * previously used. */ MatrixBase & - operator = (const double d); + operator = (const value_type d); /** * Release all memory and return * to a state just like after diff --git a/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h b/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h index 38128e8440..c58591aea4 100644 --- a/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h +++ b/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h @@ -34,7 +34,7 @@ namespace PETScWrappers { namespace MPI { - + /** * Implementation of a parallel sparse matrix class based on PETSC, with rows * of the matrix distributed across an MPI network. All the functionality is @@ -104,7 +104,7 @@ namespace PETScWrappers * which. Thus, it is not sufficient to just count the number of degrees of * freedom that belong to a particular process, but you have to have the * numbers for all processes available at all processes. - * + * * @ingroup PETScWrappers * @ingroup Matrix1 * @author Wolfgang Bangerth, 2004 @@ -147,7 +147,7 @@ namespace PETScWrappers * empty matrix. */ SparseMatrix (); - + /** * Create a sparse matrix of * dimensions @p m times @p n, with @@ -164,7 +164,7 @@ namespace PETScWrappers * local_row and @p local_columns * parameters, see the class * documentation. - * + * * The @p is_symmetric flag determines * whether we should tell PETSc that * the matrix is going to be symmetric @@ -197,7 +197,7 @@ namespace PETScWrappers * local_row and @p local_columns * parameters, see the class * documentation. - * + * * Just as for the other * constructors: PETSc is able to * cope with the situation that more @@ -226,7 +226,7 @@ namespace PETScWrappers const unsigned int local_columns, const std::vector &row_lengths, const bool is_symmetric = false); - + /** * Initialize using the given * sparsity pattern with @@ -238,7 +238,7 @@ namespace PETScWrappers * local_columns_per_process * parameters, see the class * documentation. - * + * * Note that PETSc can be very slow * if you do not provide it with a * good estimate of the lengths of @@ -291,7 +291,7 @@ namespace PETScWrappers * but keep the sparsity pattern * previously used. */ - SparseMatrix & operator = (const double d); + SparseMatrix & operator = (const value_type d); /** * Throw away the present matrix and @@ -366,14 +366,14 @@ namespace PETScWrappers const std::vector &local_columns_per_process, const unsigned int this_process, const bool preset_nonzero_locations = true); - + /** * Return a reference to the MPI * communicator object in use with * this matrix. */ virtual const MPI_Comm & get_mpi_communicator () const; - + /** @addtogroup Exceptions * @{ */ /** @@ -383,7 +383,7 @@ namespace PETScWrappers int, int, << "The number of local rows " << arg1 << " must be larger than the total number of rows " << arg2); - //@} + //@} private: /** @@ -431,7 +431,7 @@ namespace PETScWrappers * To allow calling protected * prepare_add() and * prepare_set(). - */ + */ friend class BlockMatrixBase; }; @@ -440,7 +440,7 @@ namespace PETScWrappers // -------- template and inline functions ---------- inline - const MPI_Comm & + const MPI_Comm & SparseMatrix::get_mpi_communicator () const { return communicator; diff --git a/deal.II/lac/include/lac/petsc_solver.h b/deal.II/lac/include/lac/petsc_solver.h index 380652cfe4..45deb8ea3b 100644 --- a/deal.II/lac/include/lac/petsc_solver.h +++ b/deal.II/lac/include/lac/petsc_solver.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 2004, 2005, 2006, 2007, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -32,8 +32,8 @@ namespace PETScWrappers class MatrixBase; class VectorBase; class PreconditionerBase; - - + + /** * Base class for solver classes using the PETSc solvers. Since solvers in * PETSc are selected based on flags passed to a generic solver object, @@ -62,7 +62,7 @@ namespace PETScWrappers * fact that the default argument to the solver classes, @p PETSC_COMM_SELF, * is the appropriate argument for the sequential case (which is why it is the * default argument), so this error only shows up in parallel mode. - * + * * @ingroup PETScWrappers * @author Wolfgang Bangerth, 2004 */ @@ -87,7 +87,7 @@ namespace PETScWrappers */ SolverBase (SolverControl &cn, const MPI_Comm &mpi_communicator); - + /** * Destructor. */ @@ -165,7 +165,7 @@ namespace PETScWrappers int convergence_test (KSP ksp, const int iteration, - const PetscScalar residual_norm, + const PetscReal residual_norm, KSPConvergedReason *reason, void *solver_control); @@ -204,7 +204,7 @@ namespace PETScWrappers * Destructor */ ~SolverData (); - + #if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2) /** * A PETSc solver object. @@ -245,7 +245,7 @@ namespace PETScWrappers * pipe additional data to the * solver. */ - struct AdditionalData + struct AdditionalData { /** * Constructor. By default, @@ -253,13 +253,13 @@ namespace PETScWrappers * to one. */ AdditionalData (const double omega = 1); - + /** * Relaxation parameter. */ double omega; }; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -327,7 +327,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -395,7 +395,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -463,7 +463,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -538,7 +538,7 @@ namespace PETScWrappers */ AdditionalData (const unsigned int restart_parameter = 30, const bool right_preconditioning = false); - + /** * Maximum number of * tmp vectors. @@ -551,7 +551,7 @@ namespace PETScWrappers */ bool right_preconditioning; }; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -619,7 +619,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -685,7 +685,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -753,7 +753,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -826,7 +826,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -894,7 +894,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -962,7 +962,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -1010,14 +1010,14 @@ namespace PETScWrappers */ virtual void set_solver_type (KSP &ksp) const; }; - + /** * An implementation of the solver interface using the PETSc PREONLY * solver. Actually this is NOT a real solution algorithm. Its only * purpose is to provide a solver object, when the preconditioner * should be used as real solver. It is very useful in conjunction with - * the complete LU decomposition preconditioner PreconditionLU , + * the complete LU decomposition preconditioner PreconditionLU , * which in conjunction with this solver class becomes a direct solver. * * @ingroup PETScWrappers @@ -1033,7 +1033,7 @@ namespace PETScWrappers */ struct AdditionalData {}; - + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -1065,7 +1065,7 @@ namespace PETScWrappers SolverPreOnly (SolverControl &cn, const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); - + protected: /** * Store a copy of the flags for this diff --git a/deal.II/lac/include/lac/petsc_vector_base.h b/deal.II/lac/include/lac/petsc_vector_base.h index 5cf061e7cd..b2f920801f 100644 --- a/deal.II/lac/include/lac/petsc_vector_base.h +++ b/deal.II/lac/include/lac/petsc_vector_base.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 2004, 2005, 2006, 2007, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -47,7 +47,7 @@ namespace PETScWrappers /** * @cond internal */ - + /** * A namespace for internal implementation details of the PETScWrapper * members. @@ -84,7 +84,7 @@ namespace PETScWrappers */ VectorReference (const VectorBase &vector, const unsigned int index); - + public: /** * This looks like a copy operator, @@ -107,7 +107,7 @@ namespace PETScWrappers * the assignment operator const. */ const VectorReference & operator = (const VectorReference &r) const; - + /** * Set the referenced element of the * vector to s. @@ -162,7 +162,7 @@ namespace PETScWrappers << " of a distributed vector, but only elements " << arg2 << " through " << arg3 << " are stored locally and can be accessed."); - + private: /** * Point to the vector we are @@ -187,7 +187,7 @@ namespace PETScWrappers /** * @endcond */ - + /** * Base class for all vector classes that are implemented on top of the PETSc @@ -224,7 +224,7 @@ namespace PETScWrappers * class. */ typedef PetscScalar value_type; - typedef PetscScalar real_type; + typedef PetscReal real_type; typedef size_t size_type; typedef internal::VectorReference reference; typedef const internal::VectorReference const_reference; @@ -235,7 +235,7 @@ namespace PETScWrappers * to initialize the data. */ VectorBase (); - + /** * Copy constructor. Sets the dimension * to that of the given vector, and @@ -279,8 +279,8 @@ namespace PETScWrappers * other values is deprecated and may * be disallowed in the future. */ - VectorBase & operator = (const PetscScalar s); - + VectorBase & operator = (const PetscScalar s); + /** * Test for equality. This function * assumes that the present vector and @@ -290,7 +290,7 @@ namespace PETScWrappers * sense anyway. */ bool operator == (const VectorBase &v) const; - + /** * Test for inequality. This function * assumes that the present vector and @@ -383,7 +383,7 @@ namespace PETScWrappers */ void add (const std::vector &indices, const std::vector &values); - + /** * Return the scalar product of two * vectors. The vectors must have the @@ -421,7 +421,7 @@ namespace PETScWrappers * powers of the absolute values * of the elements. */ - real_type lp_norm (const PetscScalar p) const; + real_type lp_norm (const real_type p) const; /** * Maximum absolute value of the @@ -461,7 +461,7 @@ namespace PETScWrappers * The model for this function * should be similer to add (). */ - VectorBase & mult (); + VectorBase & mult (); /** * Return whether the vector contains @@ -482,13 +482,13 @@ namespace PETScWrappers * really all positive (or zero). */ bool is_non_negative () const; - + /** * Multiply the entire vector by a * fixed factor. */ VectorBase & operator *= (const PetscScalar factor); - + /** * Divide the entire vector by a * fixed factor. @@ -513,26 +513,26 @@ namespace PETScWrappers * scalar and not a vector. */ void add (const PetscScalar s); - + /** * Simple vector addition, equal to the * operator +=. */ void add (const VectorBase &V); - + /** * Simple addition of a multiple of a * vector, i.e. *this += a*V. */ void add (const PetscScalar a, const VectorBase &V); - + /** * Multiple addition of scaled vectors, * i.e. *this += a*V+b*W. */ void add (const PetscScalar a, const VectorBase &V, const PetscScalar b, const VectorBase &W); - + /** * Scaling and simple vector addition, * i.e. @@ -540,7 +540,7 @@ namespace PETScWrappers */ void sadd (const PetscScalar s, const VectorBase &V); - + /** * Scaling and simple addition, i.e. * *this = s*(*this)+a*V. @@ -548,7 +548,7 @@ namespace PETScWrappers void sadd (const PetscScalar s, const PetscScalar a, const VectorBase &V); - + /** * Scaling and multiple addition. */ @@ -557,7 +557,7 @@ namespace PETScWrappers const VectorBase &V, const PetscScalar b, const VectorBase &W); - + /** * Scaling and multiple addition. * *this = s*(*this)+a*V + b*W + c*X. @@ -566,10 +566,10 @@ namespace PETScWrappers const PetscScalar a, const VectorBase &V, const PetscScalar b, - const VectorBase &W, + const VectorBase &W, const PetscScalar c, const VectorBase &X); - + /** * Scale each element of this * vector by the corresponding @@ -580,12 +580,12 @@ namespace PETScWrappers * diagonal scaling matrix. */ void scale (const VectorBase &scaling_factors); - + /** * Assignment *this = a*V. */ void equ (const PetscScalar a, const VectorBase &V); - + /** * Assignment *this = a*V + b*W. */ @@ -654,7 +654,7 @@ namespace PETScWrappers * to standard functions. */ void swap (VectorBase &v); - + /** * Conversion operator to gain access * to the underlying PETSc type. If you @@ -674,7 +674,7 @@ namespace PETScWrappers * for this class). */ unsigned int memory_consumption () const; - + protected: /** * A generic vector object in @@ -741,7 +741,7 @@ namespace PETScWrappers -// ------------------- inline and template functions -------------- +// ------------------- inline and template functions -------------- /** * Global function @p swap which overloads the default implementation @@ -784,7 +784,7 @@ namespace PETScWrappers } - + inline const VectorReference & VectorReference::operator = (const PetscScalar &value) const @@ -798,15 +798,15 @@ namespace PETScWrappers ierr = VecAssemblyEnd (vector); AssertThrow (ierr == 0, ExcPETScError(ierr)); } - + const signed int petsc_i = index; - + const int ierr = VecSetValues (vector, 1, &petsc_i, &value, INSERT_VALUES); AssertThrow (ierr == 0, ExcPETScError(ierr)); vector.last_action = VectorBase::LastAction::insert; - + return *this; } @@ -837,14 +837,14 @@ namespace PETScWrappers // zero if (value == PetscScalar()) return *this; - + // use the PETSc function to add something const signed int petsc_i = index; const int ierr = VecSetValues (vector, 1, &petsc_i, &value, ADD_VALUES); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return *this; } @@ -865,7 +865,7 @@ namespace PETScWrappers vector.last_action = VectorBase::LastAction::add; } - + // we have to do above actions in any // case to be consistent with the MPI // communication model (see the @@ -882,7 +882,7 @@ namespace PETScWrappers const int ierr = VecSetValues (vector, 1, &petsc_i, &subtractand, ADD_VALUES); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return *this; } @@ -903,7 +903,7 @@ namespace PETScWrappers vector.last_action = VectorBase::LastAction::insert; } - + // we have to do above actions in any // case to be consistent with the MPI // communication model (see the @@ -918,11 +918,11 @@ namespace PETScWrappers const PetscScalar new_value = static_cast(*this) * value; - + const int ierr = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return *this; } @@ -943,7 +943,7 @@ namespace PETScWrappers vector.last_action = VectorBase::LastAction::insert; } - + // we have to do above actions in any // case to be consistent with the MPI // communication model (see the @@ -958,15 +958,15 @@ namespace PETScWrappers const PetscScalar new_value = static_cast(*this) / value; - + const int ierr = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return *this; } } - + inline @@ -977,7 +977,7 @@ namespace PETScWrappers const int ierr = VecGetOwnershipRange (static_cast(vector), &begin, &end); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return ((index >= static_cast(begin)) && (index < static_cast(end))); } @@ -999,7 +999,7 @@ namespace PETScWrappers { return static_cast(internal::VectorReference (*this, index)); } - + #endif // DOXYGEN } diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index c1d9f089ee..9e817471ee 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -49,7 +49,7 @@ namespace internal { Assert (false, ExcMessage ("Complex numbers do not have an ordering.")); - + return false; } @@ -63,7 +63,7 @@ namespace internal else std::printf (" %5.2f", double(t)); } - + template @@ -75,10 +75,41 @@ namespace internal else std::printf (" %5.2f+%5.2fi", double(t.real()), double(t.imag())); + } + + // call std::copy, except for in + // the case where we want to copy + // from std::complex to a + // non-complex type + template + void copy (const T *begin, + const T *end, + U *dest) + { + std::copy (begin, end, dest); + } + + template + void copy (const std::complex *begin, + const std::complex *end, + std::complex *dest) + { + std::copy (begin, end, dest); + } + + template + void copy (const std::complex *, + const std::complex *, + U *) + { + Assert (false, ExcMessage ("Can't convert a vector of complex numbers " + "into a vector of reals/doubles")); } } + + template Vector::Vector (const Vector& v) : @@ -119,6 +150,7 @@ Vector::Vector (const Vector& v) #ifdef DEAL_II_USE_PETSC + template Vector::Vector (const PETScWrappers::Vector &v) : @@ -137,8 +169,8 @@ Vector::Vector (const PETScWrappers::Vector &v) PetscScalar *start_ptr; int ierr = VecGetArray (static_cast(v), &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - - std::copy (start_ptr, start_ptr+vec_size, begin()); + + internal::copy (start_ptr, start_ptr+vec_size, begin()); // restore the representation of the // vector @@ -198,7 +230,7 @@ Vector::Vector (const TrilinosWrappers::MPI::Vector &v) 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()); } } @@ -224,7 +256,7 @@ Vector::Vector (const TrilinosWrappers::Vector &v) int ierr = v.trilinos_vector().ExtractView (&start_ptr); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - + std::copy (start_ptr[0], start_ptr[0]+vec_size, begin()); } } @@ -272,7 +304,7 @@ bool Vector::is_non_negative () const { Assert (vec_size!=0, ExcEmptyObject()); - + for (unsigned int i=0; i Number Vector::operator * (const Vector& v) const { Assert (vec_size!=0, ExcEmptyObject()); - + if (PointerComparison::equal (this, &v)) return norm_sqr(); - + Assert (vec_size == v.size(), ExcDimensionMismatch(vec_size, v.size())); - + Number sum = 0; // multiply the two vectors. we have to @@ -304,7 +336,7 @@ Number Vector::operator * (const Vector& v) const // is not defined by default for (unsigned int i=0; i::conjugate(v.val[i])); - + return sum; } @@ -319,7 +351,7 @@ Vector::norm_sqr () const for (unsigned int i=0; i::abs_square(val[i]); - + return sum; } @@ -333,7 +365,7 @@ Number Vector::mean_value () const for (unsigned int i=0; i::lp_norm (const real_type p) const for (unsigned int i=0; i::abs(val[i]), p); - + return std::pow(sum, static_cast(1./p)); } @@ -454,9 +486,9 @@ template void Vector::add (const Number a, const Vector& v, const Number b, const Vector& w) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(b), + Assert (numbers::is_finite(b), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); Assert (vec_size!=0, ExcEmptyObject()); @@ -477,7 +509,7 @@ template void Vector::sadd (const Number x, const Vector& v) { - Assert (numbers::is_finite(x), + Assert (numbers::is_finite(x), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); Assert (vec_size!=0, ExcEmptyObject()); @@ -498,11 +530,11 @@ void Vector::sadd (const Number x, const Number a, const Vector& v, const Number b, const Vector& w) { - Assert (numbers::is_finite(x), + Assert (numbers::is_finite(x), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(b), + Assert (numbers::is_finite(b), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); Assert (vec_size!=0, ExcEmptyObject()); @@ -564,7 +596,7 @@ template void Vector::equ (const Number a, const Vector& u) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); Assert (vec_size!=0, ExcEmptyObject()); @@ -584,7 +616,7 @@ template void Vector::equ (const Number a, const Vector& u) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); Assert (vec_size!=0, ExcEmptyObject()); @@ -606,9 +638,9 @@ template void Vector::equ (const Number a, const Vector& u, const Number b, const Vector& v) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(b), + Assert (numbers::is_finite(b), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); Assert (vec_size!=0, ExcEmptyObject()); @@ -677,7 +709,7 @@ Vector::operator = (const BlockVector& v) for (unsigned int b=0; b::operator = (const PETScWrappers::Vector &v) PetscScalar *start_ptr; int ierr = VecGetArray (static_cast(v), &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - - std::copy (start_ptr, start_ptr+vec_size, begin()); + + internal::copy (start_ptr, start_ptr+vec_size, begin()); // restore the representation of the // vector @@ -758,7 +790,7 @@ Vector::operator = (const TrilinosWrappers::Vector &v) TrilinosScalar **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()); } @@ -813,7 +845,7 @@ void Vector::print (std::ostream &out, std::ios::fmtflags old_flags = out.flags(); unsigned int old_precision = out.precision (precision); - + out.precision (precision); if (scientific) out.setf (std::ios::scientific, std::ios::floatfield); @@ -827,7 +859,7 @@ void Vector::print (std::ostream &out, for (unsigned int i=0; i::block_write (std::ostream &out) const // environment const unsigned int sz = size(); char buf[16]; - + std::sprintf(buf, "%d", sz); std::strcat(buf, "\n["); - + out.write(buf, std::strlen(buf)); out.write (reinterpret_cast(begin()), reinterpret_cast(end()) - reinterpret_cast(begin())); - + // out << ']'; const char outro = ']'; out.write (&outro, 1); - + AssertThrow (out, ExcIO()); } @@ -875,24 +907,24 @@ void Vector::block_read (std::istream &in) unsigned int sz; char buf[16]; - + in.getline(buf,16,'\n'); sz=std::atoi(buf); - + // fast initialization, since the // data elements are overwritten anyway - reinit (sz, true); + reinit (sz, true); char c; // in >> c; in.read (&c, 1); AssertThrow (c=='[', ExcIO()); - + in.read (reinterpret_cast(begin()), reinterpret_cast(end()) - reinterpret_cast(begin())); - + // in >> c; in.read (&c, 1); AssertThrow (c==']', ExcIO()); diff --git a/deal.II/lac/source/petsc_matrix_base.cc b/deal.II/lac/source/petsc_matrix_base.cc index b6d27f8648..2508eaf432 100644 --- a/deal.II/lac/source/petsc_matrix_base.cc +++ b/deal.II/lac/source/petsc_matrix_base.cc @@ -41,7 +41,7 @@ namespace PETScWrappers return; } - + // otherwise first flush PETSc caches matrix->compress (); @@ -58,7 +58,7 @@ namespace PETScWrappers const int *colnums; const PetscScalar *values; #endif - + int ierr; ierr = MatGetRow(*matrix, this->a_row, &ncols, &colnums, &values); AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr)); @@ -73,7 +73,7 @@ namespace PETScWrappers colnum_cache.reset (new std::vector (colnums, colnums+ncols)); value_cache.reset (new std::vector (values, values+ncols)); - + // and finally restore the matrix ierr = MatRestoreRow(*matrix, this->a_row, &ncols, &colnums, &values); AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr)); @@ -81,18 +81,18 @@ namespace PETScWrappers } - + MatrixBase::MatrixBase () : last_action (LastAction::none) {} - + MatrixBase::~MatrixBase () { const int ierr = MatDestroy (matrix); - AssertThrow (ierr == 0, ExcPETScError(ierr)); + AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -102,7 +102,7 @@ namespace PETScWrappers { // destroy the matrix... int ierr = MatDestroy (matrix); - AssertThrow (ierr == 0, ExcPETScError(ierr)); + AssertThrow (ierr == 0, ExcPETScError(ierr)); // ...and replace it by an empty // sequential matrix const int m=0, n=0, n_nonzero_per_row=0; @@ -110,20 +110,20 @@ namespace PETScWrappers 0, &matrix); AssertThrow (ierr == 0, ExcPETScError(ierr)); } - + MatrixBase & - MatrixBase::operator = (const double d) + MatrixBase::operator = (const value_type d) { - Assert (d==0, ExcScalarAssignmentOnlyForZeroValue()); + Assert (d==value_type(), ExcScalarAssignmentOnlyForZeroValue()); // flush previously cached elements. this // seems to be necessary since petsc // 2.2.1, at least for parallel vectors // (see test bits/petsc_64) compress (); - + const int ierr = MatZeroEntries (matrix); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -147,10 +147,10 @@ namespace PETScWrappers #else const PetscInt petsc_row = row; #endif - + IS index_set; ISCreateGeneral (get_mpi_communicator(), 1, &petsc_row, &index_set); - + #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR <= 2) const int ierr = MatZeroRows(matrix, index_set, &new_diag_value); @@ -158,11 +158,11 @@ namespace PETScWrappers const int ierr = MatZeroRowsIS(matrix, index_set, new_diag_value); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); ISDestroy (index_set); - + compress (); } @@ -190,7 +190,7 @@ namespace PETScWrappers IS index_set; ISCreateGeneral (get_mpi_communicator(), rows.size(), &petsc_rows[0], &index_set); - + #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR <= 2) const int ierr = MatZeroRows(matrix, index_set, &new_diag_value); @@ -198,15 +198,15 @@ namespace PETScWrappers const int ierr = MatZeroRowsIS(matrix, index_set, new_diag_value); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); ISDestroy (index_set); - + compress (); } - - + + PetscScalar MatrixBase::el (const unsigned int i, @@ -215,7 +215,7 @@ namespace PETScWrappers const signed int petsc_i = i; const signed int petsc_j = j; PetscScalar value; - + const int ierr = MatGetValues (matrix, 1, &petsc_i, 1, &petsc_j, &value); @@ -230,14 +230,14 @@ namespace PETScWrappers MatrixBase::diag_element (const unsigned int i) const { Assert (m() == n(), ExcNotQuadratic()); - + // this doesn't seem to work any // different than any other element return el(i,i); } - - + + void MatrixBase::compress () { @@ -265,7 +265,7 @@ namespace PETScWrappers return n_rows; } - + unsigned int @@ -291,7 +291,7 @@ namespace PETScWrappers } - + std::pair MatrixBase::local_range () const { @@ -317,7 +317,7 @@ namespace PETScWrappers } - + unsigned int MatrixBase:: row_length (const unsigned int row) const @@ -364,21 +364,21 @@ namespace PETScWrappers MatrixBase::l1_norm () const { PetscReal result; - + const int ierr = MatNorm (matrix, NORM_1, &result); AssertThrow (ierr == 0, ExcPETScError(ierr)); return result; } - - + + PetscReal MatrixBase::linfty_norm () const { PetscReal result; - + const int ierr = MatNorm (matrix, NORM_INFINITY, &result); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -392,7 +392,7 @@ namespace PETScWrappers MatrixBase::frobenius_norm () const { PetscReal result; - + const int ierr = MatNorm (matrix, NORM_FROBENIUS, &result); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -421,7 +421,7 @@ namespace PETScWrappers MatrixBase::operator /= (const PetscScalar a) { const PetscScalar factor = 1./a; - + #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) const int ierr = MatScale (&factor, matrix); #else @@ -439,13 +439,13 @@ namespace PETScWrappers const VectorBase &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); - + const int ierr = MatMult (matrix, src, dst); AssertThrow (ierr == 0, ExcPETScError(ierr)); } - - + + void MatrixBase::Tvmult (VectorBase &dst, const VectorBase &src) const @@ -455,15 +455,15 @@ namespace PETScWrappers const int ierr = MatMultTranspose (matrix, src, dst); AssertThrow (ierr == 0, ExcPETScError(ierr)); } - - + + void MatrixBase::vmult_add (VectorBase &dst, const VectorBase &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); - + const int ierr = MatMultAdd (matrix, src, dst, dst); AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -475,7 +475,7 @@ namespace PETScWrappers const VectorBase &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); - + const int ierr = MatMultTransposeAdd (matrix, src, dst, dst); AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -489,9 +489,9 @@ namespace PETScWrappers vmult (tmp, v); return tmp*v; } - - + + PetscScalar MatrixBase::matrix_scalar_product (const VectorBase &u, const VectorBase &v) const @@ -500,9 +500,9 @@ namespace PETScWrappers vmult (tmp, v); return u*tmp; } - - + + PetscScalar MatrixBase::residual (VectorBase &dst, const VectorBase &x, @@ -514,23 +514,23 @@ namespace PETScWrappers vmult (dst, x); dst -= b; dst *= -1; - + return dst.l2_norm(); } - - + + MatrixBase::operator const Mat () const { return matrix; - } + } void MatrixBase::transpose () { int ierr; -#if (PETSC_VERSION_MAJOR <= 2) +#if (PETSC_VERSION_MAJOR <= 2) ierr = MatTranspose(matrix, PETSC_NULL); #else ierr = MatTranspose(matrix, MAT_REUSE_MATRIX, &matrix); @@ -540,17 +540,17 @@ namespace PETScWrappers } PetscTruth - MatrixBase::is_symmetric (const double tolerance) + MatrixBase::is_symmetric (const double tolerance) { PetscTruth truth; // First flush PETSc caches compress (); MatIsSymmetric (matrix, tolerance, &truth); return truth; - } + } PetscTruth - MatrixBase::is_hermitian (const double tolerance) + MatrixBase::is_hermitian (const double tolerance) { PetscTruth truth; // First flush PETSc caches @@ -559,13 +559,13 @@ namespace PETScWrappers #if (PETSC_VERSION_MAJOR <= 2) // avoid warning about unused variables (void) tolerance; - + MatIsHermitian (matrix, &truth); #else MatIsHermitian (matrix, tolerance, &truth); #endif return truth; - } + } void MatrixBase::write_ascii () @@ -574,7 +574,7 @@ namespace PETScWrappers compress (); // Set options -#if (PETSC_VERSION_MAJOR <= 2) +#if (PETSC_VERSION_MAJOR <= 2) PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DEFAULT); #else diff --git a/deal.II/lac/source/petsc_parallel_sparse_matrix.cc b/deal.II/lac/source/petsc_parallel_sparse_matrix.cc index 0b9c1c2a76..0d804eb214 100644 --- a/deal.II/lac/source/petsc_parallel_sparse_matrix.cc +++ b/deal.II/lac/source/petsc_parallel_sparse_matrix.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2004, 2005, 2006, 2008 by the deal.II authors +// Copyright (C) 2004, 2005, 2006, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -27,7 +27,7 @@ namespace PETScWrappers { namespace MPI { - + SparseMatrix::SparseMatrix () { // just like for vectors: since we @@ -89,17 +89,17 @@ namespace PETScWrappers local_columns_per_process, this_process, preset_nonzero_locations); } - + SparseMatrix & - SparseMatrix::operator = (const double d) + SparseMatrix::operator = (const value_type d) { MatrixBase::operator = (d); return *this; } - + void SparseMatrix::reinit (const MPI_Comm &communicator, @@ -115,8 +115,8 @@ namespace PETScWrappers // get rid of old matrix and generate a // new one const int ierr = MatDestroy (matrix); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - + AssertThrow (ierr == 0, ExcPETScError(ierr)); + do_reinit (m, n, local_rows, local_columns, n_nonzero_per_row, is_symmetric); } @@ -133,14 +133,14 @@ namespace PETScWrappers const bool is_symmetric) { this->communicator = communicator; - + // get rid of old matrix and generate a // new one const int ierr = MatDestroy (matrix); - AssertThrow (ierr == 0, ExcPETScError(ierr)); + AssertThrow (ierr == 0, ExcPETScError(ierr)); do_reinit (m, n, local_rows, local_columns, row_lengths, is_symmetric); - } + } @@ -155,18 +155,18 @@ namespace PETScWrappers const bool preset_nonzero_locations) { this->communicator = communicator; - + // get rid of old matrix and generate a // new one const int ierr = MatDestroy (matrix); - AssertThrow (ierr == 0, ExcPETScError(ierr)); + AssertThrow (ierr == 0, ExcPETScError(ierr)); do_reinit (sparsity_pattern, local_rows_per_process, local_columns_per_process, this_process, preset_nonzero_locations); } - + void SparseMatrix::do_reinit (const unsigned int m, @@ -177,7 +177,7 @@ namespace PETScWrappers const bool is_symmetric) { Assert (local_rows <= m, ExcLocalRowsTooLarge (local_rows, m)); - + // use the call sequence indicating only // a maximal number of elements per row // for all rows globally @@ -192,7 +192,7 @@ namespace PETScWrappers // set symmetric flag, if so requested if (is_symmetric == true) { -#if (PETSC_VERSION_MAJOR <= 2) +#if (PETSC_VERSION_MAJOR <= 2) const int ierr = MatSetOption (matrix, MAT_SYMMETRIC); #else @@ -228,7 +228,7 @@ namespace PETScWrappers for (unsigned int i=0; i= local_col_start) && (column < local_col_end)) ++row_lengths_in_window[row-local_row_start]; @@ -359,7 +359,7 @@ namespace PETScWrappers // the whole thing: #if (PETSC_VERSION_MAJOR <= 2) && \ ((PETSC_VERSION_MINOR < 2) || \ - ((PETSC_VERSION_MINOR == 2) && (PETSC_VERSION_SUBMINOR == 0))) + ((PETSC_VERSION_MINOR == 2) && (PETSC_VERSION_SUBMINOR == 0))) std::vector row_entries; std::vector row_values; @@ -369,7 +369,7 @@ namespace PETScWrappers row_values.resize (sparsity_pattern.row_length(i), 0.0); for (unsigned int j=0; j values (sparsity_pattern.max_entries_per_row(), 1.); - + for (unsigned int i=local_row_start; i(&solver_control)); - + // then do the real work: set up solver // internal data and solve the // system. this could be joined in one @@ -120,7 +120,7 @@ namespace PETScWrappers int iterations = 0; ierr = SLESSetUp (sles, b, x); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + ierr = SLESSolve (sles, b, x, &iterations); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -135,7 +135,7 @@ namespace PETScWrappers if (solver_data == 0) { solver_data.reset (new SolverData()); - + ierr = KSPCreate (mpi_communicator, &solver_data->ksp); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -152,7 +152,7 @@ namespace PETScWrappers // object set the type of // preconditioner set_solver_type (solver_data->ksp); - + ierr = KSPGetPC (solver_data->ksp, &solver_data->pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -163,7 +163,7 @@ namespace PETScWrappers // checks with the solver_control // object we have in this object for // convergence -#if (PETSC_VERSION_MAJOR <= 2) +#if (PETSC_VERSION_MAJOR <= 2) KSPSetConvergenceTest (solver_data->ksp, &convergence_test, reinterpret_cast(&solver_control)); #else @@ -173,7 +173,7 @@ namespace PETScWrappers #endif } - + // then do the real work: set up solver // internal data and solve the // system. unfortunately, the call @@ -186,7 +186,7 @@ namespace PETScWrappers AssertThrow (ierr == 0, ExcPETScError(ierr)); ierr = KSPSetUp (solver_data->ksp); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + ierr = KSPSolve (solver_data->ksp); AssertThrow (ierr == 0, ExcPETScError(ierr)); # else @@ -195,7 +195,7 @@ namespace PETScWrappers ierr = KSPSolve (solver_data->ksp, b, x); AssertThrow (ierr == 0, ExcPETScError(ierr)); -# endif +# endif #endif // destroy solver object @@ -216,38 +216,38 @@ namespace PETScWrappers { return solver_control; } - + int SolverBase::convergence_test (KSP /*ksp*/, const int iteration, - const PetscScalar residual_norm, + const PetscReal residual_norm, KSPConvergedReason *reason, void *solver_control_x) { SolverControl &solver_control = *reinterpret_cast(solver_control_x); - + const SolverControl::State state = solver_control.check (iteration, residual_norm); - + switch (state) { case ::dealii::SolverControl::iterate: *reason = KSP_CONVERGED_ITERATING; break; - + case ::dealii::SolverControl::success: *reason = static_cast(1); break; - + case ::dealii::SolverControl::failure: if (solver_control.last_step() > solver_control.max_steps()) *reason = KSP_DIVERGED_ITS; else *reason = KSP_DIVERGED_DTOL; break; - + default: Assert (false, ExcNotImplemented()); } @@ -255,8 +255,8 @@ namespace PETScWrappers // return without failure return 0; } - - + + /* ---------------------- SolverRichardson ------------------------ */ @@ -266,8 +266,8 @@ namespace PETScWrappers omega (omega) {} - - + + SolverRichardson::SolverRichardson (SolverControl &cn, const MPI_Comm &mpi_communicator, const AdditionalData &data) @@ -313,7 +313,7 @@ namespace PETScWrappers KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(), PETSC_DEFAULT, this->solver_control.max_steps()+1); } - + /* ---------------------- SolverChebychev ------------------------ */ @@ -343,7 +343,7 @@ namespace PETScWrappers KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverCG ------------------------ */ SolverCG::SolverCG (SolverControl &cn, @@ -371,7 +371,7 @@ namespace PETScWrappers // solution vector. do so here as well: KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverBiCG ------------------------ */ @@ -401,7 +401,7 @@ namespace PETScWrappers KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverGMRES ------------------------ */ SolverGMRES::AdditionalData:: @@ -412,8 +412,8 @@ namespace PETScWrappers right_preconditioning (right_preconditioning) {} - - + + SolverGMRES::SolverGMRES (SolverControl &cn, const MPI_Comm &mpi_communicator, const AdditionalData &data) @@ -454,7 +454,7 @@ namespace PETScWrappers // so rather expand their macros by hand, // and do some equally nasty stuff that at // least doesn't yield warnings... - int (*fun_ptr)(KSP,int); + int (*fun_ptr)(KSP,int); ierr = PetscObjectQueryFunction((PetscObject)(ksp), "KSPGMRESSetRestart_C", (void (**)())&fun_ptr); @@ -476,7 +476,7 @@ namespace PETScWrappers // solution vector. do so here as well: KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverBicgstab ------------------------ */ @@ -506,7 +506,7 @@ namespace PETScWrappers KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverCGS ------------------------ */ SolverCGS::SolverCGS (SolverControl &cn, @@ -534,7 +534,7 @@ namespace PETScWrappers // solution vector. do so here as well: KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverTFQMR ------------------------ */ @@ -564,7 +564,7 @@ namespace PETScWrappers KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverTCQMR ------------------------ */ SolverTCQMR::SolverTCQMR (SolverControl &cn, @@ -593,7 +593,7 @@ namespace PETScWrappers KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverCR ------------------------ */ SolverCR::SolverCR (SolverControl &cn, @@ -622,7 +622,7 @@ namespace PETScWrappers KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverLSQR ------------------------ */ SolverLSQR::SolverLSQR (SolverControl &cn, @@ -650,7 +650,7 @@ namespace PETScWrappers // solution vector. do so here as well: KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); } - + /* ---------------------- SolverPreOnly ------------------------ */ @@ -676,11 +676,11 @@ namespace PETScWrappers // The KSPPREONLY solver of // PETSc never calls the convergence - // monitor, which leads to failure + // monitor, which leads to failure // even when everything was ok. - // Therefore the SolverControl status + // Therefore the SolverControl status // is set to some nice values, which - // guarantee a nice result at the end + // guarantee a nice result at the end // of the solution process. solver_control.check (1, 0.0); diff --git a/deal.II/lac/source/petsc_vector_base.cc b/deal.II/lac/source/petsc_vector_base.cc index 95436b78e8..a359b72f87 100644 --- a/deal.II/lac/source/petsc_vector_base.cc +++ b/deal.II/lac/source/petsc_vector_base.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -30,7 +30,7 @@ namespace PETScWrappers { Assert (index < vector.size(), ExcIndexRange (index, 0, vector.size())); - + // this is clumsy: there is no simple // way in PETSc to read an element from // a vector, i.e. there is no function @@ -49,12 +49,12 @@ namespace PETScWrappers int ierr = VecGetArray (static_cast(vector), &ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + const PetscScalar value = *(ptr+index); - + ierr = VecRestoreArray (static_cast(vector), &ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return value; } else if (dynamic_cast(&vector) != 0) @@ -76,12 +76,12 @@ namespace PETScWrappers PetscScalar *ptr; ierr = VecGetArray (static_cast(vector), &ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + const PetscScalar value = *(ptr+index-begin); - + ierr = VecRestoreArray (static_cast(vector), &ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return value; } else @@ -89,16 +89,16 @@ namespace PETScWrappers // exists there? Assert (false, ExcInternalError()); return -1e20; - } + } } - + VectorBase::VectorBase () : last_action (LastAction::none) {} - - + + VectorBase::VectorBase (const VectorBase &v) : last_action (LastAction::none) @@ -110,7 +110,7 @@ namespace PETScWrappers AssertThrow (ierr == 0, ExcPETScError(ierr)); } - + VectorBase::~VectorBase () { const int ierr = VecDestroy (vector); @@ -123,7 +123,7 @@ namespace PETScWrappers VectorBase::operator = (const PetscScalar s) { - Assert (numbers::is_finite(s), + Assert (numbers::is_finite(s), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); // flush previously cached elements. this @@ -131,13 +131,13 @@ namespace PETScWrappers // 2.2.1, at least for parallel vectors // (see test bits/petsc_65) compress (); - + #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) const int ierr = VecSet (&s, vector); #else const int ierr = VecSet (vector, s); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); return *this; @@ -150,9 +150,9 @@ namespace PETScWrappers { Assert (size() == v.size(), ExcDimensionMismatch(size(), v.size())); - + PetscTruth flag; - + const int ierr = VecEqual (vector, v.vector, &flag); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -160,7 +160,7 @@ namespace PETScWrappers } - + bool VectorBase::operator != (const VectorBase &v) const { @@ -168,7 +168,7 @@ namespace PETScWrappers ExcDimensionMismatch(size(), v.size())); PetscTruth flag; - + const int ierr = VecEqual (vector, v.vector, &flag); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -176,7 +176,7 @@ namespace PETScWrappers } - + unsigned int VectorBase::size () const { @@ -229,7 +229,7 @@ namespace PETScWrappers { do_set_add_operation(indices, values, true); } - + PetscScalar @@ -245,7 +245,7 @@ namespace PETScWrappers return result; } - + void @@ -263,7 +263,7 @@ namespace PETScWrappers VectorBase::real_type VectorBase::norm_sqr () const { - const PetscScalar d = l2_norm(); + const real_type d = l2_norm(); return d*d; } @@ -300,48 +300,48 @@ namespace PETScWrappers // add up remaining elements while (ptr != start_ptr+size()) sum0 += *ptr++; - - mean = (sum0+sum1+sum2+sum3)/size(); + + mean = (sum0+sum1+sum2+sum3)/static_cast(size()); } - + // restore the representation of the // vector ierr = VecRestoreArray (vector, &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return mean; } - + VectorBase::real_type VectorBase::l1_norm () const { - PetscScalar d; + real_type d; const int ierr = VecNorm (vector, NORM_1, &d); AssertThrow (ierr == 0, ExcPETScError(ierr)); return d; } - + VectorBase::real_type VectorBase::l2_norm () const { - PetscScalar d; + real_type d; const int ierr = VecNorm (vector, NORM_2, &d); AssertThrow (ierr == 0, ExcPETScError(ierr)); return d; } - + VectorBase::real_type - VectorBase::lp_norm (const PetscScalar p) const + VectorBase::lp_norm (const real_type p) const { // get a representation of the vector and // loop over all the elements @@ -349,12 +349,12 @@ namespace PETScWrappers int ierr = VecGetArray (vector, &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - PetscScalar norm = 0; + real_type norm = 0; { - PetscScalar sum0 = 0, - sum1 = 0, - sum2 = 0, - sum3 = 0; + real_type sum0 = 0, + sum1 = 0, + sum2 = 0, + sum3 = 0; // use modern processors better by // allowing pipelined commands to be @@ -363,24 +363,24 @@ namespace PETScWrappers const PetscScalar * eptr = ptr + (size()/4)*4; while (ptr!=eptr) { - sum0 += std::pow(std::fabs(*ptr++), p); - sum1 += std::pow(std::fabs(*ptr++), p); - sum2 += std::pow(std::fabs(*ptr++), p); - sum3 += std::pow(std::fabs(*ptr++), p); - }; + sum0 += std::pow(numbers::NumberTraits::abs(*ptr++), p); + sum1 += std::pow(numbers::NumberTraits::abs(*ptr++), p); + sum2 += std::pow(numbers::NumberTraits::abs(*ptr++), p); + sum3 += std::pow(numbers::NumberTraits::abs(*ptr++), p); + } // add up remaining elements while (ptr != start_ptr+size()) - sum0 += std::pow(std::fabs(*ptr++), p); - + sum0 += std::pow(numbers::NumberTraits::abs(*ptr++), p); + norm = std::pow(sum0+sum1+sum2+sum3, - static_cast(1./p)); + 1./p); } - + // restore the representation of the // vector ierr = VecRestoreArray (vector, &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return norm; } @@ -389,7 +389,7 @@ namespace PETScWrappers VectorBase::real_type VectorBase::linfty_norm () const { - PetscScalar d; + real_type d; const int ierr = VecNorm (vector, NORM_INFINITY, &d); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -402,7 +402,7 @@ namespace PETScWrappers VectorBase::real_type VectorBase::normalize () const { - PetscScalar d; + real_type d; const int ierr = VecNormalize (vector, &d); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -415,8 +415,8 @@ namespace PETScWrappers VectorBase::real_type VectorBase::max () const { - PetscInt p; - PetscScalar d; + PetscInt p; + real_type d; const int ierr = VecMax (vector, &p, &d); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -427,7 +427,7 @@ namespace PETScWrappers VectorBase & - VectorBase::abs () + VectorBase::abs () { const int ierr = VecAbs (vector); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -438,7 +438,7 @@ namespace PETScWrappers VectorBase & - VectorBase::conjugate () + VectorBase::conjugate () { const int ierr = VecConjugate (vector); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -449,7 +449,7 @@ namespace PETScWrappers VectorBase & - VectorBase::mult () + VectorBase::mult () { const int ierr = VecPointwiseMult (vector,vector,vector); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -473,23 +473,44 @@ namespace PETScWrappers bool flag = true; while (ptr != eptr) { - if (*ptr != 0) + if (*ptr != value_type()) { flag = false; break; } ++ptr; } - + // restore the representation of the // vector ierr = VecRestoreArray (vector, &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return flag; } + namespace internal + { + template + bool is_non_negative (const T &t) + { + return t >= 0; + } + + + + template + bool is_non_negative (const std::complex &) + { + Assert (false, + ExcMessage ("You can't ask a complex value " + "whether it is non-negative.")) + return true; + } + } + + bool VectorBase::is_non_negative () const @@ -505,29 +526,29 @@ namespace PETScWrappers bool flag = true; while (ptr != eptr) { - if (*ptr < 0.0) + if (! internal::is_non_negative(*ptr)) { flag = false; break; } ++ptr; } - + // restore the representation of the // vector ierr = VecRestoreArray (vector, &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + return flag; } - + VectorBase & VectorBase::operator *= (const PetscScalar a) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) @@ -535,7 +556,7 @@ namespace PETScWrappers #else const int ierr = VecScale (vector, a); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); return *this; @@ -547,12 +568,12 @@ namespace PETScWrappers VectorBase::operator /= (const PetscScalar a) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); const PetscScalar factor = 1./a; - - Assert (numbers::is_finite(factor), + + Assert (numbers::is_finite(factor), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) @@ -572,12 +593,12 @@ namespace PETScWrappers VectorBase::operator += (const VectorBase &v) { #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) - const PetscScalar one = 1.0; + const PetscScalar one = 1.0; const int ierr = VecAXPY (&one, v, vector); #else const int ierr = VecAXPY (vector, 1, v); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); return *this; @@ -589,12 +610,12 @@ namespace PETScWrappers VectorBase::operator -= (const VectorBase &v) { #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) - const PetscScalar minus_one = -1.0; + const PetscScalar minus_one = -1.0; const int ierr = VecAXPY (&minus_one, v, vector); #else const int ierr = VecAXPY (vector, -1, v); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); return *this; @@ -606,7 +627,7 @@ namespace PETScWrappers VectorBase::add (const PetscScalar s) { - Assert (numbers::is_finite(s), + Assert (numbers::is_finite(s), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) @@ -614,10 +635,10 @@ namespace PETScWrappers #else const int ierr = VecShift (vector, s); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); } - + void @@ -625,7 +646,7 @@ namespace PETScWrappers { *this += v; } - + void @@ -633,7 +654,7 @@ namespace PETScWrappers const VectorBase &v) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) @@ -641,10 +662,10 @@ namespace PETScWrappers #else const int ierr = VecAXPY (vector, a, v); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); } - + void @@ -654,14 +675,14 @@ namespace PETScWrappers const VectorBase &w) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(b), + Assert (numbers::is_finite(b), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); const PetscScalar weights[2] = {a,b}; Vec addends[2] = {v.vector, w.vector}; - + #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) const int ierr = VecMAXPY (2, weights, vector, addends); #else @@ -670,7 +691,7 @@ namespace PETScWrappers AssertThrow (ierr == 0, ExcPETScError(ierr)); } - + void @@ -678,7 +699,7 @@ namespace PETScWrappers const VectorBase &v) { - Assert (numbers::is_finite(s), + Assert (numbers::is_finite(s), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) @@ -686,10 +707,10 @@ namespace PETScWrappers #else const int ierr = VecAYPX (vector, s, v); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); } - + void @@ -698,9 +719,9 @@ namespace PETScWrappers const VectorBase &v) { - Assert (numbers::is_finite(s), + Assert (numbers::is_finite(s), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); // there is nothing like a AXPAY @@ -709,7 +730,7 @@ namespace PETScWrappers *this *= s; add (a,v); } - + void @@ -719,27 +740,27 @@ namespace PETScWrappers const PetscScalar b, const VectorBase &w) { - - Assert (numbers::is_finite(s), + + Assert (numbers::is_finite(s), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(b), + Assert (numbers::is_finite(b), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); // there is no operation like MAXPAY, so // do it in two steps *this *= s; - + const PetscScalar weights[2] = {a,b}; Vec addends[2] = {v.vector,w.vector}; - + #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) const int ierr = VecMAXPY (2, weights, vector, addends); #else const int ierr = VecMAXPY (vector, 2, weights, addends); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -754,14 +775,14 @@ namespace PETScWrappers const PetscScalar c, const VectorBase &x) { - - Assert (numbers::is_finite(s), + + Assert (numbers::is_finite(s), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(b), + Assert (numbers::is_finite(b), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(c), + Assert (numbers::is_finite(c), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); // there is no operation like MAXPAY, so @@ -770,13 +791,13 @@ namespace PETScWrappers const PetscScalar weights[3] = {a,b,c}; Vec addends[3] = {v.vector, w.vector, x.vector}; - + #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3) const int ierr = VecMAXPY (3, weights, vector, addends); #else const int ierr = VecMAXPY (vector, 3, weights, addends); #endif - + AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -797,7 +818,7 @@ namespace PETScWrappers const VectorBase &v) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); Assert (size() == v.size(), @@ -811,7 +832,7 @@ namespace PETScWrappers *this *= a; } - + void @@ -821,9 +842,9 @@ namespace PETScWrappers const VectorBase &w) { - Assert (numbers::is_finite(a), + Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); - Assert (numbers::is_finite(b), + Assert (numbers::is_finite(b), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); Assert (size() == v.size(), @@ -877,29 +898,29 @@ namespace PETScWrappers if (across) for (unsigned int i=0; i(val[i]) << ' '; + out << val[i] << ' '; else for (unsigned int i=0; i(val[i]) << std::endl; + out << val[i] << std::endl; out << std::endl; // restore the representation of the // vector ierr = VecRestoreArray (vector, &val); AssertThrow (ierr == 0, ExcPETScError(ierr)); - + AssertThrow (out, ExcIO()); } - + void VectorBase::swap (VectorBase &v) { const int ierr = VecSwap (vector, v.vector); AssertThrow (ierr == 0, ExcPETScError(ierr)); } - + VectorBase::operator const Vec & () const @@ -915,7 +936,7 @@ namespace PETScWrappers AssertThrow(false, ExcNotImplemented() ); return 0; } - + void @@ -966,7 +987,7 @@ namespace PETScWrappers last_action = LastAction::insert; } - + } DEAL_II_NAMESPACE_CLOSE -- 2.39.5