DOLFINx
DOLFINx C++ interface
petsc.h
1// Copyright (C) 2004-2018 Johan Hoffman, Johan Jansson, Anders Logg and
2// Garth N. Wells
3//
4// This file is part of DOLFINx (https://www.fenicsproject.org)
5//
6// SPDX-License-Identifier: LGPL-3.0-or-later
7
8#pragma once
9
10#include "Vector.h"
11#include "utils.h"
12#include <boost/lexical_cast.hpp>
13#include <functional>
14#include <petscksp.h>
15#include <petscmat.h>
16#include <petscoptions.h>
17#include <petscvec.h>
18#include <string>
19#include <vector>
20#include <xtl/xspan.hpp>
21
22namespace dolfinx::common
23{
24class IndexMap;
25}
26
27namespace dolfinx::la
28{
29class SparsityPattern;
30
31namespace petsc
32{
34void error(int error_code, std::string filename, std::string petsc_function);
35
43std::vector<Vec>
44create_vectors(MPI_Comm comm,
45 const std::vector<xtl::span<const PetscScalar>>& x);
46
52Vec create_vector(const common::IndexMap& map, int bs);
53
62Vec create_vector(MPI_Comm comm, std::array<std::int64_t, 2> range,
63 const xtl::span<const std::int64_t>& ghosts, int bs);
64
74Vec create_vector_wrap(const common::IndexMap& map, int bs,
75 const xtl::span<const PetscScalar>& x);
76
80template <typename Allocator>
81Vec create_vector_wrap(const la::Vector<PetscScalar, Allocator>& x)
82{
83 assert(x.map());
84 return create_vector_wrap(*x.map(), x.bs(), x.array());
85}
86
98std::vector<IS> create_index_sets(
99 const std::vector<
100 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
101
103std::vector<std::vector<PetscScalar>> get_local_vectors(
104 const Vec x,
105 const std::vector<
106 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
107
109void scatter_local_vectors(
110 Vec x, const std::vector<xtl::span<const PetscScalar>>& x_b,
111 const std::vector<
112 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
115Mat create_matrix(MPI_Comm comm, const SparsityPattern& sp,
116 const std::string& type = std::string());
117
123MatNullSpace create_nullspace(MPI_Comm comm, const xtl::span<const Vec>& basis);
124
130namespace options
131{
133void set(std::string option);
134
136template <typename T>
137void set(std::string option, const T value)
138{
139 if (option[0] != '-')
140 option = '-' + option;
141
142 PetscErrorCode ierr;
143 ierr = PetscOptionsSetValue(nullptr, option.c_str(),
144 boost::lexical_cast<std::string>(value).c_str());
145 if (ierr != 0)
146 petsc::error(ierr, __FILE__, "PetscOptionsSetValue");
147}
148
150void clear(std::string option);
151
153void clear();
154} // namespace options
155
162{
163public:
168 Vector(const common::IndexMap& map, int bs);
169
170 // Delete copy constructor to avoid accidental copying of 'heavy' data
171 Vector(const Vector& x) = delete;
172
174 Vector(Vector&& x);
175
187 Vector(Vec x, bool inc_ref_count);
188
190 virtual ~Vector();
191
192 // Assignment operator (disabled)
193 Vector& operator=(const Vector& x) = delete;
194
196 Vector& operator=(Vector&& x);
197
200 Vector copy() const;
201
203 std::int64_t size() const;
204
206 std::int32_t local_size() const;
207
209 std::array<std::int64_t, 2> local_range() const;
210
212 MPI_Comm comm() const;
213
215 void set_options_prefix(std::string options_prefix);
216
219 std::string get_options_prefix() const;
220
222 void set_from_options();
223
225 Vec vec() const;
226
227private:
228 // PETSc Vec pointer
229 Vec _x;
230};
231
235{
236public:
238 Operator(Mat A, bool inc_ref_count);
239
240 // Copy constructor (deleted)
241 Operator(const Operator& A) = delete;
242
244 Operator(Operator&& A);
245
247 virtual ~Operator();
248
250 Operator& operator=(const Operator& A) = delete;
251
254
257 std::array<std::int64_t, 2> size() const;
258
264 Vec create_vector(std::size_t dim) const;
265
267 Mat mat() const;
268
269protected:
270 // PETSc Mat pointer
271 Mat _matA;
272};
273
279class Matrix : public Operator
280{
281public:
286 static auto set_fn(Mat A, InsertMode mode)
287 {
288 return [A, mode, cache = std::vector<PetscInt>()](
289 const xtl::span<const std::int32_t>& rows,
290 const xtl::span<const std::int32_t>& cols,
291 const xtl::span<const PetscScalar>& vals) mutable -> int
292 {
293 PetscErrorCode ierr;
294#ifdef PETSC_USE_64BIT_INDICES
295 cache.resize(rows.size() + cols.size());
296 std::copy(rows.begin(), rows.end(), cache.begin());
297 std::copy(cols.begin(), cols.end(),
298 std::next(cache.begin(), rows.size()));
299 const PetscInt* _rows = cache.data();
300 const PetscInt* _cols = cache.data() + rows.size();
301 ierr = MatSetValuesLocal(A, rows.size(), _rows, cols.size(), _cols,
302 vals.data(), mode);
303#else
304 ierr = MatSetValuesLocal(A, rows.size(), rows.data(), cols.size(),
305 cols.data(), vals.data(), mode);
306#endif
307
308#ifndef NDEBUG
309 if (ierr != 0)
310 petsc::error(ierr, __FILE__, "MatSetValuesLocal");
311#endif
312 return ierr;
313 };
314 }
315
321 static auto set_block_fn(Mat A, InsertMode mode)
322 {
323 return [A, mode, cache = std::vector<PetscInt>()](
324 const xtl::span<const std::int32_t>& rows,
325 const xtl::span<const std::int32_t>& cols,
326 const xtl::span<const PetscScalar>& vals) mutable -> int
327 {
328 PetscErrorCode ierr;
329#ifdef PETSC_USE_64BIT_INDICES
330 cache.resize(rows.size() + cols.size());
331 std::copy(rows.begin(), rows.end(), cache.begin());
332 std::copy(cols.begin(), cols.end(),
333 std::next(cache.begin(), rows.size()));
334 const PetscInt* _rows = cache.data();
335 const PetscInt* _cols = cache.data() + rows.size();
336 ierr = MatSetValuesBlockedLocal(A, rows.size(), _rows, cols.size(), _cols,
337 vals.data(), mode);
338#else
339 ierr = MatSetValuesBlockedLocal(A, rows.size(), rows.data(), cols.size(),
340 cols.data(), vals.data(), mode);
341#endif
342
343#ifndef NDEBUG
344 if (ierr != 0)
345 petsc::error(ierr, __FILE__, "MatSetValuesBlockedLocal");
346#endif
347 return ierr;
348 };
349 }
350
359 static auto set_block_expand_fn(Mat A, int bs0, int bs1, InsertMode mode)
360 {
361 return [A, bs0, bs1, mode, cache0 = std::vector<PetscInt>(),
362 cache1 = std::vector<PetscInt>()](
363 const xtl::span<const std::int32_t>& rows,
364 const xtl::span<const std::int32_t>& cols,
365 const xtl::span<const PetscScalar>& vals) mutable -> int
366 {
367 PetscErrorCode ierr;
368 cache0.resize(bs0 * rows.size());
369 cache1.resize(bs1 * cols.size());
370 for (std::size_t i = 0; i < rows.size(); ++i)
371 for (int k = 0; k < bs0; ++k)
372 cache0[bs0 * i + k] = bs0 * rows[i] + k;
373
374 for (std::size_t i = 0; i < cols.size(); ++i)
375 for (int k = 0; k < bs1; ++k)
376 cache1[bs1 * i + k] = bs1 * cols[i] + k;
377
378 ierr = MatSetValuesLocal(A, cache0.size(), cache0.data(), cache1.size(),
379 cache1.data(), vals.data(), mode);
380#ifndef NDEBUG
381 if (ierr != 0)
382 petsc::error(ierr, __FILE__, "MatSetValuesLocal");
383#endif
384 return ierr;
385 };
386 }
387
389 Matrix(MPI_Comm comm, const SparsityPattern& sp,
390 const std::string& type = std::string());
391
396 Matrix(Mat A, bool inc_ref_count);
397
398 // Copy constructor (deleted)
399 Matrix(const Matrix& A) = delete;
400
402 Matrix(Matrix&& A) = default;
403
405 ~Matrix() = default;
406
408 Matrix& operator=(const Matrix& A) = delete;
409
411 Matrix& operator=(Matrix&& A) = default;
412
416 enum class AssemblyType : std::int32_t
417 {
418 FINAL,
419 FLUSH
420 };
421
427 void apply(AssemblyType type);
428
430 double norm(Norm norm_type) const;
431
432 //--- Special PETSc Functions ---
433
436 void set_options_prefix(std::string options_prefix);
437
440 std::string get_options_prefix() const;
441
443 void set_from_options();
444};
445
449{
450public:
453 explicit KrylovSolver(MPI_Comm comm);
454
458 KrylovSolver(KSP ksp, bool inc_ref_count);
459
460 // Copy constructor (deleted)
461 KrylovSolver(const KrylovSolver& solver) = delete;
462
464 KrylovSolver(KrylovSolver&& solver);
465
468
469 // Assignment operator (deleted)
470 KrylovSolver& operator=(const KrylovSolver&) = delete;
471
473 KrylovSolver& operator=(KrylovSolver&& solver);
474
476 void set_operator(const Mat A);
477
479 void set_operators(const Mat A, const Mat P);
480
483 int solve(Vec x, const Vec b, bool transpose = false) const;
484
487 void set_options_prefix(std::string options_prefix);
488
491 std::string get_options_prefix() const;
492
494 void set_from_options() const;
495
497 KSP ksp() const;
498
499private:
500 // PETSc solver pointer
501 KSP _ksp;
502};
503} // namespace petsc
504} // namespace dolfinx::la
This class represents the distribution index arrays across processes. An index array is a contiguous ...
Definition: IndexMap.h:60
This class provides a sparsity pattern data structure that can be used to initialize sparse matrices....
Definition: SparsityPattern.h:34
This class implements Krylov methods for linear systems of the form Ax = b. It is a wrapper for the K...
Definition: petsc.h:449
void set_operator(const Mat A)
Set operator (Mat)
Definition: petsc.cpp:652
KSP ksp() const
Return PETSc KSP pointer.
Definition: petsc.cpp:780
std::string get_options_prefix() const
Returns the prefix used by PETSc when searching the PETSc options database.
Definition: petsc.cpp:762
void set_operators(const Mat A, const Mat P)
Set operator and preconditioner matrix (Mat)
Definition: petsc.cpp:654
KrylovSolver(MPI_Comm comm)
Create Krylov solver for a particular method and named preconditioner.
Definition: petsc.cpp:615
void set_from_options() const
Set options from PETSc options database.
Definition: petsc.cpp:772
~KrylovSolver()
Destructor.
Definition: petsc.cpp:640
int solve(Vec x, const Vec b, bool transpose=false) const
Solve linear system Ax = b and return number of iterations (A^t x = b if transpose is true)
Definition: petsc.cpp:663
void set_options_prefix(std::string options_prefix)
Sets the prefix used by PETSc when searching the PETSc options database.
Definition: petsc.cpp:753
It is a simple wrapper for a PETSc matrix pointer (Mat). Its main purpose is to assist memory managem...
Definition: petsc.h:280
static auto set_block_fn(Mat A, InsertMode mode)
Return a function with an interface for adding or inserting values into the matrix A using blocked in...
Definition: petsc.h:321
void set_from_options()
Call PETSc function MatSetFromOptions on the PETSc Mat object.
Definition: petsc.cpp:608
double norm(Norm norm_type) const
Return norm of matrix.
Definition: petsc.cpp:552
static auto set_fn(Mat A, InsertMode mode)
Return a function with an interface for adding or inserting values into the matrix A (calls MatSetVal...
Definition: petsc.h:286
Matrix(MPI_Comm comm, const SparsityPattern &sp, const std::string &type=std::string())
Create holder for a PETSc Mat object from a sparsity pattern.
Definition: petsc.cpp:540
Matrix(Matrix &&A)=default
Move constructor (falls through to base class move constructor)
std::string get_options_prefix() const
Returns the prefix used by PETSc when searching the options database.
Definition: petsc.cpp:600
Matrix & operator=(Matrix &&A)=default
Move assignment operator.
static auto set_block_expand_fn(Mat A, int bs0, int bs1, InsertMode mode)
Return a function with an interface for adding or inserting blocked values to the matrix A using non-...
Definition: petsc.h:359
~Matrix()=default
Destructor.
Matrix & operator=(const Matrix &A)=delete
Assignment operator (deleted)
AssemblyType
Assembly type FINAL - corresponds to PETSc MAT_FINAL_ASSEMBLY FLUSH - corresponds to PETSc MAT_FLUSH_...
Definition: petsc.h:417
void apply(AssemblyType type)
Finalize assembly of tensor. The following values are recognized for the mode parameter:
Definition: petsc.cpp:577
void set_options_prefix(std::string options_prefix)
Sets the prefix used by PETSc when searching the options database.
Definition: petsc.cpp:594
This class is a base class for matrices that can be used in petsc::KrylovSolver.
Definition: petsc.h:235
std::array< std::int64_t, 2 > size() const
Return number of rows and columns (num_rows, num_cols). PETSc returns -1 if size has not been set.
Definition: petsc.cpp:499
Operator & operator=(const Operator &A)=delete
Assignment operator (deleted)
Operator(Mat A, bool inc_ref_count)
Constructor.
Definition: petsc.cpp:474
Vec create_vector(std::size_t dim) const
Initialize vector to be compatible with the matrix-vector product y = Ax. In the parallel case,...
Definition: petsc.cpp:509
Mat mat() const
Return PETSc Mat pointer.
Definition: petsc.cpp:537
virtual ~Operator()
Destructor.
Definition: petsc.cpp:485
A simple wrapper for a PETSc vector pointer (Vec). Its main purpose is to assist with memory/lifetime...
Definition: petsc.h:162
void set_from_options()
Call PETSc function VecSetFromOptions on the underlying Vec object.
Definition: petsc.cpp:464
std::int64_t size() const
Return global size of the vector.
Definition: petsc.cpp:411
std::string get_options_prefix() const
Returns the prefix used by PETSc when searching the options database.
Definition: petsc.cpp:455
Vector copy() const
Create a copy of the vector.
Definition: petsc.cpp:401
virtual ~Vector()
Destructor.
Definition: petsc.cpp:389
Vector(const common::IndexMap &map, int bs)
Create a vector.
Definition: petsc.cpp:374
MPI_Comm comm() const
Return MPI communicator.
Definition: petsc.cpp:439
std::array< std::int64_t, 2 > local_range() const
Return ownership range for calling rank.
Definition: petsc.cpp:429
std::int32_t local_size() const
Return local size of vector (belonging to the call rank)
Definition: petsc.cpp:420
Vec vec() const
Return pointer to PETSc Vec object.
Definition: petsc.cpp:471
void set_options_prefix(std::string options_prefix)
Sets the prefix used by PETSc when searching the options database.
Definition: petsc.cpp:448
Miscellaneous classes, functions and types.
Mat create_matrix(const Form< PetscScalar > &a, const std::string &type=std::string())
Create a matrix.
Definition: petsc.cpp:19
void clear(std::string option)
Clear a PETSc option.
Definition: petsc.cpp:354
void set(std::string option)
Set PETSc option that takes no value.
Definition: petsc.cpp:349
Linear algebra interface.
Definition: sparsitybuild.h:13
Norm
Norm types.
Definition: utils.h:13