Flow123d  master-f44eb46
linsys_PERMON.cc
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.cc
15  * @brief Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construction
16  * @author Jakub Sistek
17  */
18 
19 // derived from base linsys
20 #include "la/linsys_PERMON.hh"
21 #include "petscvec.h"
22 #include "petscksp.h"
23 #include "petscmat.h"
24 #include "system/sys_profiler.hh"
25 #include "system/system.hh"
26 
27 
28 //#include <boost/bind.hpp>
29 
30 namespace it = Input::Type;
31 
33  return it::Record("Petsc", "PETSc solver settings.\n It provides interface to various PETSc solvers. The convergence criteria is:\n"
34  "```\n"
35  "norm( res_i ) < max( norm( res_0 ) * r_tol, a_tol )\n"
36  "```\n"
37  "where ```res_i``` is the residuum vector after i-th iteration of the solver and ```res_0``` is the estimate of the norm of the initial residual. "
38  "If the initial guess of the solution is provided (usually only for transient equations) the residual of this estimate is used, "
39  "otherwise the norm of preconditioned RHS is used. "
40  "The default norm is (($L_2$)) norm of preconditioned residual: (($ P^{-1}(Ax-b)$)), usage of other norm may be prescribed using the 'option' key. "
41  "See also PETSc documentation for KSPSetNormType.")
43  .declare_key("r_tol", it::Double(0.0, 1.0), it::Default::read_time("Default value is set by the nonlinear solver or the equation. "
44  "If not, we use the value 1.0e-7."),
45  "Residual tolerance relative to the initial error.")
46  .declare_key("a_tol", it::Double(0.0), it::Default::read_time("Default value is set by the nonlinear solver or the equation. "
47  "If not, we use the value 1.0e-11."),
48  "Absolute residual tolerance.")
49  .declare_key("max_it", it::Integer(0), it::Default::read_time("Default value is set by the nonlinear solver or the equation. "
50  "If not, we use the value 1000."),
51  "Maximum number of outer iterations of the linear solver.")
52  .declare_key("options", it::String(), it::Default("\"\""), "This options is passed to PETSC to create a particular KSP (Krylov space method).\n"
53  "If the string is left empty (by default), the internal default options is used.")
54  .close();
55 }
56 
57 
59 
60 
61 LinSys_PERMON::LinSys_PERMON( const Distribution * rows_ds, const std::string &params)
62  : LinSys_PETSC( rows_ds, params )
63 {
64  matrix_ineq_ = NULL;
65  ineq_ = NULL;
66 }
67 
69  : LinSys_PETSC(other)
70 {
71  MatCopy(other.matrix_ineq_, matrix_ineq_, DIFFERENT_NONZERO_PATTERN);
72  VecCopy(other.ineq_, ineq_);
73 }
74 
75 void LinSys_PERMON::set_inequality(Mat matrix_ineq, Vec ineq)
76 {
77  // TODO ref count?
78  matrix_ineq_ = matrix_ineq;
79  ineq_ = ineq;
80 }
81 
83 {
84 
85  const char *petsc_dflt_opt;
86  int nits;
87 
88  // -mat_no_inode ... inodes are usefull only for
89  // vector problems e.g. MH without Schur complement reduction
90 
91  /* Comment to PETSc options:
92  *
93  * -ksp_diagonal_scale scales the matrix before solution, while -ksp_diagonal_scale_fix just fixes the scaling after solution
94  * -pc_asm_type basic enforces classical Schwartz method, which seems more stable for positive definite systems.
95  * The default 'restricted' probably violates s.p.d. structure, many tests fail.
96  */
97  if (rows_ds_->np() > 1) {
98  // parallel setting
99  if (this->is_positive_definite())
100  petsc_dflt_opt="-ksp_type cg -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_type asm -pc_asm_type basic -pc_asm_overlap 4 -sub_pc_type icc -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
101  //petsc_dflt_opt="-ksp_type bcgs -ksp_diagonal_scale_fix -pc_type asm -pc_asm_overlap 4 -sub_pc_type ilu -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
102  else
103  petsc_dflt_opt="-ksp_type bcgs -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_type asm -pc_asm_overlap 4 -sub_pc_type ilu -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
104 
105  }
106  else {
107  // serial setting
108  if (this->is_positive_definite())
109  petsc_dflt_opt="-ksp_type cg -pc_type icc -pc_factor_levels 3 -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
110  //petsc_dflt_opt="-ksp_type bcgs -pc_type ilu -pc_factor_levels 5 -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
111  else
112  petsc_dflt_opt="-ksp_type bcgs -pc_type ilu -pc_factor_levels 5 -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
113  }
114 
115  if (params_ == "") params_ = petsc_dflt_opt;
116  LogOut().fmt("inserting petsc options: {}\n",params_.c_str());
117 
118  // now takes an optional PetscOptions object as the first argument
119  // value NULL will preserve previous behaviour previous behavior.
120  PetscOptionsInsertString(NULL, params_.c_str()); // overwrites previous options values
121 
122  MatSetOption( matrix_, MAT_USE_INODES, PETSC_FALSE );
123 
124  chkerr(QPCreate(comm_, &system));
125  chkerr(QPSetOperator(system, matrix_));
126  chkerr(QPSetRhs(system, rhs_));
127  chkerr(QPSetInitialVector(system, solution_));
128  if (ineq_) {
129  //chkerr(MatScale(matrix_ineq_,-1.));
130  //chkerr(VecScale(ineq_,-1.));
131  // Bx<=c
132  chkerr(QPSetIneq(system, matrix_ineq_, ineq_));
133  chkerr(QPSetIneq(system, matrix_ineq_, ineq_));
134  chkerr(QPTDualize(system, MAT_INV_MONOLITHIC, MAT_REG_NONE));
135  }
136  // Set runtime options, e.g -qp_chain_view_kkt
137  chkerr(QPSetFromOptions(system));
138 
139 
140  chkerr(QPSCreate(comm_, &solver));
141  chkerr(QPSSetQP(solver, system));
142 
143  // TODO take care of tolerances - shall we support both input file and command line petsc setting
144  chkerr(QPSSetTolerances(solver, r_tol_, a_tol_, d_tol_, PETSC_DEFAULT));
145  chkerr(QPSSetTolerances(solver, r_tol_, a_tol_, d_tol_, max_it_));
146  chkerr(QPSSetFromOptions(solver));
147 
148  {
149  START_TIMER("PERMON linear solver");
150  START_TIMER("PERMON linear iteration");
151  chkerr(QPSSolve(solver));
152  QPSGetConvergedReason(solver,&reason);
153  QPSGetIterationNumber(solver,&nits);
154  ADD_CALLS(nits);
155  }
156  // substitute by PETSc call for residual
157  //VecNorm(rhs_, NORM_2, &residual_norm_);
158 
159  LogOut().fmt("convergence reason {}, number of iterations is {}\n", reason, nits);
160 
161  // get residual norm
162  //KSPGetResidualNorm(system, &solution_precision_);
163 
164  // TODO: I do not understand this
165  //Profiler::instance()->set_timer_subframes("SOLVING MH SYSTEM", nits);
166 
167  chkerr(QPSDestroy(&solver));
168  chkerr(QPDestroy(&system));
169 
170  return LinSys::SolveInfo(static_cast<int>(reason), static_cast<int>(nits));
171 
172 }
173 
174 void LinSys_PERMON::view(string text )
175 {
176  FilePath matFileName(text + "_flow123d_matrix.m",FilePath::FileType::output_file);
177  FilePath rhsFileName(text + "_flow123d_rhs.m",FilePath::FileType::output_file);
178  FilePath solFileName(text + "_flow123d_sol.m",FilePath::FileType::output_file);
179  FilePath mat_ineqFileName(text + "_flow123d_matrix_ineq.m",FilePath::FileType::output_file);
180  FilePath ineqFileName(text + "_flow123d_ineq.m",FilePath::FileType::output_file);
181 
182  PetscViewer myViewer;
183 
184  if ( matrix_ != NULL ) {
185  PetscViewerASCIIOpen(comm_,((string)matFileName).c_str(),&myViewer);
186  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
187  MatView( matrix_, myViewer );
188  PetscViewerDestroy(&myViewer);
189  }
190  else
191  WarningOut() << "PetscViewer: the matrix of LinSys is not set.\n";
192 
193  if ( rhs_ != NULL ) {
194  PetscViewerASCIIOpen(comm_,((string)rhsFileName).c_str(),&myViewer);
195  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
196  VecView( rhs_, myViewer );
197  PetscViewerDestroy(&myViewer);
198  }
199  else
200  WarningOut() << "PetscViewer: the rhs of LinSys is not set.\n";
201 
202  if ( solution_ != NULL ) {
203  PetscViewerASCIIOpen(comm_,((string)solFileName).c_str(),&myViewer);
204  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
205  VecView( solution_, myViewer );
206  PetscViewerDestroy(&myViewer);
207  }
208  else
209  WarningOut() << "PetscViewer: the solution of LinSys is not set.\n";
210  if ( matrix_ineq_ != NULL ) {
211  PetscViewerASCIIOpen(comm_,((string)mat_ineqFileName).c_str(),&myViewer);
212  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
213  MatView( matrix_ineq_, myViewer );
214  PetscViewerDestroy(&myViewer);
215  }
216  else
217  WarningOut() << "PetscViewer: the inequality matrix of LinSys is not set.\n";
218  if ( ineq_ != NULL ) {
219  PetscViewerASCIIOpen(comm_,((string)ineqFileName).c_str(),&myViewer);
220  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
221  VecView( ineq_, myViewer );
222  PetscViewerDestroy(&myViewer);
223  }
224  else
225  WarningOut() << "PetscViewer: the inequality vector of LinSys is not set.\n";
226 }
227 
229 {
230  // TODO cleanup ineq objects
231 }
232 
233 
234 
235 
237 {
238  return solution_precision_;
239 }
240 
241 
243 {
244  // TODO return ||A*x - b + B'*lambda||?
245  MatMult(matrix_, solution_, residual_);
246  VecAXPY(residual_,-1.0, rhs_);
247  double residual_norm;
248  VecNorm(residual_, NORM_2, &residual_norm);
249  return residual_norm;
250 }
unsigned int np() const
get num of processors
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default read_time(const std::string &description)
The factory function to make an default value that will be specified at the time when a key will be r...
Definition: type_record.hh:97
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Class for declaration of the integral input data.
Definition: type_base.hh:483
Record type proxy class.
Definition: type_record.hh:182
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Class for declaration of the input data that are in string format.
Definition: type_base.hh:582
static const int registrar
Registrar of class to factory.
void view(string text="") override
LinSys_PERMON(const Distribution *rows_ds, const std::string &params="")
Mat matrix_ineq_
PETSc matrix of inequality constraint.
double get_solution_precision() override
LinSys_PETSC::SolveInfo solve() override
static const Input::Type::Record & get_input_type()
Vec ineq_
PETSc vector of inequality constraint.
double compute_residual() override
void set_inequality(Mat matrix_ineq, Vec ineq)
Mat matrix_
Petsc matrix of the problem.
double solution_precision_
std::string params_
command-line-like options for the PETSc solver
KSPConvergedReason reason
Vec rhs_
PETSc vector constructed with vx array.
unsigned int max_it_
maximum number of iterations of linear solver
Definition: linsys.hh:668
double d_tol_
tolerance for divergence of linear solver
Definition: linsys.hh:667
MPI_Comm comm_
Definition: linsys.hh:670
bool is_positive_definite()
Definition: linsys.hh:597
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
Vec solution_
PETSc vector constructed with vb array.
Definition: linsys.hh:686
double r_tol_
relative tolerance of linear solver
Definition: linsys.hh:665
static Input::Type::Abstract & get_input_type()
Definition: linsys.cc:29
PERMON QP solvers and FETI.
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
#define LogOut()
Macro defining 'log' record of log.
Definition: logger.hh:281
#define ADD_CALLS(n_calls)
Increase number of calls in actual timer.
#define START_TIMER(tag)
Starts a timer with specified tag.
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:142