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_LASPACK_MATRIX_H 21 #define LIBMESH_LASPACK_MATRIX_H 22 23 #include "libmesh/libmesh_config.h" 24 25 #ifdef LIBMESH_HAVE_LASPACK 26 27 // Local includes 28 #include "libmesh/sparse_matrix.h" 29 30 // Laspack includes 31 #include <qmatrix.h> 32 33 // C++ includes 34 #include <algorithm> 35 #include <cstddef> 36 37 namespace libMesh 38 { 39 40 // Forward declarations 41 template <typename T> class DenseMatrix; 42 template <typename T> class LaspackVector; 43 template <typename T> class LaspackLinearSolver; 44 45 /** 46 * The LaspackMatrix class wraps a QMatrix object from the Laspack 47 * library. Currently Laspack only supports real datatypes, so this 48 * class is a full specialization of \p SparseMatrix<T> with \p T = \p 49 * Real. All overridden virtual functions are documented in 50 * sparse_matrix.h. 51 * 52 * \author Benjamin S. Kirk 53 * \date 2003 54 */ 55 template <typename T> 56 class LaspackMatrix final : public SparseMatrix<T> 57 { 58 59 public: 60 /** 61 * Constructor; initializes the matrix to be empty, without any 62 * structure, i.e. the matrix is not usable at all. This 63 * constructor is therefore only useful for matrices which are 64 * members of a class. All other matrices should be created at a 65 * point in the data flow where all necessary information is 66 * available. 67 * 68 * You have to initialize the matrix before usage with \p init(...). 69 */ 70 LaspackMatrix (const Parallel::Communicator & comm); 71 72 /** 73 * This class manages a C-style struct (QMatrix) manually, so we 74 * don't want to allow any automatic copy/move functions to be 75 * generated, and we can't default the destructor. 76 */ 77 LaspackMatrix (LaspackMatrix &&) = delete; 78 LaspackMatrix (const LaspackMatrix &) = delete; 79 LaspackMatrix & operator= (const LaspackMatrix &) = delete; 80 LaspackMatrix & operator= (LaspackMatrix &&) = delete; 81 virtual ~LaspackMatrix (); 82 83 /** 84 * The \p LaspackMatrix needs the full sparsity pattern. 85 */ need_full_sparsity_pattern()86 virtual bool need_full_sparsity_pattern() const override 87 { return true; } 88 89 /** 90 * Updates the matrix sparsity pattern. This will tell the 91 * underlying matrix storage scheme how to map the \f$ (i,j) \f$ 92 * elements. 93 */ 94 virtual void update_sparsity_pattern (const SparsityPattern::Graph &) override; 95 96 virtual void init (const numeric_index_type m, 97 const numeric_index_type n, 98 const numeric_index_type m_l, 99 const numeric_index_type n_l, 100 const numeric_index_type nnz=30, 101 const numeric_index_type noz=10, 102 const numeric_index_type blocksize=1) override; 103 104 virtual void init (ParallelType = PARALLEL) override; 105 106 virtual void clear () override; 107 108 virtual void zero () override; 109 110 virtual std::unique_ptr<SparseMatrix<T>> zero_clone () const override; 111 112 virtual std::unique_ptr<SparseMatrix<T>> clone () const override; 113 114 virtual void close () override; 115 116 virtual numeric_index_type m () const override; 117 118 virtual numeric_index_type n () const override; 119 120 virtual numeric_index_type row_start () const override; 121 122 virtual numeric_index_type row_stop () const override; 123 124 virtual void set (const numeric_index_type i, 125 const numeric_index_type j, 126 const T value) override; 127 128 virtual void add (const numeric_index_type i, 129 const numeric_index_type j, 130 const T value) override; 131 132 virtual void add_matrix (const DenseMatrix<T> & dm, 133 const std::vector<numeric_index_type> & rows, 134 const std::vector<numeric_index_type> & cols) override; 135 136 virtual void add_matrix (const DenseMatrix<T> & dm, 137 const std::vector<numeric_index_type> & dof_indices) override; 138 139 /** 140 * Compute A += a*X for scalar \p a, matrix \p X. 141 * 142 * \note LASPACK does not provide a true \p axpy for matrices, 143 * so a hand-coded version with hopefully acceptable performance 144 * is provided. 145 */ 146 virtual void add (const T a, const SparseMatrix<T> & X) override; 147 148 virtual T operator () (const numeric_index_type i, 149 const numeric_index_type j) const override; 150 l1_norm()151 virtual Real l1_norm () const override { libmesh_not_implemented(); return 0.; } 152 linfty_norm()153 virtual Real linfty_norm () const override { libmesh_not_implemented(); return 0.; } 154 closed()155 virtual bool closed() const override { return _closed; } 156 157 virtual void print_personal(std::ostream & os=libMesh::out) const override { this->print(os); } 158 159 virtual void get_diagonal (NumericVector<T> & dest) const override; 160 161 virtual void get_transpose (SparseMatrix<T> & dest) const override; 162 163 private: 164 165 /** 166 * \returns The position in the compressed row 167 * storage scheme of the \f$ (i,j) \f$ element. 168 */ 169 numeric_index_type pos (const numeric_index_type i, 170 const numeric_index_type j) const; 171 172 /** 173 * The Laspack sparse matrix pointer. 174 */ 175 QMatrix _QMat; 176 177 /** 178 * The compressed row indices. 179 */ 180 std::vector<numeric_index_type> _csr; 181 182 /** 183 * The start of each row in the compressed 184 * row index data structure. 185 */ 186 std::vector<std::vector<numeric_index_type>::const_iterator> _row_start; 187 188 /** 189 * Flag indicating if the matrix has been closed yet. 190 */ 191 bool _closed; 192 193 /** 194 * Make other Laspack datatypes friends 195 */ 196 friend class LaspackVector<T>; 197 friend class LaspackLinearSolver<T>; 198 }; 199 200 } // namespace libMesh 201 202 #endif // #ifdef LIBMESH_HAVE_LASPACK 203 #endif // #ifdef LIBMESH_LASPACK_MATRIX_H 204