75 #include "petscmath.h"
86 namespace Input {
namespace Type {
class Abstract; } }
133 int lsizeInt =
static_cast<int>( rows_ds->
lsize() );
136 size_ =
static_cast<unsigned>( sizeInt );
156 ASSERT_PERMANENT(
false ).error(
"Using copy constructor of LinSys is not allowed!");
189 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function get_matrix is not implemented for linsys type\n.");
205 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function get_rhs is not implemented for linsys type\n.");
228 virtual void set_tolerances(
double r_tol,
double a_tol,
double d_tol,
unsigned int max_it) = 0;
248 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function set_matrix is not implemented for linsys type \n.");
257 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function set_rhs is not implemented for linsys type \n.");
266 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function mat_zero_entries is not implemented for linsys type \n.");
275 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function vec_zero_entries is not implemented for linsys type \n.");
297 VecRestoreArray(
solution_, &out_array );
335 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function start_allocation is not implemented for linsys type \n.");
343 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function start_add_assembly is not implemented for linsys type \n.");
351 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function start_insert_assembly is not implemented for linsys type \n.");
386 inline void set_values(
int nrow,
int *
rows,
int ncol,
int *
cols,PetscScalar *mat_vals, PetscScalar *rhs_vals)
421 for (
uint i=0; i<m; i++)
422 row_dofs[i]= local_to_global_map[local.
row_dofs[i]];
423 for (
uint i=0; i<n; i++)
424 col_dofs[i]= local_to_global_map[local.
col_dofs[i]];
459 bool negative_row =
false;
460 bool negative_col =
false;
462 for(
unsigned int l_row = 0; l_row < row_dofs.size(); l_row++)
463 if (row_dofs[l_row] < 0) {
465 tmp.col(l_row).zeros();
469 for(
unsigned int l_col = 0; l_col < col_dofs.size(); l_col++)
470 if (col_dofs[l_col] < 0) {
471 tmp_rhs -= matrix.col(l_col) * col_solution[l_col];
472 tmp.row(l_col).zeros();
477 if (negative_row && negative_col) {
479 for(
unsigned int l_row = 0; l_row < row_dofs.size(); l_row++)
480 if (row_dofs[l_row] < 0)
481 for(
unsigned int l_col = 0; l_col < col_dofs.size(); l_col++)
482 if (col_dofs[l_col] < 0 && row_dofs[l_row] == col_dofs[l_col]) {
483 double new_diagonal = fabs(matrix.at(l_row,l_col));
484 if (new_diagonal == 0.0) {
485 if (matrix.is_square()) {
486 new_diagonal = arma::sum( abs(matrix.diag())) / matrix.n_rows;
488 new_diagonal = arma::accu( abs(matrix) ) / matrix.n_elem;
491 tmp.at(l_col, l_row) = new_diagonal;
492 tmp_rhs(l_row) = new_diagonal * row_solution[l_row];
498 for(
int &row : row_dofs) row=abs(row);
501 for(
int &col : col_dofs) col=abs(col);
504 mat_set_values(row_dofs.size(),
const_cast<int *
>(&(row_dofs[0])),
505 col_dofs.size(),
const_cast<int *
>(&(col_dofs[0])), tmp.memptr() );
506 rhs_set_values(row_dofs.size(),
const_cast<int *
>(&(row_dofs[0])), tmp_rhs.memptr() );
635 ASSERT_PERMANENT(
false )(
typeid(*this).name()).error(
"Function view is not implemented for linsys type %s \n.");
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
unsigned int lsize(int proc) const
get local size
const Vec & get_solution()
static constexpr double default_r_tol_
unsigned int max_it_
maximum number of iterations of linear solver
virtual const Vec * get_rhs()
void set_negative_definite(bool flag=true)
void set_local_system(LocalSystem &local, const std::vector< LongIdx > &local_to_global_map)
void set_solution(Vec sol_vec)
void set_values(std::vector< int > &row_dofs, std::vector< int > &col_dofs, const arma::mat &matrix, const arma::vec &rhs, const arma::vec &row_solution, const arma::vec &col_solution)
static constexpr double default_a_tol_
double get_residual_norm()
bool spd_via_symmetric_general_
double get_relative_accuracy()
double d_tol_
tolerance for divergence of linear solver
void set_matrix_changed()
std::vector< double > globalSolution_
global solution in numbering for linear system
bool is_spd_via_symmetric_general()
double * get_solution_array()
bool is_negative_definite()
virtual void set_from_input(const Input::Record in_rec)
virtual void rhs_set_values(int nrow, int *rows, double *vals)=0
bool is_positive_definite()
virtual double get_solution_precision()=0
virtual void start_add_assembly()
virtual void finish_assembly()=0
static constexpr double default_d_tol_
void rhs_set_value(int row, double val)
virtual const Mat * get_matrix()
virtual void apply_constrains(double scalar)=0
void set_solution(double *sol_array)
virtual void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it)=0
ConstraintVec_ constraints_
unsigned size_
global number of matrix rows, i.e. problem size
virtual void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals)=0
static constexpr unsigned int default_max_it_
virtual double compute_residual()=0
bool own_solution_
Indicates if the solution array has been allocated by this class.
void set_local_system(LocalSystem &local)
void set_symmetric(bool flag=true)
virtual PetscErrorCode rhs_zero_entries()
double a_tol_
absolute tolerance of linear solver
void mat_set_value(int row, int col, double val)
void add_constraint(int row, double value)
std::pair< unsigned, double > Constraint_
SetValuesMode status_
Set value status of the linear system.
bool own_vec_
Indicates if the solution vector has been allocated by this class.
virtual void start_insert_assembly()
virtual void start_allocation()
virtual PetscErrorCode set_matrix(Mat &, MatStructure)
const Distribution * rows_ds_
final distribution of rows of MH matrix
Vec solution_
PETSc vector constructed with vb array.
virtual SolveInfo solve()=0
virtual void view(string)
void set_values(int nrow, int *rows, int ncol, int *cols, PetscScalar *mat_vals, PetscScalar *rhs_vals)
Set values in the system matrix and values in the right-hand side vector on corresponding rows.
double * v_solution_
local solution array pointing into Vec solution_
double residual_norm_
local solution array pointing into Vec solution_
virtual PetscErrorCode set_rhs(Vec &)
void set_spd_via_symmetric_general(bool flag=true)
LinSys(const Distribution *rows_ds)
void set_positive_definite(bool flag=true)
double r_tol_
relative tolerance of linear solver
virtual double get_absolute_accuracy()
bool matrix_changed_
true if the matrix was changed since the last solve
const unsigned lsize_
local number of matrix rows (non-overlapping division of rows)
virtual PetscErrorCode mat_zero_entries()
std::vector< Constraint_ > ConstraintVec_
static Input::Type::Abstract & get_input_type()
bool rhs_changed_
true if the right hand side was changed since the last solve
arma::vec rhs
local system RHS
arma::mat matrix
local system matrix
Support classes for parallel programing.
static constexpr bool value
#define MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm)
ArmaMat< double, N, M > mat
SolveInfo(int cr, int nits)
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.