update

2025.09.22_cpp_with_eigen_package/Eigen/src/MetisSupport/InternalHeaderCheck.h

@@ -0,0 +1,3 @@
+#ifndef EIGEN_METISSUPPORT_MODULE_H
+#error "Please include Eigen/MetisSupport instead of including headers inside the src directory directly."
+#endif

2025.09.22_cpp_with_eigen_package/Eigen/src/MetisSupport/MetisSupport.h

@@ -0,0 +1,125 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef METIS_SUPPORT_H
+#define METIS_SUPPORT_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+/**
+ * Get the fill-reducing ordering from the METIS package
+ *
+ * If A is the original matrix and Ap is the permuted matrix,
+ * the fill-reducing permutation is defined as follows :
+ * Row (column) i of A is the matperm(i) row (column) of Ap.
+ * WARNING: As computed by METIS, this corresponds to the vector iperm (instead of perm)
+ */
+template <typename StorageIndex>
+class MetisOrdering {
+ public:
+  typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
+  typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
+
+  template <typename MatrixType>
+  void get_symmetrized_graph(const MatrixType& A) {
+    Index m = A.cols();
+    eigen_assert((A.rows() == A.cols()) && "ONLY FOR SQUARED MATRICES");
+    // Get the transpose of the input matrix
+    MatrixType At = A.transpose();
+    // Get the number of nonzeros elements in each row/col of At+A
+    Index TotNz = 0;
+    IndexVector visited(m);
+    visited.setConstant(-1);
+    for (StorageIndex j = 0; j < m; j++) {
+      // Compute the union structure of of A(j,:) and At(j,:)
+      visited(j) = j;  // Do not include the diagonal element
+      // Get the nonzeros in row/column j of A
+      for (typename MatrixType::InnerIterator it(A, j); it; ++it) {
+        Index idx = it.index();  // Get the row index (for column major) or column index (for row major)
+        if (visited(idx) != j) {
+          visited(idx) = j;
+          ++TotNz;
+        }
+      }
+      // Get the nonzeros in row/column j of At
+      for (typename MatrixType::InnerIterator it(At, j); it; ++it) {
+        Index idx = it.index();
+        if (visited(idx) != j) {
+          visited(idx) = j;
+          ++TotNz;
+        }
+      }
+    }
+    // Reserve place for A + At
+    m_indexPtr.resize(m + 1);
+    m_innerIndices.resize(TotNz);
+
+    // Now compute the real adjacency list of each column/row
+    visited.setConstant(-1);
+    StorageIndex CurNz = 0;
+    for (StorageIndex j = 0; j < m; j++) {
+      m_indexPtr(j) = CurNz;
+
+      visited(j) = j;  // Do not include the diagonal element
+      // Add the pattern of row/column j of A to A+At
+      for (typename MatrixType::InnerIterator it(A, j); it; ++it) {
+        StorageIndex idx = it.index();  // Get the row index (for column major) or column index (for row major)
+        if (visited(idx) != j) {
+          visited(idx) = j;
+          m_innerIndices(CurNz) = idx;
+          CurNz++;
+        }
+      }
+      // Add the pattern of row/column j of At to A+At
+      for (typename MatrixType::InnerIterator it(At, j); it; ++it) {
+        StorageIndex idx = it.index();
+        if (visited(idx) != j) {
+          visited(idx) = j;
+          m_innerIndices(CurNz) = idx;
+          ++CurNz;
+        }
+      }
+    }
+    m_indexPtr(m) = CurNz;
+  }
+
+  template <typename MatrixType>
+  void operator()(const MatrixType& A, PermutationType& matperm) {
+    StorageIndex m = internal::convert_index<StorageIndex>(
+        A.cols());  // must be StorageIndex, because it is passed by address to METIS
+    IndexVector perm(m), iperm(m);
+    // First, symmetrize the matrix graph.
+    get_symmetrized_graph(A);
+    int output_error;
+
+    // Call the fill-reducing routine from METIS
+    output_error = METIS_NodeND(&m, m_indexPtr.data(), m_innerIndices.data(), NULL, NULL, perm.data(), iperm.data());
+
+    if (output_error != METIS_OK) {
+      // FIXME The ordering interface should define a class of possible errors
+      std::cerr << "ERROR WHILE CALLING THE METIS PACKAGE \n";
+      return;
+    }
+
+    // Get the fill-reducing permutation
+    // NOTE:  If Ap is the permuted matrix then perm and iperm vectors are defined as follows
+    // Row (column) i of Ap is the perm(i) row(column) of A, and row (column) i of A is the iperm(i) row(column) of Ap
+
+    matperm.resize(m);
+    for (int j = 0; j < m; j++) matperm.indices()(iperm(j)) = j;
+  }
+
+ protected:
+  IndexVector m_indexPtr;      // Pointer to the adjacenccy list of each row/column
+  IndexVector m_innerIndices;  // Adjacency list
+};
+
+}  // namespace Eigen
+#endif