Flow123d  master-f44eb46
linsys_PETSC.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file linsys_PETSC.hh
15  * @brief Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construction
16  * @author Jakub Sistek
17  */
18 
19 #ifndef LA_LINSYS_PETSC_HH_
20 #define LA_LINSYS_PETSC_HH_
21 
22 #include <functional> // for unary_function
23 #include <string> // for string
24 #include <vector> // for vector
25 #include "la/linsys.hh" // for LinSys
26 #include "petscksp.h" // for KSP, KSPConvergedReason, _p_KSP
27 #include "petscmat.h" // for Mat, MatCopy, MatZeroEntries, MatAssemblyType
28 #include "petscmath.h" // for PetscScalar
29 #include "petscsys.h" // for PetscErrorCode, PETSC_NULL
30 #include "petscvec.h" // for Vec, _p_Vec, VecCopy, VecSet
31 
32 class Distribution;
33 namespace Input {
34  class Record;
35  namespace Type {
36  class Record;
37  }
38 }
39 namespace la {
40  class BddcmlWrapper;
41 }
42 
43 class LinSys_PETSC : public LinSys
44 {
45 
46 public:
48 
49  static const Input::Type::Record & get_input_type();
50 
51  LinSys_PETSC(const Distribution * rows_ds, const std::string &params = "");
52 
53  /**
54  * Copy constructor.
55  */
56  LinSys_PETSC( LinSys_PETSC &other );
57 
58 
59  void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) override;
60 
61  /**
62  * Returns whole Distribution class for distribution of the solution.
63  */
64  inline const Distribution* get_ds( )
65  {
66  return rows_ds_;
67  }
68 
69  const Mat *get_matrix() override
70  {
71  return &matrix_;
72  }
73 
74  const Vec *get_rhs() override
75  {
76  return &rhs_;
77  }
78 
79  PetscErrorCode set_matrix(Mat &matrix, MatStructure str) override
80  {
81  matrix_changed_ = true;
82  return MatCopy(matrix, matrix_, str);
83  }
84 
85  PetscErrorCode set_rhs(Vec &rhs) override
86  {
87  rhs_changed_ = true;
88  return VecCopy(rhs, rhs_);
89  }
90 
91  PetscErrorCode mat_zero_entries() override
92  {
93  matrix_changed_ = true;
94  constraints_.clear();
95  return MatZeroEntries(matrix_);
96  }
97 
98  PetscErrorCode rhs_zero_entries() override
99  {
100  rhs_changed_ = true;
101  return VecSet(rhs_, 0);
102  }
103 
104  void start_allocation() override;
105 
106  void start_add_assembly() override;
107 
108  void start_insert_assembly() override;
109 
110  void mat_set_values( int nrow, int *rows, int ncol, int *cols, double *vals ) override;
111 
112  void rhs_set_values( int nrow, int *rows, double *vals ) override;
113 
114  void preallocate_values(int nrow,int *rows,int ncol,int *cols);
115 
116  void preallocate_matrix();
117 
118  void finish_assembly() override;
119 
120  void finish_assembly( MatAssemblyType assembly_type );
121 
122  void apply_constrains( double scalar = 1. ) override;
123 
124  void set_initial_guess_nonzero(bool set_nonzero = true);
125 
126  LinSys::SolveInfo solve() override;
127 
128  /**
129  * Returns information on absolute solver accuracy
130  */
131  inline double get_absolute_accuracy() override {
132  return a_tol_;
133  };
134 
135  void view(string text="") override;
136 
137  /**
138  * Sets specific parameters of LinSys_PETSC defined by user in input file and used to calculate
139  */
140  void set_from_input(const Input::Record in_rec) override;
141 
142  double get_solution_precision() override;
143 
144  double compute_residual() override;
145 
146  ~LinSys_PETSC( );
147 
148 private:
149  /// Registrar of class to factory
150  static const int registrar;
151 
152  // make a pointer to the data array out of a std::vector
153  template<typename T>
155  {
156  if ( array.size() ) return &(array[0]);
157  return PETSC_NULL;
158  }
159 
160  // PetscScalar to double casting functor
161  struct PetscScalar2Double_ : public std::unary_function< PetscScalar, double >
162  {
163  double operator()( PetscScalar arg )
164  {
165  return static_cast<double>( arg );
166  }
167  };
168 
169 protected:
170 
171  std::string params_; //!< command-line-like options for the PETSc solver
172 
173  bool init_guess_nonzero; //!< flag for starting from nonzero guess
174 
175  Mat matrix_; //!< Petsc matrix of the problem.
176  Vec rhs_; //!< PETSc vector constructed with vx array.
178 
179  double *v_rhs_; //!< local RHS array pointing to Vec rhs_
180 
181  Vec on_vec_; //!< Vectors for counting non-zero entries in diagonal block.
182  Vec off_vec_; //!< Vectors for counting non-zero entries in off-diagonal block.
183 
184 
185  double solution_precision_; // precision of KSP system solver
186 
187  KSP system;
188  KSPConvergedReason reason;
189 
190 
191 };
192 
193 #endif /* LA_LINSYS_PETSC_HH_ */
unsigned int size() const
Definition: armor.hh:728
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
Record type proxy class.
Definition: type_record.hh:182
void finish_assembly() override
void start_add_assembly() override
void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals) override
PetscErrorCode set_rhs(Vec &rhs) override
Definition: linsys_PETSC.hh:85
void apply_constrains(double scalar=1.) override
double compute_residual() override
const Distribution * get_ds()
Definition: linsys_PETSC.hh:64
Mat matrix_
Petsc matrix of the problem.
void preallocate_matrix()
LinSys::SolveInfo solve() override
double solution_precision_
void preallocate_values(int nrow, int *rows, int ncol, int *cols)
bool init_guess_nonzero
flag for starting from nonzero guess
void start_insert_assembly() override
void set_from_input(const Input::Record in_rec) override
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
void set_initial_guess_nonzero(bool set_nonzero=true)
LinSys FactoryBaseType
Definition: linsys_PETSC.hh:47
std::string params_
command-line-like options for the PETSc solver
LinSys_PETSC(const Distribution *rows_ds, const std::string &params="")
Definition: linsys_PETSC.cc:64
KSPConvergedReason reason
PetscErrorCode set_matrix(Mat &matrix, MatStructure str) override
Definition: linsys_PETSC.hh:79
PetscErrorCode mat_zero_entries() override
Definition: linsys_PETSC.hh:91
PetscErrorCode rhs_zero_entries() override
Definition: linsys_PETSC.hh:98
static const int registrar
Registrar of class to factory.
double get_absolute_accuracy() override
Vec rhs_
PETSc vector constructed with vx array.
Vec off_vec_
Vectors for counting non-zero entries in off-diagonal block.
Vec on_vec_
Vectors for counting non-zero entries in diagonal block.
const Mat * get_matrix() override
Definition: linsys_PETSC.hh:69
T * makePetscPointer_(std::vector< T > &array)
void view(string text="") override
void rhs_set_values(int nrow, int *rows, double *vals) override
double * v_rhs_
local RHS array pointing to Vec rhs_
void start_allocation() override
double get_solution_precision() override
const Vec * get_rhs() override
Definition: linsys_PETSC.hh:74
void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) override
Definition: linsys_PETSC.cc:93
ConstraintVec_ constraints_
Definition: linsys.hh:693
double a_tol_
absolute tolerance of linear solver
Definition: linsys.hh:666
const Distribution * rows_ds_
final distribution of rows of MH matrix
Definition: linsys.hh:676
bool matrix_changed_
true if the matrix was changed since the last solve
Definition: linsys.hh:683
bool rhs_changed_
true if the right hand side was changed since the last solve
Definition: linsys.hh:684
Wrappers for linear systems based on MPIAIJ and MATIS format.
Array< double > array
Definition: armor.hh:890
Abstract linear system class.
Definition: balance.hh:40
internal::NamedArg< char > arg(StringRef name, const T &arg)
Definition: format.h:3291
double operator()(PetscScalar arg)