HermesCommon  3.0
petsc_solver.cpp
Go to the documentation of this file.
1 // This file is part of HermesCommon
2 //
3 // Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
4 // Email: hpfem-group@unr.edu, home page: http://www.hpfem.org/.
5 //
6 // Hermes2D is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published
8 // by the Free Software Foundation; either version 2 of the License,
9 // or (at your option) any later version.
10 //
11 // Hermes2D is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with Hermes2D; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 #include "config.h"
23 #ifdef WITH_PETSC
24 #include "petsc_solver.h"
25 #include "callstack.h"
26 #include "common.h"
27 #include "util/memory_handling.h"
28 
30 
31 namespace Hermes
32 {
33  namespace Algebra
34  {
35  static int num_petsc_objects = 0;
36 
37 #ifdef PETSC_USE_COMPLEX
38  inline void vec_get_value(Vec x, PetscInt ni, const PetscInt ix[], std::complex<double> y[])
39  {
40  VecGetValues(x, ni, ix, y);
41  }
42 
43  void vec_get_value(Vec x, PetscInt ni, const PetscInt ix[], double y[])
44  {
45  PetscScalar *py = malloc_with_check<PetscScalar>(ni);
46  VecGetValues(x, ni, ix, py);
47  for (int i = 0; i < ni; i++)y[i] = py[i].real();
48  free_with_check(py);
49  }
50 #else
51  inline void vec_get_value(Vec x, PetscInt ni, const PetscInt ix[], double y[])
52  {
53  VecGetValues(x, ni, ix, y);
54  }
55  inline void vec_get_value(Vec x, PetscInt ni, const PetscInt ix[], std::complex<double> y[])
56  {
57  throw(Exceptions::Exception("PETSc with complex numbers support required."));
58  }
59 #endif
60 
61  int remove_petsc_object()
62  {
63  PetscBool petsc_initialized, petsc_finalized;
64  int ierr = PetscFinalized(&petsc_finalized); CHKERRQ(ierr);
65  ierr = PetscInitialized(&petsc_initialized); CHKERRQ(ierr);
66  if (petsc_finalized == PETSC_TRUE || petsc_initialized == PETSC_FALSE)
67  // This should never happen here.
68  return -1;
69 
70  if (--num_petsc_objects == 0)
71  {
72  int ierr = PetscFinalize();
73  CHKERRQ(ierr);
74  //FIXME this->info("PETSc finalized. No more PETSc usage allowed until application restart.");
75  }
76  }
77 
78  int add_petsc_object()
79  {
80  int ierr;
81  PetscBool petsc_initialized, petsc_finalized;
82  ierr = PetscFinalized(&petsc_finalized); CHKERRQ(ierr);
83 
84  if (petsc_finalized == PETSC_TRUE)
85  throw Hermes::Exceptions::Exception("PETSc cannot be used once it has been finalized. You must restart the application.");
86 
87  ierr = PetscInitialized(&petsc_initialized); CHKERRQ(ierr);
88 
89  if (petsc_initialized != PETSC_TRUE)
90  {
91  ierr = PetscInitializeNoArguments();
92  CHKERRQ(ierr);
93  }
94 
95  num_petsc_objects++;
96  }
97 
98  template<typename Scalar>
99  PetscMatrix<Scalar>::PetscMatrix()
100  {
101  inited = false;
102  add_petsc_object();
103  }
104 
105  template<typename Scalar>
106  PetscMatrix<Scalar>::~PetscMatrix()
107  {
108  free();
109  remove_petsc_object();
110  }
111 
112  template<typename Scalar>
113  void PetscMatrix<Scalar>::alloc()
114  {
115  assert(this->pages != nullptr);
116 
117  // calc nnz
118  int *nnz_array = malloc_with_check<int>(this->size);
119 
120  // fill in nnz_array
121  int aisize = this->get_num_indices();
122  int *ai = malloc_with_check<int>(aisize);
123 
124  // sort the indices and remove duplicities, insert into ai
125  int pos = 0;
126  for (unsigned int i = 0; i < this->size; i++)
127  {
128  nnz_array[i] = this->sort_and_store_indices(this->pages[i], ai + pos, ai + aisize);
129  pos += nnz_array[i];
130  }
131  // stote the number of nonzeros
132  nnz = pos;
133  free_with_check(this->pages);
134  free_with_check(ai);
135 
136  //
137  MatCreateSeqAIJ(PETSC_COMM_SELF, this->size, this->size, 0, nnz_array, &matrix);
138  // MatSetOption(matrix, MAT_ROW_ORIENTED);
139  // MatSetOption(matrix, MAT_ROWS_SORTED);
140  free_with_check(nnz_array);
141 
142  inited = true;
143  }
144 
145  template<typename Scalar>
146  void PetscMatrix<Scalar>::free()
147  {
148  if (inited) MatDestroy(&matrix);
149  inited = false;
150  }
151 
152  template<typename Scalar>
153  void PetscMatrix<Scalar>::finish()
154  {
155  MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY);
156  MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY);
157  }
158 
159  template<>
160  double PetscMatrix<double>::get(unsigned int m, unsigned int n) const
161  {
162  PetscScalar pv;
163  MatGetValues(matrix, 1, (PetscInt*)&m, 1, (PetscInt*)&n, &pv);
164  return PetscRealPart(pv);
165  }
166 
167 #ifdef PETSC_USE_COMPLEX
168  template<>
169  std::complex<double> PetscMatrix<std::complex<double> >::get(unsigned int m, unsigned int n) const
170  {
171  std::complex<double> v = 0.0;
172  MatGetValues(matrix, 1, (PetscInt*)&m, 1, (PetscInt*)&n, &v);
173  return v;
174  }
175 #else
176  template<>
177  std::complex<double> PetscMatrix<std::complex<double> >::get(unsigned int m, unsigned int n) const
178  {
179  throw(Exceptions::Exception("PETSc with complex numbers support required."));
180  }
181 #endif
182 
183  template<typename Scalar>
184  void PetscMatrix<Scalar>::zero()
185  {
186  MatZeroEntries(matrix);
187  }
188 
189 #ifdef PETSC_USE_COMPLEX
190  inline PetscScalar to_petsc(double x)
191  {
192  return std::complex<double>(x, 0);
193  }
194 
195  inline PetscScalar to_petsc(std::complex<double> x)
196  {
197  return x;
198  }
199 #else
200  inline PetscScalar to_petsc(double x)
201  {
202  return x;
203  }
204  inline PetscScalar to_petsc(std::complex<double> x)
205  {
206  throw(Exceptions::Exception("PETSc with complex numbers support required."));
207  }
208 #endif
209 
210  template<typename Scalar>
211  void PetscMatrix<Scalar>::add(unsigned int m, unsigned int n, Scalar v)
212  {
213  if (v != 0.0)
214  { // ignore zero values.
215 #pragma omp critical (PetscMatrix_add)
216  MatSetValue(matrix, (PetscInt)m, (PetscInt)n, to_petsc(v), ADD_VALUES);
217  }
218  }
219 
220  template<typename Scalar>
221  void PetscMatrix<Scalar>::export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char* number_format)
222  {
223  throw Exceptions::MethodNotImplementedException("PetscVector<double>::export_to_file");
224  /*
225  switch (fmt)
226  {
227  case DF_MATLAB_SPARSE: //only to stdout
228  PetscViewer viewer = PETSC_VIEWER_STDOUT_SELF;
229  PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
230  MatView(matrix, viewer);
231  return true;
232  }
233  */
234  }
235 
236  template<typename Scalar>
237  unsigned int PetscMatrix<Scalar>::get_nnz() const
238  {
239  return nnz;
240  }
241 
242  template<typename Scalar>
243  double PetscMatrix<Scalar>::get_fill_in() const
244  {
245  return (double)nnz / ((double)this->size*this->size);
246  }
247 
248  template<typename Scalar>
249  void PetscMatrix<Scalar>::add_petsc_matrix(PetscMatrix<Scalar>* mat)
250  {
251  //matrix = 1*mat + matrix (matrix and mat have different nonzero structure)
252  MatAXPY(matrix, 1, mat->matrix, DIFFERENT_NONZERO_PATTERN);
253  }
254 
255  template<typename Scalar>
256  void PetscMatrix<Scalar>::add_sparse_matrix(SparseMatrix<Scalar>* mat)
257  {
258  PetscMatrix<Scalar>* mat_petsc = (PetscMatrix<Scalar>*)mat;
259  if (mat_petsc)
260  this->add_petsc_matrix(mat_petsc);
261  else
263  }
264 
265  // Multiplies matrix with a Scalar.
266 
267  template<typename Scalar>
268  void PetscMatrix<Scalar>::multiply_with_Scalar(Scalar value)
269  {
270  MatScale(matrix, to_petsc(value));
271  }
272 
273  template<typename Scalar>
274  void PetscMatrix<Scalar>::create(unsigned int size, unsigned int nnz, int* ap, int* ai, Scalar* ax)
275  {
276  this->size = size;
277  this->nnz = nnz;
278  PetscScalar* pax = malloc_with_check<PetscScalar>(nnz, this);
279  for (unsigned i = 0; i < nnz; i++)
280  pax[i] = to_petsc(ax[i]);
281  MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, size, size, ap, ai, pax, &matrix);
282  delete pax;
283  }
284 
285  template<typename Scalar>
286  PetscMatrix<Scalar>* PetscMatrix<Scalar>::duplicate() const
287  {
288  PetscMatrix<Scalar>*ptscmatrix = new PetscMatrix<Scalar>();
289  MatDuplicate(matrix, MAT_COPY_VALUES, &(ptscmatrix->matrix));
290  ptscmatrix->size = this->size;
291  ptscmatrix->nnz = nnz;
292  return ptscmatrix;
293  }
294 
295  template<typename Scalar>
296  PetscVector<Scalar>::PetscVector()
297  {
298  inited = false;
299  add_petsc_object();
300  }
301 
302  template<typename Scalar>
303  PetscVector<Scalar>::~PetscVector()
304  {
305  free();
306  remove_petsc_object();
307  }
308 
309  template<typename Scalar>
310  void PetscVector<Scalar>::alloc(unsigned int n)
311  {
312  free();
313  this->size = n;
314  VecCreateSeq(PETSC_COMM_SELF, this->size, &vec);
315  inited = true;
316  }
317 
318  template<typename Scalar>
319  void PetscVector<Scalar>::free()
320  {
321  if (inited) VecDestroy(&vec);
322  inited = false;
323  }
324 
325  template<typename Scalar>
326  void PetscVector<Scalar>::finish()
327  {
328  VecAssemblyBegin(vec);
329  VecAssemblyEnd(vec);
330  }
331 
332  template<>
333  double PetscVector<double>::get(unsigned int idx) const
334  {
335  PetscScalar py;
336  VecGetValues(vec, 1, (PetscInt*)&idx, &py);
337  return PetscRealPart(py);
338  }
339 #ifdef PETSC_USE_COMPLEX
340  template<>
341  std::complex<double> PetscVector<std::complex<double> >::get(unsigned int idx) const
342  {
343  std::complex<double> y = 0;
344  VecGetValues(vec, 1, (PetscInt*)&idx, &y);
345  return y;
346  }
347 #else
348  template<>
349  std::complex<double> PetscVector<std::complex<double> >::get(unsigned int idx) const
350  {
351  throw(Exceptions::Exception("PETSc with complex numbers support required."));
352  }
353 #endif
354 
355  template<typename Scalar>
356  void PetscVector<Scalar>::extract(Scalar *v) const
357  {
358  int *idx = malloc_with_check<int>(this->size);
359  for (unsigned int i = 0; i < this->size; i++) idx[i] = i;
360  vec_get_value(vec, this->size, idx, v);
361  free_with_check(idx);
362  }
363 
364  template<typename Scalar>
365  void PetscVector<Scalar>::zero()
366  {
367  VecZeroEntries(vec);
368  }
369 
370  template<typename Scalar>
371  Vector<Scalar>* PetscVector<Scalar>::change_sign()
372  {
373  PetscScalar* y = malloc_with_check<PetscScalar>(this->size);
374  int *idx = malloc_with_check<int>(this->size);
375  for (unsigned int i = 0; i < this->size; i++) idx[i] = i;
376  VecGetValues(vec, this->size, idx, y);
377  for (unsigned int i = 0; i < this->size; i++) y[i] *= -1.;
378  VecSetValues(vec, this->size, idx, y, INSERT_VALUES);
379  free_with_check(y);
380  free_with_check(idx);
381  return this;
382  }
383 
384  template<typename Scalar>
385  void PetscVector<Scalar>::set(unsigned int idx, Scalar y)
386  {
387  VecSetValue(vec, idx, to_petsc(y), INSERT_VALUES);
388  }
389 
390  template<typename Scalar>
391  void PetscVector<Scalar>::add(unsigned int idx, Scalar y)
392  {
393 #pragma omp critical (PetscVector_add)
394  VecSetValue(vec, idx, to_petsc(y), ADD_VALUES);
395  }
396 
397  template<typename Scalar>
398  void PetscVector<Scalar>::add(unsigned int n, unsigned int *idx, Scalar *y)
399  {
400  PetscScalar py;
401  for (unsigned int i = 0; i < n; i++)
402  {
403  VecSetValue(vec, idx[i], to_petsc(y[i]), ADD_VALUES);
404  }
405  }
406 
407  template<typename Scalar>
408  Vector<Scalar>* PetscVector<Scalar>::add_vector(Vector<Scalar>* vec)
409  {
410  assert(this->get_size() == vec->get_size());
411  for (unsigned int i = 0; i < this->get_size(); i++)
412  this->add(i, vec->get(i));
413  return this;
414  }
415 
416  template<typename Scalar>
417  Vector<Scalar>* PetscVector<Scalar>::add_vector(Scalar* vec)
418  {
419  for (unsigned int i = 0; i < this->get_size(); i++)
420  this->add(i, vec[i]);
421  return this;
422  }
423 
424  template<typename Scalar>
425  void PetscVector<Scalar>::export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char* number_format)
426  {
427  throw Exceptions::MethodNotImplementedException("PetscVector<double>::export_to_file");
428  /*
429  switch (fmt)
430  {
431  case DF_MATLAB_SPARSE: //only to stdout
432  PetscViewer viewer = PETSC_VIEWER_STDOUT_SELF;
433  PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
434  VecView(vec, viewer);
435  return true;
436  }
437  */
438  }
439 
440  template class HERMES_API PetscMatrix < double > ;
441  template class HERMES_API PetscMatrix < std::complex<double> > ;
442  template class HERMES_API PetscVector < double > ;
443  template class HERMES_API PetscVector < std::complex<double> > ;
444  }
445  namespace Solvers
446  {
447  template<typename Scalar>
448  PetscLinearMatrixSolver<Scalar>::PetscLinearMatrixSolver(PetscMatrix<Scalar> *mat, PetscVector<Scalar> *rhs)
449  : DirectSolver<Scalar>(mat, rhs), m(mat), rhs(rhs)
450  {
451  add_petsc_object();
452  }
453 
454  template<typename Scalar>
455  PetscLinearMatrixSolver<Scalar>::~PetscLinearMatrixSolver()
456  {
457  remove_petsc_object();
458  }
459 
460  template<typename Scalar>
461  int PetscLinearMatrixSolver<Scalar>::get_matrix_size()
462  {
463  return m->get_size();
464  }
465 
466  template<typename Scalar>
467  void PetscLinearMatrixSolver<Scalar>::solve()
468  {
469  assert(m != nullptr);
470  assert(rhs != nullptr);
471 
472  PetscErrorCode ec;
473  KSP ksp;
474  Vec x;
475 
476  this->tick();
477 
478  KSPCreate(PETSC_COMM_WORLD, &ksp);
479 
480  KSPSetOperators(ksp, m->matrix, m->matrix, DIFFERENT_NONZERO_PATTERN);
481  KSPSetFromOptions(ksp);
482  VecDuplicate(rhs->vec, &x);
483 
484  ec = KSPSolve(ksp, rhs->vec, x);
485  //FIXME if (ec) return false;
486 
487  this->tick();
488  this->time = this->accumulated();
489 
490  // allocate memory for solution vector
491  free_with_check(this->sln);
492  this->sln = malloc_with_check<Scalar>(m->get_size());
493  memset(this->sln, 0, m->get_size() * sizeof(Scalar));
494 
495  // index map vector (basic serial code uses the map sln[i] = x[i] for all dofs.
496  int *idx = malloc_with_check<int>(m->get_size());
497  for (unsigned int i = 0; i < m->get_size(); i++) idx[i] = i;
498 
499  // copy solution to the output solution vector
500  vec_get_value(x, m->get_size(), idx, this->sln);
501  free_with_check(idx);
502 
503  KSPDestroy(&ksp);
504  VecDestroy(&x);
505  }
506 
507  template class HERMES_API PetscLinearMatrixSolver < double > ;
508  template class HERMES_API PetscLinearMatrixSolver < std::complex<double> > ;
509  }
510 }
511 #endif
General namespace for the Hermes library.
Exception interface Basically a std::exception, but with a constructor with string and with print_msg...
Definition: exceptions.h:49
File containing common definitions, and basic global enums etc. for HermesCommon. ...
MatrixExportFormat
Format of file matrix and vector output.
virtual void add_sparse_matrix(SparseMatrix< Scalar > *mat)
Definition: matrix.cpp:190
File containing functionality for investigating call stack.
void change_sign(T **matrix, unsigned int m, unsigned int n)
Changes the sign of a matrix.
PETSc solver interface.
File containing common definitions, and basic global enums etc. for HermesCommon. ...