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,
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);
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;
std::vector<ConstraintLine>::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<typename VectorType::value_type>
+ (vec(constraint_line->line)) *
constraint_line->entries[q].second);
vec(constraint_line->line) = 0.;
// 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());
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);
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
// 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
// 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
std::vector<ConstraintLine>::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],
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;
ExcNotQuadratic());
if (use_vectors == true)
{
- Assert (vec.size() == sparsity.n_rows(),
+ Assert (vec.size() == sparsity.n_rows(),
ExcDimensionMismatch(vec.size(), sparsity.n_rows()));
}
for (unsigned int i=0; i<uncondensed.m(); ++i)
average_diagonal += std::fabs (uncondensed.diag_element(i));
average_diagonal /= uncondensed.m();
-
+
// store for each index whether it must be
// distributed or not. If entry is
// invalid_unsigned_int, no distribution is
// handles this index
std::vector<unsigned int> distribute (sparsity.n_rows(),
numbers::invalid_unsigned_int);
-
+
for (unsigned int c=0; c<lines.size(); ++c)
distribute[lines[c].line] = c;
entry != uncondensed.end(row); ++entry)
{
const unsigned int column = entry->column();
-
+
// end of row reached?
// this should not
// happen, since we only
// matrices!
Assert (column != SparsityPattern::invalid_entry,
ExcMatrixNotClosed());
-
+
if (distribute[column] != numbers::invalid_unsigned_int)
// distribute entry at
// regular row @p row
// 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
// 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() *
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;
}
// 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);
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 ();
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()));
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,
// index
std::vector<unsigned int> distribute (sparsity.n_rows(),
numbers::invalid_unsigned_int);
-
+
for (unsigned int c=0; c<lines.size(); ++c)
distribute[lines[c].line] = c;
const std::pair<unsigned int,unsigned int>
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
{
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
// 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,
// 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.;
{
const unsigned int global_col
= index_mapping.local_to_global (block_col, entry->column());
-
+
if (distribute[global_col] ==
numbers::invalid_unsigned_int)
// distribute
// 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 *
// 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)
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;
}
// 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);
// vector that tells about whether a
// certain constraint exists. then, we
// simply copy over the data.
- const std::vector<ConstraintLine>::const_iterator position =
+ const std::vector<ConstraintLine>::const_iterator position =
find_constraint(local_dof_indices[i]);
if (position==lines.end())
if (val != 0)
for (unsigned int j=0; j<n_local_dofs; ++j)
{
- const std::vector<ConstraintLine>::const_iterator
+ const std::vector<ConstraintLine>::const_iterator
position_j = find_constraint(local_dof_indices[j]);
if (position_j == lines.end())
// feature further down.
Vector<double> dummy(0);
distribute_local_to_global (local_matrix, dummy, local_dof_indices,
- global_matrix, dummy,
+ global_matrix, dummy,
internal::bool2type<IsBlockMatrix<MatrixType>::value>());
}
// 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<IsBlockMatrix<MatrixType>::value>());
}
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);
// only evaluated so often as there are
// entries in new_line[*] which tells us
// which constraints exist
- for (unsigned int line=0; line<uncondensed.size(); ++line)
+ for (unsigned int line=0; line<uncondensed.size(); ++line)
if (old_line[line] != -1)
// line was not condensed away
uncondensed(line) = condensed(old_line[line]);
Assert (sorted == true, ExcMatrixNotClosed());
std::vector<ConstraintLine>::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; i<next_constraint->entries.size(); ++i)
- new_value += (vec(next_constraint->entries[i].first) *
+ new_value += (static_cast<typename VectorType::value_type>
+ (vec(next_constraint->entries[i].first)) *
next_constraint->entries[i].second);
vec(next_constraint->line) = new_value;
}
// 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);
}
if (&*pos1->constraints == 0)
- pos1->constraints =
+ pos1->constraints =
new std::vector<std::pair<unsigned int,double> > (1,constraint);
else
pos1->constraints->push_back (constraint);
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))
{
block_starts[0] = 0;
for (unsigned int i=1;i<num_blocks;++i)
{
- row_iterator first_block =
+ row_iterator first_block =
std::lower_bound (col_indices,
row_indices.end(),
block_object.get_row_indices().block_start(i));
// 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(),
{
have_indirect_rows = true;
internals::insert_index(my_indices, position->entries[q].first,
- std::make_pair<unsigned int,double>
+ std::make_pair<unsigned int,double>
(local_row, position->entries[q].second));
}
// 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);
typedef std::vector<std::pair<unsigned int,double> > constraint_format;
- // now do the actual job.
+ // now do the actual job.
for (unsigned int i=0; i<n_actual_dofs; ++i)
{
const unsigned int row = my_indices[i].global_row;
// 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());
{
constraint_format &constraint_j = *my_indices[j].constraints;
for (unsigned int p=0; p<constraint_j.size(); ++p)
- col_val += matrix_ptr[constraint_j[p].first]
+ col_val += matrix_ptr[constraint_j[p].first]
*
constraint_j[p].second;
}
{
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<constraint_j.size(); ++p)
add_this += local_matrix(constraint_i[q].first,
- constraint_j[p].first)
+ constraint_j[p].first)
*
constraint_j[p].second;
}
}
}
-
+
// if we got some nontrivial value,
// append it to the array of values.
if (my_indices[i].constraints != 0)
{
- std::vector<std::pair<unsigned int,double> > &constraint_i =
+ std::vector<std::pair<unsigned int,double> > &constraint_i =
*my_indices[i].constraints;
for (unsigned int q=0; q<constraint_i.size(); ++q)
// 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(),
{
have_indirect_rows = true;
internals::insert_index(my_indices, position->entries[q].first,
- std::make_pair<unsigned int,double>
+ std::make_pair<unsigned int,double>
(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);
// 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<constraint_j.size(); ++p)
col_val += local_matrix(loc_row,
- constraint_j[p].first)
+ constraint_j[p].first)
*
constraint_j[p].second;
}
for (unsigned int q=0; q<constraint_i.size(); ++q)
{
- double add_this =
+ 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 =
+ constraint_format &constraint_j =
*my_indices[j].constraints;
for (unsigned int p=0; p<constraint_j.size(); ++p)
add_this += local_matrix(constraint_i[q].first,
- constraint_j[p].first)
+ constraint_j[p].first)
*
constraint_j[p].second;
}
if (col_val != 0)
{
*col_ptr++ = localized_indices[j];
- *val_ptr++ =
+ *val_ptr++ =
static_cast<typename MatrixType::value_type>(col_val);
}
}
if (my_indices[i].constraints != 0)
{
- std::vector<std::pair<unsigned int,double> > &constraint_i =
+ std::vector<std::pair<unsigned int,double> > &constraint_i =
*my_indices[i].constraints;
for (unsigned int q=0; q<constraint_i.size(); ++q)
}
}
if (val != 0)
- global_vector(my_indices[i].global_row) +=
+ global_vector(my_indices[i].global_row) +=
static_cast<typename VectorType::value_type>(val);
}
}
}
else
{
- std::vector<unsigned int>::iterator it =
+ std::vector<unsigned int>::iterator it =
std::lower_bound(actual_dof_indices.begin(),
actual_dof_indices.end(),
new_index);
{
have_indirect_rows = true;
internals::insert_index(my_indices, position->entries[q].first,
- std::make_pair<unsigned int,double>
+ std::make_pair<unsigned int,double>
(local_row, position->entries[q].second));
}
// constraints
if (my_indices[j].constraints != 0)
{
- std::vector<std::pair<unsigned int,double> > &constraint_j =
+ std::vector<std::pair<unsigned int,double> > &constraint_j =
*my_indices[j].constraints;
for (unsigned int p=0; p<constraint_j.size(); ++p)
// given column.
if (my_indices[i].constraints != 0)
{
- std::vector<std::pair<unsigned int,double> > &constraint_i =
+ std::vector<std::pair<unsigned int,double> > &constraint_i =
*my_indices[i].constraints;
for (unsigned int q=0; q<constraint_i.size(); ++q)
{
if (dof_mask[constraint_i[q].first][loc_col] == true)
goto add_this_index;
}
-
+
if (my_indices[j].constraints != 0)
{
- std::vector<std::pair<unsigned int,double> > &constraint_j =
+ std::vector<std::pair<unsigned int,double> > &constraint_j =
*my_indices[j].constraints;
for (unsigned int p=0; p<constraint_j.size(); ++p)
}
else
{
- std::vector<unsigned int>::iterator it =
+ std::vector<unsigned int>::iterator it =
std::lower_bound(actual_dof_indices.begin(),
actual_dof_indices.end(),
new_index);
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;
}
{
have_indirect_rows = true;
internals::insert_index(my_indices, position->entries[q].first,
- std::make_pair<unsigned int,double>
+ std::make_pair<unsigned int,double>
(local_row, position->entries[q].second));
}
// constraints
if (my_indices[j].constraints != 0)
{
- std::vector<std::pair<unsigned int,double> >
+ std::vector<std::pair<unsigned int,double> >
&constraint_j = *my_indices[j].constraints;
for (unsigned int p=0; p<constraint_j.size(); ++p)
// given column.
if (my_indices[i].constraints != 0)
{
- std::vector<std::pair<unsigned int,double> >
+ std::vector<std::pair<unsigned int,double> >
&constraint_i = *my_indices[i].constraints;
for (unsigned int q=0; q<constraint_i.size(); ++q)
{
if (dof_mask[constraint_i[q].first][loc_col] == true)
goto add_this_index;
}
-
+
if (my_indices[j].constraints != 0)
{
- std::vector<std::pair<unsigned int,double> >
+ std::vector<std::pair<unsigned int,double> >
&constraint_j = *my_indices[j].constraints;
for (unsigned int p=0; p<constraint_j.size(); ++p)
// 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);
}
* previously used.
*/
MatrixBase &
- operator = (const double d);
+ operator = (const value_type d);
/**
* Release all memory and return
* to a state just like after
{
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
* 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
* empty matrix.
*/
SparseMatrix ();
-
+
/**
* Create a sparse matrix of
* dimensions @p m times @p n, with
* 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
* 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
const unsigned int local_columns,
const std::vector<unsigned int> &row_lengths,
const bool is_symmetric = false);
-
+
/**
* Initialize using the given
* sparsity pattern with
* 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
* but keep the sparsity pattern
* previously used.
*/
- SparseMatrix & operator = (const double d);
+ SparseMatrix & operator = (const value_type d);
/**
* Throw away the present matrix and
const std::vector<unsigned int> &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
* @{ */
/**
int, int,
<< "The number of local rows " << arg1
<< " must be larger than the total number of rows " << arg2);
- //@}
+ //@}
private:
/**
* To allow calling protected
* prepare_add() and
* prepare_set().
- */
+ */
friend class BlockMatrixBase<SparseMatrix>;
};
// -------- template and inline functions ----------
inline
- const MPI_Comm &
+ const MPI_Comm &
SparseMatrix::get_mpi_communicator () const
{
return communicator;
// $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
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,
* 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
*/
*/
SolverBase (SolverControl &cn,
const MPI_Comm &mpi_communicator);
-
+
/**
* Destructor.
*/
int
convergence_test (KSP ksp,
const int iteration,
- const PetscScalar residual_norm,
+ const PetscReal residual_norm,
KSPConvergedReason *reason,
void *solver_control);
* Destructor
*/
~SolverData ();
-
+
#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2)
/**
* A PETSc solver object.
* pipe additional data to the
* solver.
*/
- struct AdditionalData
+ struct AdditionalData
{
/**
* Constructor. By default,
* to one.
*/
AdditionalData (const double omega = 1);
-
+
/**
* Relaxation parameter.
*/
double omega;
};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
AdditionalData (const unsigned int restart_parameter = 30,
const bool right_preconditioning = false);
-
+
/**
* Maximum number of
* tmp vectors.
*/
bool right_preconditioning;
};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
*/
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 <tt> PreconditionLU </tt>,
+ * the complete LU decomposition preconditioner <tt> PreconditionLU </tt>,
* which in conjunction with this solver class becomes a direct solver.
*
* @ingroup PETScWrappers
*/
struct AdditionalData
{};
-
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
SolverPreOnly (SolverControl &cn,
const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
const AdditionalData &data = AdditionalData());
-
+
protected:
/**
* Store a copy of the flags for this
// $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
/**
* @cond internal
*/
-
+
/**
* A namespace for internal implementation details of the PETScWrapper
* members.
*/
VectorReference (const VectorBase &vector,
const unsigned int index);
-
+
public:
/**
* This looks like a copy operator,
* the assignment operator const.
*/
const VectorReference & operator = (const VectorReference &r) const;
-
+
/**
* Set the referenced element of the
* vector to <tt>s</tt>.
<< " of a distributed vector, but only elements "
<< arg2 << " through " << arg3
<< " are stored locally and can be accessed.");
-
+
private:
/**
* Point to the vector we are
/**
* @endcond
*/
-
+
/**
* Base class for all vector classes that are implemented on top of the PETSc
* 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;
* to initialize the data.
*/
VectorBase ();
-
+
/**
* Copy constructor. Sets the dimension
* to that of the given vector, and
* 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
* sense anyway.
*/
bool operator == (const VectorBase &v) const;
-
+
/**
* Test for inequality. This function
* assumes that the present vector and
*/
void add (const std::vector<unsigned int> &indices,
const std::vector<PetscScalar> &values);
-
+
/**
* Return the scalar product of two
* vectors. The vectors must have the
* 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
* The model for this function
* should be similer to add ().
*/
- VectorBase & mult ();
+ VectorBase & mult ();
/**
* Return whether the vector contains
* 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.
* scalar and not a vector.
*/
void add (const PetscScalar s);
-
+
/**
* Simple vector addition, equal to the
* <tt>operator +=</tt>.
*/
void add (const VectorBase &V);
-
+
/**
* Simple addition of a multiple of a
* vector, i.e. <tt>*this += a*V</tt>.
*/
void add (const PetscScalar a, const VectorBase &V);
-
+
/**
* Multiple addition of scaled vectors,
* i.e. <tt>*this += a*V+b*W</tt>.
*/
void add (const PetscScalar a, const VectorBase &V,
const PetscScalar b, const VectorBase &W);
-
+
/**
* Scaling and simple vector addition,
* i.e.
*/
void sadd (const PetscScalar s,
const VectorBase &V);
-
+
/**
* Scaling and simple addition, i.e.
* <tt>*this = s*(*this)+a*V</tt>.
void sadd (const PetscScalar s,
const PetscScalar a,
const VectorBase &V);
-
+
/**
* Scaling and multiple addition.
*/
const VectorBase &V,
const PetscScalar b,
const VectorBase &W);
-
+
/**
* Scaling and multiple addition.
* <tt>*this = s*(*this)+a*V + b*W + c*X</tt>.
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
* diagonal scaling matrix.
*/
void scale (const VectorBase &scaling_factors);
-
+
/**
* Assignment <tt>*this = a*V</tt>.
*/
void equ (const PetscScalar a, const VectorBase &V);
-
+
/**
* Assignment <tt>*this = a*V + b*W</tt>.
*/
* to standard functions.
*/
void swap (VectorBase &v);
-
+
/**
* Conversion operator to gain access
* to the underlying PETSc type. If you
* for this class).
*/
unsigned int memory_consumption () const;
-
+
protected:
/**
* A generic vector object in
-// ------------------- inline and template functions --------------
+// ------------------- inline and template functions --------------
/**
* Global function @p swap which overloads the default implementation
}
-
+
inline
const VectorReference &
VectorReference::operator = (const PetscScalar &value) const
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;
}
// 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;
}
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
const int ierr
= VecSetValues (vector, 1, &petsc_i, &subtractand, ADD_VALUES);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
return *this;
}
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
const PetscScalar new_value
= static_cast<PetscScalar>(*this) * value;
-
+
const int ierr
= VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
return *this;
}
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
const PetscScalar new_value
= static_cast<PetscScalar>(*this) / value;
-
+
const int ierr
= VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
return *this;
}
}
-
+
inline
const int ierr = VecGetOwnershipRange (static_cast<const Vec &>(vector),
&begin, &end);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
return ((index >= static_cast<unsigned int>(begin)) &&
(index < static_cast<unsigned int>(end)));
}
{
return static_cast<PetscScalar>(internal::VectorReference (*this, index));
}
-
+
#endif // DOXYGEN
}
{
Assert (false,
ExcMessage ("Complex numbers do not have an ordering."));
-
+
return false;
}
else
std::printf (" %5.2f", double(t));
}
-
+
template <typename T>
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 <typename T, typename U>
+ void copy (const T *begin,
+ const T *end,
+ U *dest)
+ {
+ std::copy (begin, end, dest);
+ }
+
+ template <typename T, typename U>
+ void copy (const std::complex<T> *begin,
+ const std::complex<T> *end,
+ std::complex<U> *dest)
+ {
+ std::copy (begin, end, dest);
+ }
+
+ template <typename T, typename U>
+ void copy (const std::complex<T> *,
+ const std::complex<T> *,
+ U *)
+ {
+ Assert (false, ExcMessage ("Can't convert a vector of complex numbers "
+ "into a vector of reals/doubles"));
}
}
+
+
template <typename Number>
Vector<Number>::Vector (const Vector<Number>& v)
:
#ifdef DEAL_II_USE_PETSC
+
template <typename Number>
Vector<Number>::Vector (const PETScWrappers::Vector &v)
:
PetscScalar *start_ptr;
int ierr = VecGetArray (static_cast<const Vec&>(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
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());
}
}
int ierr = v.trilinos_vector().ExtractView (&start_ptr);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
+
std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());
}
}
Vector<Number>::is_non_negative () const
{
Assert (vec_size!=0, ExcEmptyObject());
-
+
for (unsigned int i=0; i<vec_size; ++i)
if ( ! internal::is_non_negative (val[i]))
return false;
Number Vector<Number>::operator * (const Vector<Number2>& 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
// is not defined by default
for (unsigned int i=0; i<vec_size; ++i)
sum += val[i] * Number(numbers::NumberTraits<Number2>::conjugate(v.val[i]));
-
+
return sum;
}
for (unsigned int i=0; i<vec_size; ++i)
sum += numbers::NumberTraits<Number>::abs_square(val[i]);
-
+
return sum;
}
for (unsigned int i=0; i<vec_size; ++i)
sum += val[i];
-
+
return sum / real_type(size());
}
for (unsigned int i=0; i<vec_size; ++i)
sum += std::pow (numbers::NumberTraits<Number>::abs(val[i]), p);
-
+
return std::pow(sum, static_cast<real_type>(1./p));
}
void Vector<Number>::add (const Number a, const Vector<Number>& v,
const Number b, const Vector<Number>& 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());
void Vector<Number>::sadd (const Number x,
const Vector<Number>& 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());
const Vector<Number>& v, const Number b,
const Vector<Number>& 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());
void Vector<Number>::equ (const Number a,
const Vector<Number>& 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());
void Vector<Number>::equ (const Number a,
const Vector<Number2>& 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());
void Vector<Number>::equ (const Number a, const Vector<Number>& u,
const Number b, const Vector<Number>& 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());
for (unsigned int b=0; b<v.n_blocks(); ++b)
for (unsigned int i=0; i<v.block(b).size(); ++i, ++this_index)
val[this_index] = v.block(b)(i);
-
+
return *this;
}
PetscScalar *start_ptr;
int ierr = VecGetArray (static_cast<const Vec&>(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
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());
}
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);
for (unsigned int i=0; i<size(); ++i)
out << val[i] << std::endl;
out << std::endl;
-
+
AssertThrow (out, ExcIO());
// reset output format
out.flags (old_flags);
// 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<const char*>(begin()),
reinterpret_cast<const char*>(end())
- reinterpret_cast<const char*>(begin()));
-
+
// out << ']';
const char outro = ']';
out.write (&outro, 1);
-
+
AssertThrow (out, ExcIO());
}
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<char*>(begin()),
reinterpret_cast<const char*>(end())
- reinterpret_cast<const char*>(begin()));
-
+
// in >> c;
in.read (&c, 1);
AssertThrow (c==']', ExcIO());
return;
}
-
+
// otherwise first flush PETSc caches
matrix->compress ();
const int *colnums;
const PetscScalar *values;
#endif
-
+
int ierr;
ierr = MatGetRow(*matrix, this->a_row, &ncols, &colnums, &values);
AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr));
colnum_cache.reset (new std::vector<unsigned int> (colnums,
colnums+ncols));
value_cache.reset (new std::vector<PetscScalar> (values, values+ncols));
-
+
// and finally restore the matrix
ierr = MatRestoreRow(*matrix, this->a_row, &ncols, &colnums, &values);
AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr));
}
-
+
MatrixBase::MatrixBase ()
:
last_action (LastAction::none)
{}
-
+
MatrixBase::~MatrixBase ()
{
const int ierr = MatDestroy (matrix);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
}
{
// 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;
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));
#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);
const int ierr
= MatZeroRowsIS(matrix, index_set, new_diag_value);
#endif
-
+
AssertThrow (ierr == 0, ExcPETScError(ierr));
ISDestroy (index_set);
-
+
compress ();
}
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);
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,
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);
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 ()
{
return n_rows;
}
-
+
unsigned int
}
-
+
std::pair<unsigned int, unsigned int>
MatrixBase::local_range () const
{
}
-
+
unsigned int
MatrixBase::
row_length (const unsigned int row) const
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));
MatrixBase::frobenius_norm () const
{
PetscReal result;
-
+
const int ierr
= MatNorm (matrix, NORM_FROBENIUS, &result);
AssertThrow (ierr == 0, ExcPETScError(ierr));
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
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
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));
}
const VectorBase &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
-
+
const int ierr = MatMultTransposeAdd (matrix, src, dst, dst);
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
vmult (tmp, v);
return tmp*v;
}
-
-
+
+
PetscScalar
MatrixBase::matrix_scalar_product (const VectorBase &u,
const VectorBase &v) const
vmult (tmp, v);
return u*tmp;
}
-
-
+
+
PetscScalar
MatrixBase::residual (VectorBase &dst,
const VectorBase &x,
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);
}
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
#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 ()
compress ();
// Set options
-#if (PETSC_VERSION_MAJOR <= 2)
+#if (PETSC_VERSION_MAJOR <= 2)
PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_ASCII_DEFAULT);
#else
// $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
{
namespace MPI
{
-
+
SparseMatrix::SparseMatrix ()
{
// just like for vectors: since we
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,
// 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);
}
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);
- }
+ }
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,
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
// 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
for (unsigned int i=0; i<row_lengths.size(); ++i)
Assert(row_lengths[i]<=local_columns,
ExcIndexRange(row_lengths[i], 1, local_columns+1));
-
+
// use the call sequence indicating a
// maximal number of elements for each
// row individually. annoyingly, we
// 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
#endif
AssertThrow (ierr == 0, ExcPETScError(ierr));
- }
+ }
}
for (unsigned int c=0; c<sparsity_pattern.row_length(row); ++c)
{
const unsigned int column = sparsity_pattern.column_number(row,c);
-
+
if ((column >= local_col_start) &&
(column < local_col_end))
++row_lengths_in_window[row-local_row_start];
// 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<int> row_entries;
std::vector<PetscScalar> row_values;
row_values.resize (sparsity_pattern.row_length(i), 0.0);
for (unsigned int j=0; j<sparsity_pattern.row_length(i); ++j)
row_entries[j] = sparsity_pattern.column_number (i,j);
-
+
const int int_row = i;
MatSetValues (matrix, 1, &int_row,
sparsity_pattern.row_length(i), &row_entries[0],
}
compress ();
-
+
#else
// first set up the column number
}
colnums_in_window.resize (n_cols+1, -1);
}
-
+
// now copy over the information
// from the sparsity pattern.
{
colnums_in_window[index] = sparsity_pattern.column_number(i,j);
Assert (index == colnums_in_window.size()-1, ExcInternalError());
}
-
+
// then call the petsc function
// that summarily allocates these
// entries:
const std::vector<PetscScalar>
values (sparsity_pattern.max_entries_per_row(),
1.);
-
+
for (unsigned int i=local_row_start; i<local_row_end; ++i)
{
const int petsc_i = i;
&values[0], INSERT_VALUES);
}
}
-
+
compress ();
// set the dummy entries set above
// back to zero
*this = 0;
compress ();
-
+
#endif
// Now we won't insert any
// further entries, so PETSc can
// internally optimize some data
// structures.
-#if (PETSC_VERSION_MAJOR <= 2)
+#if (PETSC_VERSION_MAJOR <= 2)
const int ierr =
MatSetOption (matrix, MAT_NO_NEW_NONZERO_LOCATIONS);
#else
const int ierr =
MatSetOption (matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
#endif
-
+
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
}
// $Id$
// Version: $Name$
//
-// Copyright (C) 2004, 2006, 2008 by the deal.II authors
+// Copyright (C) 2004, 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
{
// destroy the solver object
int ierr = KSPDestroy (ksp);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
// and destroy the solver object if we
// are in old PETSc mode
#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2)
#endif
}
-
-
+
+
SolverBase::SolverBase (SolverControl &cn,
const MPI_Comm &mpi_communicator)
:
mpi_communicator (mpi_communicator)
{}
-
+
SolverBase::~SolverBase ()
{}
-
+
void
SolverBase::solve (const MatrixBase &A,
SLES sles;
ierr = SLESCreate (mpi_communicator, &sles);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
// set the matrices involved. the
// last argument is irrelevant here,
// since we use the solver only once
ierr = SLESSetOperators (sles, A, preconditioner,
SAME_PRECONDITIONER);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
// let derived classes set the solver
// type, and the preconditioning object
// set the type of preconditioner
KSP ksp;
ierr = SLESGetKSP (sles, &ksp);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
set_solver_type (ksp);
-
+
PC pc;
ierr = SLESGetPC (sles, &pc);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
preconditioner.set_preconditioner_type (pc);
-
+
// then a convergence monitor
// function. that function simply checks
// with the solver_control object we have
// in this object for convergence
KSPSetConvergenceTest (ksp, &convergence_test,
reinterpret_cast<void *>(&solver_control));
-
+
// then do the real work: set up solver
// internal data and solve the
// system. this could be joined in one
int iterations = 0;
ierr = SLESSetUp (sles, b, x);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
ierr = SLESSolve (sles, b, x, &iterations);
AssertThrow (ierr == 0, ExcPETScError(ierr));
if (solver_data == 0)
{
solver_data.reset (new SolverData());
-
+
ierr = KSPCreate (mpi_communicator, &solver_data->ksp);
AssertThrow (ierr == 0, ExcPETScError(ierr));
// 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));
// 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<void *>(&solver_control));
#else
#endif
}
-
+
// then do the real work: set up solver
// internal data and solve the
// system. unfortunately, the call
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
ierr = KSPSolve (solver_data->ksp, b, x);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-# endif
+# endif
#endif
// destroy solver object
{
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<SolverControl*>(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<KSPConvergedReason>(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());
}
// return without failure
return 0;
}
-
-
+
+
/* ---------------------- SolverRichardson ------------------------ */
omega (omega)
{}
-
-
+
+
SolverRichardson::SolverRichardson (SolverControl &cn,
const MPI_Comm &mpi_communicator,
const AdditionalData &data)
KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(),
PETSC_DEFAULT, this->solver_control.max_steps()+1);
}
-
+
/* ---------------------- SolverChebychev ------------------------ */
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverCG ------------------------ */
SolverCG::SolverCG (SolverControl &cn,
// solution vector. do so here as well:
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverBiCG ------------------------ */
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverGMRES ------------------------ */
SolverGMRES::AdditionalData::
right_preconditioning (right_preconditioning)
{}
-
-
+
+
SolverGMRES::SolverGMRES (SolverControl &cn,
const MPI_Comm &mpi_communicator,
const AdditionalData &data)
// 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);
// solution vector. do so here as well:
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverBicgstab ------------------------ */
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverCGS ------------------------ */
SolverCGS::SolverCGS (SolverControl &cn,
// solution vector. do so here as well:
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverTFQMR ------------------------ */
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverTCQMR ------------------------ */
SolverTCQMR::SolverTCQMR (SolverControl &cn,
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverCR ------------------------ */
SolverCR::SolverCR (SolverControl &cn,
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverLSQR ------------------------ */
SolverLSQR::SolverLSQR (SolverControl &cn,
// solution vector. do so here as well:
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
}
-
+
/* ---------------------- SolverPreOnly ------------------------ */
// 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);
// $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
{
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
int ierr
= VecGetArray (static_cast<const Vec &>(vector), &ptr);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
const PetscScalar value = *(ptr+index);
-
+
ierr = VecRestoreArray (static_cast<const Vec &>(vector), &ptr);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
return value;
}
else if (dynamic_cast<const PETScWrappers::MPI::Vector *>(&vector) != 0)
PetscScalar *ptr;
ierr = VecGetArray (static_cast<const Vec &>(vector), &ptr);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
const PetscScalar value = *(ptr+index-begin);
-
+
ierr = VecRestoreArray (static_cast<const Vec &>(vector), &ptr);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
return value;
}
else
// exists there?
Assert (false, ExcInternalError());
return -1e20;
- }
+ }
}
-
+
VectorBase::VectorBase ()
:
last_action (LastAction::none)
{}
-
-
+
+
VectorBase::VectorBase (const VectorBase &v)
:
last_action (LastAction::none)
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
-
+
VectorBase::~VectorBase ()
{
const int ierr = VecDestroy (vector);
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
// 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;
{
Assert (size() == v.size(),
ExcDimensionMismatch(size(), v.size()));
-
+
PetscTruth flag;
-
+
const int ierr = VecEqual (vector, v.vector, &flag);
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
-
+
bool
VectorBase::operator != (const VectorBase &v) const
{
ExcDimensionMismatch(size(), v.size()));
PetscTruth flag;
-
+
const int ierr = VecEqual (vector, v.vector, &flag);
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
-
+
unsigned int
VectorBase::size () const
{
{
do_set_add_operation(indices, values, true);
}
-
+
PetscScalar
return result;
}
-
+
void
VectorBase::real_type
VectorBase::norm_sqr () const
{
- const PetscScalar d = l2_norm();
+ const real_type d = l2_norm();
return d*d;
}
// add up remaining elements
while (ptr != start_ptr+size())
sum0 += *ptr++;
-
- mean = (sum0+sum1+sum2+sum3)/size();
+
+ mean = (sum0+sum1+sum2+sum3)/static_cast<PetscScalar>(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
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
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<value_type>::abs(*ptr++), p);
+ sum1 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
+ sum2 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
+ sum3 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
+ }
// add up remaining elements
while (ptr != start_ptr+size())
- sum0 += std::pow(std::fabs(*ptr++), p);
-
+ sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
+
norm = std::pow(sum0+sum1+sum2+sum3,
- static_cast<PetscScalar>(1./p));
+ 1./p);
}
-
+
// restore the representation of the
// vector
ierr = VecRestoreArray (vector, &start_ptr);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
+
return norm;
}
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));
VectorBase::real_type
VectorBase::normalize () const
{
- PetscScalar d;
+ real_type d;
const int ierr = VecNormalize (vector, &d);
AssertThrow (ierr == 0, ExcPETScError(ierr));
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));
VectorBase &
- VectorBase::abs ()
+ VectorBase::abs ()
{
const int ierr = VecAbs (vector);
AssertThrow (ierr == 0, ExcPETScError(ierr));
VectorBase &
- VectorBase::conjugate ()
+ VectorBase::conjugate ()
{
const int ierr = VecConjugate (vector);
AssertThrow (ierr == 0, ExcPETScError(ierr));
VectorBase &
- VectorBase::mult ()
+ VectorBase::mult ()
{
const int ierr = VecPointwiseMult (vector,vector,vector);
AssertThrow (ierr == 0, ExcPETScError(ierr));
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 <typename T>
+ bool is_non_negative (const T &t)
+ {
+ return t >= 0;
+ }
+
+
+
+ template <typename T>
+ bool is_non_negative (const std::complex<T> &)
+ {
+ Assert (false,
+ ExcMessage ("You can't ask a complex value "
+ "whether it is non-negative."))
+ return true;
+ }
+ }
+
+
bool
VectorBase::is_non_negative () const
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)
#else
const int ierr = VecScale (vector, a);
#endif
-
+
AssertThrow (ierr == 0, ExcPETScError(ierr));
return *this;
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)
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;
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;
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)
#else
const int ierr = VecShift (vector, s);
#endif
-
+
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
-
+
void
{
*this += v;
}
-
+
void
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)
#else
const int ierr = VecAXPY (vector, a, v);
#endif
-
+
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
-
+
void
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
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
-
+
void
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)
#else
const int ierr = VecAYPX (vector, s, v);
#endif
-
+
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
-
+
void
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
*this *= s;
add (a,v);
}
-
+
void
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));
}
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
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));
}
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(),
*this *= a;
}
-
+
void
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(),
if (across)
for (unsigned int i=0; i<size(); ++i)
- out << static_cast<double>(val[i]) << ' ';
+ out << val[i] << ' ';
else
for (unsigned int i=0; i<size(); ++i)
- out << static_cast<double>(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
AssertThrow(false, ExcNotImplemented() );
return 0;
}
-
+
void
last_action = LastAction::insert;
}
-
+
}
DEAL_II_NAMESPACE_CLOSE