// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced // at the Lawrence Livermore National Laboratory. All Rights reserved. See files // LICENSE and NOTICE for details. LLNL-CODE-806117. // // This file is part of the MFEM library. For more information and source code // availability visit https://mfem.org. // // MFEM is free software; you can redistribute it and/or modify it under the // terms of the BSD-3 license. We welcome feedback and contributions, see file // CONTRIBUTING.md for details. #ifndef MFEM_BLOCKMATRIX #define MFEM_BLOCKMATRIX #include "../config/config.hpp" #include "../general/array.hpp" #include "../general/globals.hpp" #include "vector.hpp" #include "sparsemat.hpp" namespace mfem { class BlockMatrix : public AbstractSparseMatrix { public: //! Constructor for square block matrices /** @param offsets offsets that mark the start of each row/column block (size nRowBlocks+1). @note BlockMatrix will not own/copy the data contained in offsets. */ BlockMatrix(const Array & offsets); //! Constructor for rectangular block matrices /** @param row_offsets offsets that mark the start of each row block (size nRowBlocks+1). @param col_offsets offsets that mark the start of each column block (size nColBlocks+1). @note BlockMatrix will not own/copy the data contained in offsets. */ BlockMatrix(const Array & row_offsets, const Array & col_offsets); //! Set A(i,j) = mat void SetBlock(int i, int j, SparseMatrix * mat); //! Return the number of row blocks int NumRowBlocks() const {return nRowBlocks; } //! Return the number of column blocks int NumColBlocks() const {return nColBlocks; } //! Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL SparseMatrix & GetBlock(int i, int j); //! Return a reference to block (i,j). Reference may be invalid if Aij(i,j) //! == NULL. (const version) const SparseMatrix & GetBlock(int i, int j) const; //! Check if block (i,j) is a zero block. int IsZeroBlock(int i, int j) const {return (Aij(i,j)==NULL) ? 1 : 0; } //! Return the row offsets for block starts Array & RowOffsets() { return row_offsets; } //! Return the columns offsets for block starts Array & ColOffsets() { return col_offsets; } //! Return the row offsets for block starts (const version) const Array & RowOffsets() const { return row_offsets; } //! Return the row offsets for block starts (const version) const Array & ColOffsets() const { return col_offsets; } //! Return the number of non zeros in row i int RowSize(const int i) const; /// Eliminate the row and column @a rc from the matrix. /** Eliminates the column and row @a rc, replacing the element (rc,rc) with 1.0. Assumes that element (i,rc) is assembled if and only if the element (rc,i) is assembled. If @a dpolicy is specified, the element (rc,rc) is treated according to that policy. */ void EliminateRowCol(int rc, DiagonalPolicy dpolicy = DIAG_ONE); /** @brief Eliminate the rows and columns corresponding to the entries in @a vdofs + save the eliminated entries into @a Ae so that (*this) + Ae is equal to the original matrix. */ void EliminateRowCols(const Array & vdofs, BlockMatrix *Ae, DiagonalPolicy dpolicy = DIAG_ONE); //! Symmetric elimination of the marked degree of freedom. /** @param ess_bc_dofs marker of the degree of freedom to be eliminated dof i is eliminated if @a ess_bc_dofs[i] = 1. @param sol vector that stores the values of the degree of freedom that need to be eliminated @param rhs vector that stores the rhs of the system. */ void EliminateRowCol(Array & ess_bc_dofs, Vector & sol, Vector & rhs); /// Finalize all the submatrices virtual void Finalize(int skip_zeros = 1) { Finalize(skip_zeros, false); } /// A slightly more general version of the Finalize(int) method. void Finalize(int skip_zeros, bool fix_empty_rows); //! Returns a monolithic CSR matrix that represents this operator. SparseMatrix * CreateMonolithic() const; //! Export the monolithic matrix to file. virtual void PrintMatlab(std::ostream & os = mfem::out) const; /// @name Matrix interface ///@{ /// Returns reference to a_{ij}. virtual double& Elem (int i, int j); /// Returns constant reference to a_{ij}. virtual const double& Elem (int i, int j) const; /// Returns a pointer to (approximation) of the matrix inverse. virtual MatrixInverse * Inverse() const { mfem_error("BlockMatrix::Inverse not implemented \n"); return static_cast(NULL); } ///@} ///@name AbstractSparseMatrix interface ///@{ //! Returns the total number of non zeros in the matrix. virtual int NumNonZeroElems() const; /// Gets the columns indexes and values for row *row*. /** The return value is always 0 since @a cols and @a srow are copies of the values in the matrix. */ virtual int GetRow(const int row, Array &cols, Vector &srow) const; /** @brief If the matrix is square, this method will place 1 on the diagonal (i,i) if row i has "almost" zero l1-norm. If entry (i,i) does not belong to the sparsity pattern of A, then a error will occur. */ virtual void EliminateZeroRows(const double threshold = 1e-12); /// Matrix-Vector Multiplication y = A*x virtual void Mult(const Vector & x, Vector & y) const; /// Matrix-Vector Multiplication y = y + val*A*x virtual void AddMult(const Vector & x, Vector & y, const double val = 1.) const; /// MatrixTranspose-Vector Multiplication y = A'*x virtual void MultTranspose(const Vector & x, Vector & y) const; /// MatrixTranspose-Vector Multiplication y = y + val*A'*x virtual void AddMultTranspose(const Vector & x, Vector & y, const double val = 1.) const; ///@} /** @brief Partial matrix vector multiplication of (*this) with @a x involving only the rows given by @a rows. The result is given in @a y */ void PartMult(const Array &rows, const Vector &x, Vector &y) const; /** @brief Partial matrix vector multiplication of (*this) with @a x involving only the rows given by @a rows. The result is multiplied by @a a and added to @a y */ void PartAddMult(const Array &rows, const Vector &x, Vector &y, const double a=1.0) const; //! Destructor virtual ~BlockMatrix(); //! If owns_blocks the SparseMatrix objects Aij will be deallocated. int owns_blocks; virtual Type GetType() const { return MFEM_Block_Matrix; } private: //! Given a global row iglobal finds to which row iloc in block iblock belongs to. inline void findGlobalRow(int iglobal, int & iblock, int & iloc) const; //! Given a global column jglobal finds to which column jloc in block jblock belongs to. inline void findGlobalCol(int jglobal, int & jblock, int & jloc) const; //! Number of row blocks int nRowBlocks; //! Number of columns blocks int nColBlocks; //! row offsets for each block start (length nRowBlocks+1). Array row_offsets; //! column offsets for each block start (length nColBlocks+1). Array col_offsets; //! 2D array that stores each block of the BlockMatrix. Aij(iblock, jblock) //! == NULL if block (iblock, jblock) is all zeros. Array2D Aij; }; //! Transpose a BlockMatrix: result = A' BlockMatrix * Transpose(const BlockMatrix & A); //! Multiply BlockMatrix matrices: result = A*B BlockMatrix * Mult(const BlockMatrix & A, const BlockMatrix & B); inline void BlockMatrix::findGlobalRow(int iglobal, int & iblock, int & iloc) const { if (iglobal > row_offsets[nRowBlocks]) { mfem_error("BlockMatrix::findGlobalRow"); } for (iblock = 0; iblock < nRowBlocks; ++iblock) { if (row_offsets[iblock+1] > iglobal) { break; } } iloc = iglobal - row_offsets[iblock]; } inline void BlockMatrix::findGlobalCol(int jglobal, int & jblock, int & jloc) const { if (jglobal > col_offsets[nColBlocks]) { mfem_error("BlockMatrix::findGlobalCol"); } for (jblock = 0; jblock < nColBlocks; ++jblock) { if (col_offsets[jblock+1] > jglobal) { break; } } jloc = jglobal - col_offsets[jblock]; } } #endif /* MFEM_BLOCKMATRIX */