gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
linearSystemPETSc.hpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #if defined(HAVE_PETSC)
7 
8 #include <petsc.h>
9 #include <petscksp.h>
10 #include "linearSystemPETSc.h"
11 
12 #if((PETSC_VERSION_RELEASE == 0) || \
13  ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 7)))
14 #define PetscOptionsGetInt(A, B, C, D) PetscOptionsGetInt(NULL, A, B, C, D)
15 #define PetscOptionsInsertString(A) PetscOptionsInsertString(NULL, A)
16 #endif
17 
18 static void _check(int ierr) { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
19 
20 template <class scalar>
21 int linearSystemPETSc<scalar>::_getBlockSizeFromParameters() const
22 {
23  return 1;
24 }
25 
26 template <class scalar> void linearSystemPETSc<scalar>::_kspCreate()
27 {
28  // Set option given by the user in its (python script) without using argc,argv
29  // or .petscrc
30  if(this->_parameters.count("petsc_solver_options"))
31  _check(PetscOptionsInsertString(
32  this->_parameters["petsc_solver_options"].c_str()));
33  _check(KSPCreate(_comm, &_ksp));
34  PC pc;
35  _check(KSPGetPC(_ksp, &pc));
36  // set some default options
37  //_check(PCSetType(pc, PCLU));//LU for direct solver and PCILU for indirect
38  //solver
39  /* _check(PCFactorSetMatOrderingType(pc, MATORDERING_RCM));
40  _check(PCFactorSetLevels(pc, 1));*/
41  _check(
42  KSPSetTolerances(_ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
43  // override the default options with the ones from the option
44  // database (if any)
45  if(this->_parameters.count("petscPrefix"))
46  _check(
47  KSPAppendOptionsPrefix(_ksp, this->_parameters["petscPrefix"].c_str()));
48  _check(KSPSetFromOptions(_ksp));
49  _check(PCSetFromOptions(pc));
50  _kspAllocated = true;
51 }
52 
53 template <class scalar>
54 int linearSystemPETSc<scalar>::getNumKspIteration() const
55 {
56  int n;
57  _check(KSPGetIterationNumber(_ksp, &n));
58  return n;
59 }
60 
61 template <class scalar>
62 linearSystemPETSc<scalar>::linearSystemPETSc(MPI_Comm com)
63 {
64  _comm = com;
65  _isAllocated = false;
66  _kspAllocated = false;
67  _matrixChangedSinceLastSolve = true;
68  _valuesNotAssembled = false;
69  _entriesPreAllocated = false;
70 }
71 
72 template <class scalar> linearSystemPETSc<scalar>::linearSystemPETSc()
73 {
74  _comm = PETSC_COMM_WORLD;
75  _isAllocated = false;
76  _kspAllocated = false;
77  _matrixChangedSinceLastSolve = true;
78  _valuesNotAssembled = false;
79  _entriesPreAllocated = false;
80 }
81 
82 template <class scalar> linearSystemPETSc<scalar>::~linearSystemPETSc()
83 {
84  clear();
85  if(_kspAllocated) _check(KSPDestroy(&_ksp));
86 }
87 
88 template <class scalar>
89 void linearSystemPETSc<scalar>::_assembleMatrixIfNeeded()
90 {
91 #if defined(HAVE_MPI)
92  if(_comm == PETSC_COMM_WORLD) {
93  if(Msg::GetCommSize() > 1) {
94  int value = _valuesNotAssembled ? 1 : 0;
95  int sumValue = 0;
96  MPI_Allreduce((void *)&value, (void *)&sumValue, 1, MPI_INT, MPI_SUM,
97  _comm);
98  if((sumValue > 0) && (sumValue < Msg::GetCommSize())) {
99  _valuesNotAssembled = 1;
100  }
101  }
102  }
103 #endif
104  if(_valuesNotAssembled) {
105  _check(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
106  _check(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
107  _matrixChangedSinceLastSolve = true;
108  _valuesNotAssembled = false;
109  }
110 }
111 
112 template <class scalar>
113 void linearSystemPETSc<scalar>::insertInSparsityPattern(int i, int j)
114 {
115  i -= _localRowStart;
116  if(i < 0 || i >= _localSize) return;
117  _sparsity.insertEntry(i, j);
118 }
119 
120 template <class scalar> void linearSystemPETSc<scalar>::preAllocateEntries()
121 {
122  if(_entriesPreAllocated) return;
123  if(!_isAllocated){
124  Msg::Error("System must be allocated first");
125  return;
126  }
127  int blockSize = _getBlockSizeFromParameters();
128  std::vector<int> nByRowDiag(_localSize), nByRowOffDiag(_localSize);
129  if(_sparsity.getNbRows() == 0) {
130  PetscInt prealloc = 300; // 8*27 = 216 for 8 2nd order hexas
131  PetscBool set;
132  PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
133  prealloc = std::min(prealloc, _localSize);
134  nByRowDiag.resize(0);
135  nByRowDiag.resize(_localSize, prealloc);
136  }
137  else {
138  for(int i = 0; i < _localSize; i++) {
139  int n;
140  const int *r = _sparsity.getRow(i, n);
141  for(int j = 0; j < n; j++) {
142  if(r[j] >= _localRowStart && r[j] < _localRowEnd)
143  nByRowDiag[i]++;
144  else
145  nByRowOffDiag[i]++;
146  }
147  }
148  _sparsity.clear();
149  }
150  // MatXAIJSetPreallocation is not available in petsc < 3.3
151  int commSize = 1;
152  MPI_Comm_size(_comm, &commSize);
153  if(commSize == 1) {
154  if(blockSize == 1)
155  _check(MatSeqAIJSetPreallocation(_a, 0, &nByRowDiag[0]));
156  else
157  _check(MatSeqBAIJSetPreallocation(_a, blockSize, 0, &nByRowDiag[0]));
158  }
159  else {
160  if(blockSize == 1)
161  _check(
162  MatMPIAIJSetPreallocation(_a, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]));
163  else
164  _check(MatMPIBAIJSetPreallocation(_a, blockSize, 0, &nByRowDiag[0], 0,
165  &nByRowOffDiag[0]));
166  }
167  if(blockSize > 1) _check(MatSetOption(_a, MAT_ROW_ORIENTED, PETSC_FALSE));
168  _entriesPreAllocated = true;
169 
170 #if ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 3))
171  // Preallocation routines automatically set now
172  // MAT_NEW_NONZERO_ALLOCATION_ERR, which causes a problem when the mask of the
173  // matrix changes. We must disable the error generation and allow new
174  // allocation (if needed)
175  _check(MatSetOption(_a, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
176 #endif
177 }
178 
179 template <class scalar> void linearSystemPETSc<scalar>::allocate(int nbRows)
180 {
181  int commSize;
182  MPI_Comm_size(_comm, &commSize);
183  int blockSize = _getBlockSizeFromParameters();
184  clear();
185  _check(MatCreate(_comm, &_a));
186  _check(MatSetSizes(_a, blockSize * nbRows, blockSize * nbRows,
187  PETSC_DETERMINE, PETSC_DETERMINE));
188  if(blockSize > 1) {
189  if(commSize > 1) {
190  _check(MatSetType(_a, MATMPIBAIJ));
191  }
192  else {
193  _check(MatSetType(_a, MATSEQBAIJ));
194  }
195  }
196  // override the default options with the ones from the option
197  // database (if any)
198  if(this->_parameters.count("petscOptions"))
199  _check(PetscOptionsInsertString(this->_parameters["petscOptions"].c_str()));
200  if(this->_parameters.count("petscPrefix"))
201  _check(
202  MatAppendOptionsPrefix(_a, this->_parameters["petscPrefix"].c_str()));
203  _check(MatSetFromOptions(_a));
204  // since PETSc 3.3 GetOwnershipRange and MatGetSize cannot be called before
205  // MatXXXSetPreallocation
206  _localSize = nbRows;
207 #if defined(HAVE_MPI)
208  if(commSize > 1) {
209  _localRowStart = 0;
210  if(Msg::GetCommRank() != 0) {
211  MPI_Status status;
212  MPI_Recv((void *)&_localRowStart, 1, MPI_INT, Msg::GetCommRank() - 1, 1,
213  MPI_COMM_WORLD, &status);
214  }
215  _localRowEnd = _localRowStart + nbRows;
216  if(Msg::GetCommRank() != Msg::GetCommSize() - 1) {
217  MPI_Send((void *)&_localRowEnd, 1, MPI_INT, Msg::GetCommRank() + 1, 1,
218  MPI_COMM_WORLD);
219  }
220  MPI_Allreduce((void *)&_localSize, (void *)&_globalSize, 1, MPI_INT,
221  MPI_SUM, MPI_COMM_WORLD);
222  }
223  else {
224  _localRowStart = 0;
225  _localRowEnd = nbRows;
226  _globalSize = _localSize;
227  }
228 #else
229  _localRowStart = 0;
230  _localRowEnd = nbRows;
231  _globalSize = _localSize;
232 #endif
233 
234  // preallocation option must be set after other options
235  _check(VecCreate(_comm, &_x));
236  _check(VecSetSizes(_x, blockSize * nbRows, PETSC_DETERMINE));
237  // override the default options with the ones from the option
238  // database (if any)
239  if(this->_parameters.count("petscPrefix"))
240  _check(
241  VecAppendOptionsPrefix(_x, this->_parameters["petscPrefix"].c_str()));
242  _check(VecSetFromOptions(_x));
243  _check(VecDuplicate(_x, &_b));
244  _isAllocated = true;
245  _entriesPreAllocated = false;
246 }
247 
248 template <class scalar> void linearSystemPETSc<scalar>::print()
249 {
250  _assembleMatrixIfNeeded();
251  _check(VecAssemblyBegin(_b));
252  _check(VecAssemblyEnd(_b));
253 
254  /*
255  PetscViewer fd;
256  _check(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "mat.m", &fd));
257  _check(PetscViewerSetFormat(fd, PETSC_VIEWER_ASCII_MATLAB));
258  _check(PetscObjectSetName((PetscObject)_a, "A"));
259  _check(MatView(_a, fd));
260  */
261 
262  if(Msg::GetCommRank() == 0) printf("a :\n");
263  MatView(_a, PETSC_VIEWER_STDOUT_WORLD);
264  if(Msg::GetCommRank() == 0) printf("b :\n");
265  VecView(_b, PETSC_VIEWER_STDOUT_WORLD);
266  if(Msg::GetCommRank() == 0) printf("x :\n");
267  VecView(_x, PETSC_VIEWER_STDOUT_WORLD);
268 }
269 
270 template <class scalar> void linearSystemPETSc<scalar>::clear()
271 {
272  if(_isAllocated) {
273  _check(MatDestroy(&_a));
274  _check(VecDestroy(&_x));
275  _check(VecDestroy(&_b));
276  }
277  _isAllocated = false;
278 }
279 
280 template <class scalar>
281 void linearSystemPETSc<scalar>::getFromMatrix(int row, int col,
282  scalar &val) const
283 {
284  Msg::Error("getFromMatrix not implemented for PETSc");
285 }
286 
287 template <class scalar>
288 void linearSystemPETSc<scalar>::addToRightHandSide(int row, const scalar &val,
289  int ith)
290 {
291  PetscInt i = row;
292  PetscScalar s = val;
293  _check(VecSetValues(_b, 1, &i, &s, ADD_VALUES));
294 }
295 #if defined(PETSC_USE_COMPLEX)
296 // this specialization will cast to a double (take the real part) if "val" is a
297 // double wheras Petsc is built in complex mode.
298 template <>
299 void linearSystemPETSc<double>::getFromRightHandSide(int row,
300  double &val) const;
301 #endif
302 
303 template <class scalar>
304 void linearSystemPETSc<scalar>::getFromRightHandSide(int row, scalar &val) const
305 {
306  _check(VecGetValues(_b, 1, &row, &val));
307 }
308 
309 template <class scalar>
310 double linearSystemPETSc<scalar>::normInfRightHandSide() const
311 {
312  PetscReal nor;
313  VecAssemblyBegin(_b);
314  VecAssemblyEnd(_b);
315  _check(VecNorm(_b, NORM_INFINITY, &nor));
316  return nor;
317 }
318 
319 template <class scalar>
320 void linearSystemPETSc<scalar>::addToMatrix(int row, int col, const scalar &val)
321 {
322  if(!_entriesPreAllocated) preAllocateEntries();
323  PetscInt i = row, j = col;
324  PetscScalar s = val;
325  _check(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES));
326  _valuesNotAssembled = true;
327 }
328 
329 template <class scalar>
330 void linearSystemPETSc<scalar>::addToSolution(int row, const scalar &val)
331 {
332  PetscInt i = row;
333  PetscScalar s = val;
334  _check(VecSetValues(_x, 1, &i, &s, ADD_VALUES));
335 }
336 #if defined(PETSC_USE_COMPLEX)
337 // this specialization will cast to a double (take the real part) if "val" is a
338 // double wheras Petsc is built in complex mode.
339 template <>
340 void linearSystemPETSc<double>::getFromSolution(int row, double &val) const;
341 #endif
342 
343 template <class scalar>
344 void linearSystemPETSc<scalar>::getFromSolution(int row, scalar &val) const
345 {
346  _check(VecGetValues(_x, 1, &row, &val));
347 }
348 
349 template <class scalar> void linearSystemPETSc<scalar>::zeroMatrix()
350 {
351 #if defined(HAVE_MPI)
352  if(_comm == PETSC_COMM_WORLD) {
353  if(Msg::GetCommSize() > 1) {
354  int value = _entriesPreAllocated ? 1 : 0;
355  int sumValue = 0;
356  MPI_Allreduce((void *)&value, (void *)&sumValue, 1, MPI_INT, MPI_SUM,
357  _comm);
358  if((sumValue >= 0) && (sumValue < Msg::GetCommSize()) &&
359  !_entriesPreAllocated) {
360  preAllocateEntries();
361  }
362  }
363  }
364 #endif
365  if(_isAllocated && _entriesPreAllocated) {
366  _assembleMatrixIfNeeded();
367  _check(MatZeroEntries(_a));
368  }
369 }
370 
371 template <class scalar> void linearSystemPETSc<scalar>::zeroRightHandSide()
372 {
373  if(_isAllocated) {
374  _check(VecAssemblyBegin(_b));
375  _check(VecAssemblyEnd(_b));
376  _check(VecZeroEntries(_b));
377  }
378 }
379 
380 template <class scalar> void linearSystemPETSc<scalar>::zeroSolution()
381 {
382  if(_isAllocated) {
383  _check(VecAssemblyBegin(_x));
384  _check(VecAssemblyEnd(_x));
385  _check(VecZeroEntries(_x));
386  }
387 }
388 
389 template <class scalar> int linearSystemPETSc<scalar>::systemSolve()
390 {
391  if(!_kspAllocated) _kspCreate();
392  _assembleMatrixIfNeeded();
393  if(!_matrixChangedSinceLastSolve ||
394  linearSystem<scalar>::_parameters["matrix_reuse"] == "same_matrix")
395  _check(KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER));
396  else if(linearSystem<scalar>::_parameters["matrix_reuse"] == "same_sparsity")
397  _check(KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN));
398  else
399  _check(KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN));
400  _matrixChangedSinceLastSolve = false;
401  /*MatInfo info;
402  MatGetInfo(_a, MAT_LOCAL, &info);
403  printf("mallocs %.0f nz_allocated %.0f nz_used %.0f nz_unneeded
404  %.0f\n", info.mallocs, info.nz_allocated, info.nz_used, info.nz_unneeded);*/
405  _check(VecAssemblyBegin(_b));
406  _check(VecAssemblyEnd(_b));
407  _check(KSPSolve(_ksp, _b, _x));
408  //_check(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF));
409  // PetscInt its;
410  //_check(KSPGetIterationNumber(ksp, &its));
411  // Msg::Info("%d iterations", its);
412  return 1;
413 }
414 
415 template <class scalar> int linearSystemPETSc<scalar>::matMult()
416 {
417  _check(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
418  _check(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
419  _check(VecAssemblyBegin(_b));
420  _check(VecAssemblyEnd(_b));
421  _check(MatMult(_a, _b, _x));
422  return 1;
423 }
424 
425 template <class scalar>
426 void linearSystemPETSc<scalar>::printMatlab(const char *filename) const
427 {
428  _check(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
429  _check(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
430  _check(VecAssemblyBegin(_b));
431  _check(VecAssemblyEnd(_b));
432 
433  PetscViewer viewer;
434  PetscViewerASCIIOpen(PETSC_COMM_WORLD, filename, &viewer);
435  PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
436  MatView(_a, viewer);
437  PetscViewerDestroy(&viewer);
438  return;
439 }
440 
441 /*template <class scalar>
442 std::vector<scalar> linearSystemPETSc<scalar>::getData()
443 {
444  _assembleMatrixIfNeeded();
445  PetscScalar *v;
446  _check(MatGetArray(_a,&v));
447  MatInfo info;
448  _check(MatGetInfo(_a,MAT_LOCAL,&info));
449  std::vector<scalar> data; // Maybe I should reserve or resize (SAM)
450 
451 #if defined(PETSC_USE_COMPLEX)
452  for (int i = 0; i < info.nz_allocated; i++)
453  data.push_back(v[i].real());
454 #else
455  for (int i = 0; i < info.nz_allocated; i++)
456  data.push_back(v[i]);
457 #endif
458  _check(MatRestoreArray(_a,&v));
459  return data;
460 }*/
461 /*
462 template <class scalar>
463 std::vector<int> linearSystemPETSc<scalar>::getRowPointers()
464 {
465  _assembleMatrixIfNeeded();
466  const PetscInt *rows;
467  const PetscInt *columns;
468  PetscInt n;
469  PetscBool done;
470  _check(MatGetRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done));
471 //case done == PETSC_FALSE should be handled std::vector<int> rowPointers; //
472 Maybe I should reserve or resize (SAM) for (int i = 0; i <= n; i++)
473  rowPointers.push_back(rows[i]);
474  _check(MatRestoreRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done));
475  return rowPointers;
476 }
477 
478 template <class scalar>
479 std::vector<int> linearSystemPETSc<scalar>::getColumnsIndices()
480 {
481  _assembleMatrixIfNeeded();
482  const PetscInt *rows;
483  const PetscInt *columns;
484  PetscInt n;
485  PetscBool done;
486  _check(MatGetRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done));
487 //case done == PETSC_FALSE should be handled MatInfo info;
488  _check(MatGetInfo(_a,MAT_LOCAL,&info));
489  std::vector<int> columnIndices; // Maybe I should reserve or resize (SAM)
490  for (int i = 0; i < info.nz_allocated; i++)
491  columnIndices.push_back(columns[i]);
492  _check(MatRestoreRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done));
493  return columnIndices;
494 }
495 */
496 #endif
linearSystemPETSc::getNumKspIteration
int getNumKspIteration() const
Definition: linearSystemPETSc.h:165
linearSystemPETSc::linearSystemPETSc
linearSystemPETSc()
Definition: linearSystemPETSc.h:152
linearSystemPETSc
Definition: linearSystemPETSc.h:150
linearSystemPETSc::zeroRightHandSide
void zeroRightHandSide()
Definition: linearSystemPETSc.h:167
linearSystemPETSc::getFromMatrix
void getFromMatrix(int row, int col, scalar &val) const
Definition: linearSystemPETSc.h:160
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
linearSystemPETSc.h
linearSystemPETSc::addToMatrix
void addToMatrix(int row, int col, const scalar &val)
Definition: linearSystemPETSc.h:159
linearSystemPETSc::addToRightHandSide
void addToRightHandSide(int row, const scalar &val, int ith=0)
Definition: linearSystemPETSc.h:161
linearSystemPETSc::printMatlab
void printMatlab(const char *filename) const
Definition: linearSystemPETSc.h:169
linearSystemPETSc::getFromRightHandSide
void getFromRightHandSide(int row, scalar &val) const
Definition: linearSystemPETSc.h:163
linearSystemPETSc::allocate
void allocate(int nbRows)
Definition: linearSystemPETSc.h:157
linearSystemBase::preAllocateEntries
virtual void preAllocateEntries()
Definition: linearSystem.h:22
Msg::GetCommRank
static int GetCommRank()
Definition: GmshMessage.cpp:219
linearSystemPETSc::normInfRightHandSide
double normInfRightHandSide() const
Definition: linearSystemPETSc.h:171
linearSystemPETSc::zeroSolution
void zeroSolution()
Definition: linearSystemPETSc.h:168
linearSystemPETSc::clear
void clear()
Definition: linearSystemPETSc.h:158
linearSystemBase::insertInSparsityPattern
virtual void insertInSparsityPattern(int _row, int _col)
Definition: linearSystem.h:33
Msg::GetCommSize
static int GetCommSize()
Definition: GmshMessage.cpp:224
linearSystemPETSc::matMult
virtual int matMult()
Definition: linearSystemPETSc.h:172
linearSystemPETSc::systemSolve
virtual int systemSolve()
Definition: linearSystemPETSc.h:170
linearSystemPETSc::addToSolution
void addToSolution(int row, const scalar &val)
Definition: linearSystemPETSc.h:162
linearSystemPETSc::zeroMatrix
void zeroMatrix()
Definition: linearSystemPETSc.h:166
linearSystemPETSc::getFromSolution
void getFromSolution(int row, scalar &val) const
Definition: linearSystemPETSc.h:164
linearSystem
Definition: linearSystem.h:38