1 // The libMesh Finite Element Library. 2 // Copyright (C) 2002-2020 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 3 4 // This library is free software; you can redistribute it and/or 5 // modify it under the terms of the GNU Lesser General Public 6 // License as published by the Free Software Foundation; either 7 // version 2.1 of the License, or (at your option) any later version. 8 9 // This library is distributed in the hope that it will be useful, 10 // but WITHOUT ANY WARRANTY; without even the implied warranty of 11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12 // Lesser General Public License for more details. 13 14 // You should have received a copy of the GNU Lesser General Public 15 // License along with this library; if not, write to the Free Software 16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 17 18 19 20 #ifndef LIBMESH_PETSC_MATRIX_H 21 #define LIBMESH_PETSC_MATRIX_H 22 23 #include "libmesh/libmesh_common.h" 24 25 #ifdef LIBMESH_HAVE_PETSC 26 27 // Local includes 28 #include "libmesh/sparse_matrix.h" 29 #include "libmesh/petsc_macro.h" 30 31 // C++ includes 32 #include <algorithm> 33 #ifdef LIBMESH_HAVE_CXX11_THREAD 34 #include <atomic> 35 #include <mutex> 36 #endif 37 38 // Macro to identify and debug functions which should be called in 39 // parallel on parallel matrices but which may be called in serial on 40 // serial matrices. This macro will only be valid inside non-static 41 // PetscMatrix methods 42 #undef semiparallel_only 43 #ifndef NDEBUG 44 #include <cstring> 45 46 #define semiparallel_only() do { if (this->initialized()) { const char * mytype; \ 47 MatGetType(_mat,&mytype); \ 48 if (!strcmp(mytype, MATSEQAIJ)) \ 49 parallel_object_only(); } } while (0) 50 #else 51 #define semiparallel_only() 52 #endif 53 54 55 // Petsc include files. 56 #ifdef I 57 # define LIBMESH_SAW_I 58 #endif 59 #include <petscmat.h> 60 #ifndef LIBMESH_SAW_I 61 # undef I // Avoid complex.h contamination 62 #endif 63 64 65 66 namespace libMesh 67 { 68 69 // Forward Declarations 70 template <typename T> class DenseMatrix; 71 72 enum PetscMatrixType : int { 73 AIJ=0, 74 HYPRE}; 75 76 77 /** 78 * This class provides a nice interface to the PETSc C-based data 79 * structures for parallel, sparse matrices. All overridden virtual 80 * functions are documented in sparse_matrix.h. 81 * 82 * \author Benjamin S. Kirk 83 * \date 2002 84 * \brief SparseMatrix interface to PETSc Mat. 85 */ 86 template <typename T> 87 class PetscMatrix final : public SparseMatrix<T> 88 { 89 public: 90 /** 91 * Constructor; initializes the matrix to be empty, without any 92 * structure, i.e. the matrix is not usable at all. This 93 * constructor is therefore only useful for matrices which are 94 * members of a class. All other matrices should be created at a 95 * point in the data flow where all necessary information is 96 * available. 97 * 98 * You have to initialize the matrix before usage with \p init(...). 99 */ 100 explicit 101 PetscMatrix (const Parallel::Communicator & comm_in); 102 103 /** 104 * Constructor. Creates a PetscMatrix assuming you already have a 105 * valid Mat object. In this case, m is NOT destroyed by the 106 * PetscMatrix destructor when this object goes out of scope. This 107 * allows ownership of m to remain with the original creator, and to 108 * simply provide additional functionality with the PetscMatrix. 109 */ 110 explicit 111 PetscMatrix (Mat m, 112 const Parallel::Communicator & comm_in); 113 114 /** 115 * This class manages a C-style struct (Mat) manually, so we 116 * don't want to allow any automatic copy/move functions to be 117 * generated, and we can't default the destructor. 118 */ 119 PetscMatrix (PetscMatrix &&) = delete; 120 PetscMatrix (const PetscMatrix &) = delete; 121 PetscMatrix & operator= (const PetscMatrix &) = delete; 122 PetscMatrix & operator= (PetscMatrix &&) = delete; 123 virtual ~PetscMatrix (); 124 125 void set_matrix_type(PetscMatrixType mat_type); 126 127 virtual void init (const numeric_index_type m, 128 const numeric_index_type n, 129 const numeric_index_type m_l, 130 const numeric_index_type n_l, 131 const numeric_index_type nnz=30, 132 const numeric_index_type noz=10, 133 const numeric_index_type blocksize=1) override; 134 135 /** 136 * Initialize a PETSc matrix. 137 * 138 * \param m The global number of rows. 139 * \param n The global number of columns. 140 * \param m_l The local number of rows. 141 * \param n_l The local number of columns. 142 * \param n_nz array containing the number of nonzeros in each row of the DIAGONAL portion of the local submatrix. 143 * \param n_oz Array containing the number of nonzeros in each row of the OFF-DIAGONAL portion of the local submatrix. 144 * \param blocksize Optional value indicating dense coupled blocks for systems with multiple variables all of the same type. 145 */ 146 void init (const numeric_index_type m, 147 const numeric_index_type n, 148 const numeric_index_type m_l, 149 const numeric_index_type n_l, 150 const std::vector<numeric_index_type> & n_nz, 151 const std::vector<numeric_index_type> & n_oz, 152 const numeric_index_type blocksize=1); 153 154 virtual void init (ParallelType = PARALLEL) override; 155 156 /** 157 * Update the sparsity pattern based on \p dof_map, and set the matrix 158 * to zero. This is useful in cases where the sparsity pattern changes 159 * during a computation. 160 */ 161 void update_preallocation_and_zero(); 162 163 /** 164 * Reset matrix to use the original nonzero pattern provided by users 165 */ 166 void reset_preallocation(); 167 168 virtual void clear () override; 169 170 virtual void zero () override; 171 172 virtual std::unique_ptr<SparseMatrix<T>> zero_clone () const override; 173 174 virtual std::unique_ptr<SparseMatrix<T>> clone () const override; 175 176 virtual void zero_rows (std::vector<numeric_index_type> & rows, T diag_value = 0.0) override; 177 178 virtual void close () override; 179 180 virtual void flush () override; 181 182 virtual numeric_index_type m () const override; 183 184 numeric_index_type local_m () const final; 185 186 virtual numeric_index_type n () const override; 187 188 /** 189 * Get the number of columns owned by this process 190 */ 191 numeric_index_type local_n () const; 192 193 /** 194 * Get the number of rows and columns owned by this process 195 * @param i Row size 196 * @param j Column size 197 */ 198 void get_local_size (numeric_index_type & m, numeric_index_type & n) const; 199 200 virtual numeric_index_type row_start () const override; 201 202 virtual numeric_index_type row_stop () const override; 203 204 virtual void set (const numeric_index_type i, 205 const numeric_index_type j, 206 const T value) override; 207 208 virtual void add (const numeric_index_type i, 209 const numeric_index_type j, 210 const T value) override; 211 212 virtual void add_matrix (const DenseMatrix<T> & dm, 213 const std::vector<numeric_index_type> & rows, 214 const std::vector<numeric_index_type> & cols) override; 215 216 virtual void add_matrix (const DenseMatrix<T> & dm, 217 const std::vector<numeric_index_type> & dof_indices) override; 218 219 virtual void add_block_matrix (const DenseMatrix<T> & dm, 220 const std::vector<numeric_index_type> & brows, 221 const std::vector<numeric_index_type> & bcols) override; 222 add_block_matrix(const DenseMatrix<T> & dm,const std::vector<numeric_index_type> & dof_indices)223 virtual void add_block_matrix (const DenseMatrix<T> & dm, 224 const std::vector<numeric_index_type> & dof_indices) override 225 { this->add_block_matrix (dm, dof_indices, dof_indices); } 226 227 /** 228 * Compute A += a*X for scalar \p a, matrix \p X. 229 * 230 * \note The matrices A and X need to have the same nonzero pattern, 231 * otherwise \p PETSc will crash! 232 * 233 * \note It is advisable to not only allocate appropriate memory 234 * with \p init(), but also explicitly zero the terms of \p this 235 * whenever you add a non-zero value to \p X. 236 * 237 * \note \p X will be closed, if not already done, before performing 238 * any work. 239 */ 240 virtual void add (const T a, const SparseMatrix<T> & X) override; 241 242 /** 243 * Compute Y = A*X for matrix \p X. 244 */ 245 virtual void matrix_matrix_mult (SparseMatrix<T> & X, SparseMatrix<T> & Y) override; 246 247 /** 248 * Add \p scalar* \p spm to the rows and cols of this matrix (A): 249 * A(rows[i], cols[j]) += scalar * spm(i,j) 250 */ 251 void add_sparse_matrix (const SparseMatrix<T> & spm, 252 const std::map<numeric_index_type,numeric_index_type> & row_ltog, 253 const std::map<numeric_index_type,numeric_index_type> & col_ltog, 254 const T scalar) override; 255 256 virtual T operator () (const numeric_index_type i, 257 const numeric_index_type j) const override; 258 259 virtual Real l1_norm () const override; 260 261 virtual Real linfty_norm () const override; 262 263 virtual bool closed() const override; 264 265 /** 266 * If set to false, we don't delete the Mat on destruction and allow 267 * instead for \p PETSc to manage it. 268 */ 269 virtual void set_destroy_mat_on_exit(bool destroy = true); 270 271 /** 272 * Print the contents of the matrix to the screen with the PETSc 273 * viewer. This function only allows printing to standard out since 274 * we have limited ourselves to one PETSc implementation for 275 * writing. 276 */ 277 virtual void print_personal(std::ostream & os=libMesh::out) const override; 278 279 virtual void print_matlab(const std::string & name = "") const override; 280 281 virtual void get_diagonal (NumericVector<T> & dest) const override; 282 283 virtual void get_transpose (SparseMatrix<T> & dest) const override; 284 285 /** 286 * Swaps the internal data pointers of two PetscMatrices, no actual 287 * values are swapped. 288 */ 289 void swap (PetscMatrix<T> &); 290 291 /** 292 * \returns The raw PETSc matrix pointer. 293 * 294 * \note This is generally not required in user-level code. 295 * 296 * \note Don't do anything crazy like calling MatDestroy() on 297 * it, or very bad things will likely happen! 298 */ mat()299 Mat mat () { libmesh_assert (_mat); return _mat; } 300 301 void get_row(numeric_index_type i, 302 std::vector<numeric_index_type> & indices, 303 std::vector<T> & values) const override; 304 protected: 305 306 /** 307 * This function either creates or re-initializes a matrix called \p 308 * submatrix which is defined by the indices given in the \p rows 309 * and \p cols vectors. 310 * 311 * This function is implemented in terms of MatGetSubMatrix(). The 312 * \p reuse_submatrix parameter determines whether or not PETSc will 313 * treat \p submatrix as one which has already been used (had memory 314 * allocated) or as a new matrix. 315 */ 316 virtual void _get_submatrix(SparseMatrix<T> & submatrix, 317 const std::vector<numeric_index_type> & rows, 318 const std::vector<numeric_index_type> & cols, 319 const bool reuse_submatrix) const override; 320 321 private: 322 323 /** 324 * PETSc matrix datatype to store values. 325 */ 326 Mat _mat; 327 328 /** 329 * This boolean value should only be set to \p false for the 330 * constructor which takes a PETSc Mat object. 331 */ 332 bool _destroy_mat_on_exit; 333 334 PetscMatrixType _mat_type; 335 336 #ifdef LIBMESH_HAVE_CXX11_THREAD 337 mutable std::mutex _petsc_matrix_mutex; 338 #else 339 mutable Threads::spin_mutex _petsc_matrix_mutex; 340 #endif 341 }; 342 343 } // namespace libMesh 344 345 #endif // #ifdef LIBMESH_HAVE_PETSC 346 #endif // LIBMESH_PETSC_MATRIX_H 347