diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
index 013c6ff66c50c922f38862d6a78e54908fd05b40..3b788ee481b6deea5a27973b454472b49ff23cfa 100644
--- a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
@@ -38,7 +38,7 @@ PETScLinearSolver::PETScLinearSolver(const std::string /*prefix*/,
         if (auto const subtree = option->getConfigSubtreeOptional("petsc"))
         {
             if (auto const parameters =
-                //! \ogs_file_param{prj__linear_solvers__linear_solver__petsc__parameters}
+                    //! \ogs_file_param{prj__linear_solvers__linear_solver__petsc__parameters}
                 subtree->getConfigParameterOptional<std::string>("parameters"))
             {
                 petsc_options = *parameters;
@@ -46,7 +46,7 @@ PETScLinearSolver::PETScLinearSolver(const std::string /*prefix*/,
 
             if (auto const pre =
                     //! \ogs_file_param{prj__linear_solvers__linear_solver__petsc__prefix}
-                    subtree->getConfigParameterOptional<std::string>("prefix"))
+                subtree->getConfigParameterOptional<std::string>("prefix"))
             {
                 if (!pre->empty())
                     prefix = *pre + "_";
@@ -118,7 +118,7 @@ bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector& b, PETScVector& x)
                 PetscPrintf(PETSC_COMM_WORLD,
                             " (relative convergence criterion fulfilled).");
                 break;
-            case  KSP_CONVERGED_ATOL:
+            case KSP_CONVERGED_ATOL:
                 PetscPrintf(PETSC_COMM_WORLD,
                             " (absolute convergence criterion fulfilled).");
                 break;
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.h b/MathLib/LinAlg/PETSc/PETScMatrix.h
index bca4e15069108226537647faae9ded15ae99c436..2b29072f6a09ad28078fc76397b7c2620829b403 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.h
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.h
@@ -127,11 +127,30 @@ public:
         MatSetValue(_A, i, j, value, ADD_VALUES);
     }
 
-    /// Add sub-matrix at positions given by \c indices.
+    /*!
+       \brief Add sub-matrix at positions given by global \c indices, in which
+       negative index indicates ghost entry.
+
+       In order to use MatZeroRows to apply Dirichlet boundary condition,
+       entries in the rows with the negative global indices are skipped to added
+       to the global matrix, meanwhile entries in the columns with the negative
+       global indices are added the global matrix. By using MatZeroRows to apply
+       Dirichlet boundary condition, the off diagonal entries in ghost rows of
+       the global matrix are set to zero, while the off diagonal entries in
+       ghost rows of the global matrix are assembled and kept for linear solver.
+
+       For the setting of Dirichlet boundary condition
+       in PETSc, please refer to
+       <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html#redistribute">PETSc:Documentation:FAQ</a>
+     */
     template <class T_DENSE_MATRIX>
     void add(RowColumnIndices<PetscInt> const& indices,
              const T_DENSE_MATRIX& sub_matrix)
     {
+        // Set global column indices to positive to allow all entries of columns
+        // to be added to the global matrix. For the ghost columns, only the
+        // off diagonal entries are added due to the negative indices of the
+        // corresponding rows.
         std::vector<PetscInt> cols;
         cols.reserve(indices.columns.size());
         for (auto col : indices.columns)