DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
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#ifdef HAS_PETSC
11
12#include "Vector.h"
13#include "utils.h"
14#include <boost/lexical_cast.hpp>
15#include <functional>
16#include <petscksp.h>
17#include <petscmat.h>
18#include <petscoptions.h>
19#include <petscvec.h>
20#include <span>
21#include <string>
22#include <vector>
23
24namespace dolfinx::common
25{
26class IndexMap;
27} // namespace dolfinx::common
28
29namespace dolfinx::la
30{
31class SparsityPattern;
32
34namespace petsc
35{
37void error(int error_code, std::string filename, std::string petsc_function);
38
46std::vector<Vec>
47create_vectors(MPI_Comm comm,
48 const std::vector<std::span<const PetscScalar>>& x);
49
55Vec create_vector(const common::IndexMap& map, int bs);
56
65Vec create_vector(MPI_Comm comm, std::array<std::int64_t, 2> range,
66 std::span<const std::int64_t> ghosts, int bs);
67
77Vec create_vector_wrap(const common::IndexMap& map, int bs,
78 std::span<const PetscScalar> x);
79
83template <class V>
85{
86 assert(x.index_map());
87 return create_vector_wrap(*x.index_map(), x.bs(), x.array());
88}
89
102std::vector<IS> create_index_sets(
103 const std::vector<
104 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
105
107std::vector<std::vector<PetscScalar>> get_local_vectors(
108 const Vec x,
109 const std::vector<
110 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
111
114 Vec x, const std::vector<std::span<const PetscScalar>>& x_b,
115 const std::vector<
116 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
117
120Mat create_matrix(MPI_Comm comm, const SparsityPattern& sp,
121 std::string type = std::string());
122
128MatNullSpace create_nullspace(MPI_Comm comm, std::span<const Vec> basis);
129
135namespace options
136{
138void set(std::string option);
139
141template <typename T>
142void set(std::string option, const T value)
143{
144 if (option[0] != '-')
145 option = '-' + option;
146
147 PetscErrorCode ierr;
148 ierr = PetscOptionsSetValue(nullptr, option.c_str(),
149 boost::lexical_cast<std::string>(value).c_str());
150 if (ierr != 0)
151 petsc::error(ierr, __FILE__, "PetscOptionsSetValue");
152}
153
155void clear(std::string option);
156
158void clear();
159} // namespace options
160
167{
168public:
173 Vector(const common::IndexMap& map, int bs);
174
175 // Delete copy constructor to avoid accidental copying of 'heavy' data
176 Vector(const Vector& x) = delete;
177
179 Vector(Vector&& x);
180
192 Vector(Vec x, bool inc_ref_count);
193
195 virtual ~Vector();
196
197 // Assignment operator (disabled)
198 Vector& operator=(const Vector& x) = delete;
199
201 Vector& operator=(Vector&& x);
202
205 Vector copy() const;
206
208 std::int64_t size() const;
209
211 std::int32_t local_size() const;
212
214 std::array<std::int64_t, 2> local_range() const;
215
217 MPI_Comm comm() const;
218
220 void set_options_prefix(std::string options_prefix);
221
224 std::string get_options_prefix() const;
225
227 void set_from_options();
228
230 Vec vec() const;
231
232private:
233 // PETSc Vec pointer
234 Vec _x;
235};
236
240{
241public:
243 Operator(Mat A, bool inc_ref_count);
244
245 // Copy constructor (deleted)
246 Operator(const Operator& A) = delete;
247
249 Operator(Operator&& A);
250
252 virtual ~Operator();
253
255 Operator& operator=(const Operator& A) = delete;
256
259
262 std::array<std::int64_t, 2> size() const;
263
269 Vec create_vector(std::size_t dim) const;
270
272 Mat mat() const;
273
274protected:
275 // PETSc Mat pointer
276 Mat _matA;
277};
278
284class Matrix : public Operator
285{
286public:
291 static auto set_fn(Mat A, InsertMode mode)
292 {
293 return [A, mode, cache = std::vector<PetscInt>()](
294 std::span<const std::int32_t> rows,
295 std::span<const std::int32_t> cols,
296 std::span<const PetscScalar> vals) mutable -> int
297 {
298 PetscErrorCode ierr;
299#ifdef PETSC_USE_64BIT_INDICES
300 cache.resize(rows.size() + cols.size());
301 std::ranges::copy(rows, cache.begin());
302 std::ranges::copy(cols, std::next(cache.begin(), rows.size()));
303 const PetscInt* _rows = cache.data();
304 const PetscInt* _cols = cache.data() + rows.size();
305 ierr = MatSetValuesLocal(A, rows.size(), _rows, cols.size(), _cols,
306 vals.data(), mode);
307#else
308 ierr = MatSetValuesLocal(A, rows.size(), rows.data(), cols.size(),
309 cols.data(), vals.data(), mode);
310#endif
311
312#ifndef NDEBUG
313 if (ierr != 0)
314 petsc::error(ierr, __FILE__, "MatSetValuesLocal");
315#endif
316 return ierr;
317 };
318 }
319
325 static auto set_block_fn(Mat A, InsertMode mode)
326 {
327 return [A, mode, cache = std::vector<PetscInt>()](
328 std::span<const std::int32_t> rows,
329 std::span<const std::int32_t> cols,
330 std::span<const PetscScalar> vals) mutable -> int
331 {
332 PetscErrorCode ierr;
333#ifdef PETSC_USE_64BIT_INDICES
334 cache.resize(rows.size() + cols.size());
335 std::ranges::copy(rows, cache.begin());
336 std::ranges::copy(cols, std::next(cache.begin(), rows.size()));
337 const PetscInt* _rows = cache.data();
338 const PetscInt* _cols = cache.data() + rows.size();
339 ierr = MatSetValuesBlockedLocal(A, rows.size(), _rows, cols.size(), _cols,
340 vals.data(), mode);
341#else
342 ierr = MatSetValuesBlockedLocal(A, rows.size(), rows.data(), cols.size(),
343 cols.data(), vals.data(), mode);
344#endif
345
346#ifndef NDEBUG
347 if (ierr != 0)
348 petsc::error(ierr, __FILE__, "MatSetValuesBlockedLocal");
349#endif
350 return ierr;
351 };
352 }
353
362 static auto set_block_expand_fn(Mat A, int bs0, int bs1, InsertMode mode)
363 {
364 return [A, bs0, bs1, mode, cache0 = std::vector<PetscInt>(),
365 cache1 = std::vector<PetscInt>()](
366 std::span<const std::int32_t> rows,
367 std::span<const std::int32_t> cols,
368 std::span<const PetscScalar> vals) mutable -> int
369 {
370 PetscErrorCode ierr;
371 cache0.resize(bs0 * rows.size());
372 cache1.resize(bs1 * cols.size());
373 for (std::size_t i = 0; i < rows.size(); ++i)
374 for (int k = 0; k < bs0; ++k)
375 cache0[bs0 * i + k] = bs0 * rows[i] + k;
376
377 for (std::size_t i = 0; i < cols.size(); ++i)
378 for (int k = 0; k < bs1; ++k)
379 cache1[bs1 * i + k] = bs1 * cols[i] + k;
380
381 ierr = MatSetValuesLocal(A, cache0.size(), cache0.data(), cache1.size(),
382 cache1.data(), vals.data(), mode);
383#ifndef NDEBUG
384 if (ierr != 0)
385 petsc::error(ierr, __FILE__, "MatSetValuesLocal");
386#endif
387 return ierr;
388 };
389 }
390
392 Matrix(MPI_Comm comm, const SparsityPattern& sp,
393 std::string type = std::string());
394
399 Matrix(Mat A, bool inc_ref_count);
400
401 // Copy constructor (deleted)
402 Matrix(const Matrix& A) = delete;
403
405 Matrix(Matrix&& A) = default;
406
408 ~Matrix() = default;
409
411 Matrix& operator=(const Matrix& A) = delete;
412
414 Matrix& operator=(Matrix&& A) = default;
415
419 enum class AssemblyType : std::int32_t
420 {
421 FINAL,
422 FLUSH
423 };
424
430 void apply(AssemblyType type);
431
433 double norm(Norm norm_type) const;
434
435 //--- Special PETSc Functions ---
436
439 void set_options_prefix(std::string options_prefix);
440
443 std::string get_options_prefix() const;
444
446 void set_from_options();
447};
448
452{
453public:
456 explicit KrylovSolver(MPI_Comm comm);
457
461 KrylovSolver(KSP ksp, bool inc_ref_count);
462
463 // Copy constructor (deleted)
464 KrylovSolver(const KrylovSolver& solver) = delete;
465
467 KrylovSolver(KrylovSolver&& solver);
468
471
472 // Assignment operator (deleted)
473 KrylovSolver& operator=(const KrylovSolver&) = delete;
474
476 KrylovSolver& operator=(KrylovSolver&& solver);
477
479 void set_operator(const Mat A);
480
482 void set_operators(const Mat A, const Mat P);
483
486 int solve(Vec x, const Vec b, bool transpose = false) const;
487
490 void set_options_prefix(std::string options_prefix);
491
494 std::string get_options_prefix() const;
495
497 void set_from_options() const;
498
500 KSP ksp() const;
501
502private:
503 // PETSc solver pointer
504 KSP _ksp;
505};
506} // namespace petsc
507} // namespace dolfinx::la
508
509#endif
Definition IndexMap.h:94
Definition SparsityPattern.h:26
Definition Vector.h:32
constexpr int bs() const
Get block size.
Definition Vector.h:187
std::span< const value_type > array() const
Get local part of the vector (const version)
Definition Vector.h:190
std::shared_ptr< const common::IndexMap > index_map() const
Get IndexMap.
Definition Vector.h:184
Definition petsc.h:452
void set_operator(const Mat A)
Set operator (Mat)
Definition petsc.cpp:677
KSP ksp() const
Return PETSc KSP pointer.
Definition petsc.cpp:793
std::string get_options_prefix() const
Definition petsc.cpp:775
void set_operators(const Mat A, const Mat P)
Set operator and preconditioner matrix (Mat)
Definition petsc.cpp:679
KrylovSolver(MPI_Comm comm)
Definition petsc.cpp:640
void set_from_options() const
Set options from PETSc options database.
Definition petsc.cpp:785
~KrylovSolver()
Destructor.
Definition petsc.cpp:665
int solve(Vec x, const Vec b, bool transpose=false) const
Definition petsc.cpp:688
void set_options_prefix(std::string options_prefix)
Definition petsc.cpp:766
Definition petsc.h:285
static auto set_block_fn(Mat A, InsertMode mode)
Definition petsc.h:325
void set_from_options()
Call PETSc function MatSetFromOptions on the PETSc Mat object.
Definition petsc.cpp:633
double norm(Norm norm_type) const
Return norm of matrix.
Definition petsc.cpp:577
static auto set_fn(Mat A, InsertMode mode)
Definition petsc.h:291
Matrix(MPI_Comm comm, const SparsityPattern &sp, std::string type=std::string())
Create holder for a PETSc Mat object from a sparsity pattern.
Definition petsc.cpp:565
Matrix(Matrix &&A)=default
Move constructor (falls through to base class move constructor)
std::string get_options_prefix() const
Definition petsc.cpp:625
Matrix & operator=(Matrix &&A)=default
Move assignment operator.
static auto set_block_expand_fn(Mat A, int bs0, int bs1, InsertMode mode)
Definition petsc.h:362
~Matrix()=default
Destructor.
Matrix & operator=(const Matrix &A)=delete
Assignment operator (deleted)
AssemblyType
Definition petsc.h:420
void apply(AssemblyType type)
Definition petsc.cpp:602
void set_options_prefix(std::string options_prefix)
Definition petsc.cpp:619
Definition petsc.h:240
std::array< std::int64_t, 2 > size() const
Definition petsc.cpp:523
Operator & operator=(const Operator &A)=delete
Assignment operator (deleted)
Operator(Mat A, bool inc_ref_count)
Constructor.
Definition petsc.cpp:498
Vec create_vector(std::size_t dim) const
Definition petsc.cpp:533
Mat mat() const
Return PETSc Mat pointer.
Definition petsc.cpp:562
virtual ~Operator()
Destructor.
Definition petsc.cpp:509
Definition petsc.h:167
void set_from_options()
Call PETSc function VecSetFromOptions on the underlying Vec object.
Definition petsc.cpp:488
std::int64_t size() const
Return global size of the vector.
Definition petsc.cpp:435
std::string get_options_prefix() const
Definition petsc.cpp:479
Vector copy() const
Definition petsc.cpp:425
virtual ~Vector()
Destructor.
Definition petsc.cpp:413
Vector(const common::IndexMap &map, int bs)
Definition petsc.cpp:398
MPI_Comm comm() const
Return MPI communicator.
Definition petsc.cpp:463
std::array< std::int64_t, 2 > local_range() const
Return ownership range for calling rank.
Definition petsc.cpp:453
std::int32_t local_size() const
Return local size of vector (belonging to the call rank)
Definition petsc.cpp:444
Vec vec() const
Return pointer to PETSc Vec object.
Definition petsc.cpp:495
void set_options_prefix(std::string options_prefix)
Sets the prefix used by PETSc when searching the options database.
Definition petsc.cpp:472
Miscellaneous classes, functions and types.
Definition dolfinx_common.h:8
void clear()
Clear PETSc global options database.
Definition petsc.cpp:390
void set(std::string option)
Set PETSc option that takes no value.
Definition petsc.cpp:374
Vec create_vector_wrap(const common::IndexMap &map, int bs, std::span< const PetscScalar > x)
Definition petsc.cpp:102
MatNullSpace create_nullspace(MPI_Comm comm, std::span< const Vec > basis)
Definition petsc.cpp:362
void scatter_local_vectors(Vec x, const std::vector< std::span< const PetscScalar > > &x_b, const std::vector< std::pair< std::reference_wrapper< const common::IndexMap >, int > > &maps)
Scatter local vectors to Vec.
Definition petsc.cpp:192
void error(int error_code, std::string filename, std::string petsc_function)
Print error message for PETSc calls that return an error.
Definition petsc.cpp:34
std::vector< IS > create_index_sets(const std::vector< std::pair< std::reference_wrapper< const common::IndexMap >, int > > &maps)
Compute PETSc IndexSets (IS) for a stack of index maps.
Definition petsc.cpp:129
std::vector< Vec > create_vectors(MPI_Comm comm, const std::vector< std::span< const PetscScalar > > &x)
Definition petsc.cpp:52
Mat create_matrix(MPI_Comm comm, const SparsityPattern &sp, std::string type=std::string())
Definition petsc.cpp:235
std::vector< std::vector< PetscScalar > > get_local_vectors(const Vec x, const std::vector< std::pair< std::reference_wrapper< const common::IndexMap >, int > > &maps)
Copy blocks from Vec into local vectors.
Definition petsc.cpp:149
Vec create_vector(const common::IndexMap &map, int bs)
Definition petsc.cpp:68
Linear algebra interface.
Definition sparsitybuild.h:15
Norm
Norm types.
Definition utils.h:17