MatrixMarketIterator.h

// SPDX-License-Identifier: MPL-2.0
 
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Desire 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 EIGEN_BROWSE_MATRICES_H
#define EIGEN_BROWSE_MATRICES_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
enum { SPD = 0x100, NonSymmetric = 0x0 };

template <typename Scalar>
class MatrixMarketIterator {
  typedef typename NumTraits<Scalar>::Real RealScalar;
 
 public:
  typedef Matrix<Scalar, Dynamic, 1> VectorType;
  typedef SparseMatrix<Scalar, ColMajor> MatrixType;
 
 public:
  MatrixMarketIterator(const std::string& folder)
      : m_sym(0), m_isvalid(false), m_matIsLoaded(false), m_hasRhs(false), m_hasrefX(false), m_folder(folder) {
    m_folder_id = opendir(folder.c_str());
    if (m_folder_id) Getnextvalidmatrix();
  }
 
  ~MatrixMarketIterator() {
    if (m_folder_id) closedir(m_folder_id);
  }
 
  inline MatrixMarketIterator& operator++() {
    m_matIsLoaded = false;
    m_hasrefX = false;
    m_hasRhs = false;
    Getnextvalidmatrix();
    return *this;
  }
  inline operator bool() const { return m_isvalid; }

  inline MatrixType& matrix() {
    // Read the matrix
    if (m_matIsLoaded) return m_mat;
 
    std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
    if (!loadMarket(m_mat, matrix_file)) {
      std::cerr << "Warning loadMarket failed when loading \"" << matrix_file << "\"" << std::endl;
      m_matIsLoaded = false;
      return m_mat;
    }
    m_matIsLoaded = true;
 
    if (m_sym != NonSymmetric) {
      // Check whether we need to restore a full matrix:
      RealScalar diag_norm = m_mat.diagonal().norm();
      RealScalar lower_norm = m_mat.template triangularView<Lower>().norm();
      RealScalar upper_norm = m_mat.template triangularView<Upper>().norm();
      if (lower_norm > diag_norm && upper_norm == diag_norm) {
        // only the lower part is stored
        MatrixType tmp(m_mat);
        m_mat = tmp.template selfadjointView<Lower>();
      } else if (upper_norm > diag_norm && lower_norm == diag_norm) {
        // only the upper part is stored
        MatrixType tmp(m_mat);
        m_mat = tmp.template selfadjointView<Upper>();
      }
    }
    return m_mat;
  }

  inline VectorType& rhs() {
    // Get the right hand side
    if (m_hasRhs) return m_rhs;
 
    std::string rhs_file;
    rhs_file = m_folder + "/" + m_matname + "_b.mtx";  // The pattern is matname_b.mtx
    m_hasRhs = Fileexists(rhs_file);
    if (m_hasRhs) {
      m_rhs.resize(m_mat.cols());
      m_hasRhs = loadMarketVector(m_rhs, rhs_file);
    }
    if (!m_hasRhs) {
      // Generate a random right hand side
      if (!m_matIsLoaded) this->matrix();
      m_refX.resize(m_mat.cols());
      m_refX.setRandom();
      m_rhs = m_mat * m_refX;
      m_hasrefX = true;
      m_hasRhs = true;
    }
    return m_rhs;
  }

  inline VectorType& refX() {
    // Check if a reference solution is provided
    if (m_hasrefX) return m_refX;
 
    std::string lhs_file;
    lhs_file = m_folder + "/" + m_matname + "_x.mtx";
    m_hasrefX = Fileexists(lhs_file);
    if (m_hasrefX) {
      m_refX.resize(m_mat.cols());
      m_hasrefX = loadMarketVector(m_refX, lhs_file);
    } else
      m_refX.resize(0);
    return m_refX;
  }
 
  inline std::string& matname() { return m_matname; }
 
  inline int sym() const { return m_sym; }
 
  bool hasRhs() const { return m_hasRhs; }
  bool hasrefX() const { return m_hasrefX; }
  bool isFolderValid() const { return bool(m_folder_id); }
 
 protected:
  inline bool Fileexists(std::string file) const {
    std::ifstream file_id(file.c_str());
    if (!file_id.good()) {
      return false;
    } else {
      file_id.close();
      return true;
    }
  }
 
  void Getnextvalidmatrix() {
    m_isvalid = false;
    // Here, we return with the next valid matrix in the folder
    while ((m_curs_id = readdir(m_folder_id)) != nullptr) {
      m_isvalid = false;
      std::string curfile;
      curfile = m_folder + "/" + m_curs_id->d_name;
      // Discard if it is a folder
      if (m_curs_id->d_type == DT_DIR) continue;  // FIXME: This may not be available on non-BSD systems.
 
      // Determine from the header if it is a matrix or a right hand side
      bool isvector, iscomplex = false;
      if (!getMarketHeader(curfile, m_sym, iscomplex, isvector)) continue;
      if (isvector) continue;
      if (!iscomplex) {
        EIGEN_IF_CONSTEXPR ((std::is_same<Scalar, std::complex<float> >::value ||
                             std::is_same<Scalar, std::complex<double> >::value))
          continue;
      }
      if (iscomplex) {
        EIGEN_IF_CONSTEXPR ((std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value)) {
          continue;
        }
      }
 
      // Get the matrix name
      std::string filename = m_curs_id->d_name;
      m_matname = filename.substr(0, filename.length() - 4);
 
      // Find if the matrix is SPD
      size_t found = m_matname.find("SPD");
      if ((found != std::string::npos) && (m_sym != NonSymmetric)) m_sym = SPD;
 
      m_isvalid = true;
      break;
    }
  }
  int m_sym;              // Symmetry of the matrix
  MatrixType m_mat;       // Current matrix
  VectorType m_rhs;       // Current vector
  VectorType m_refX;      // The reference solution, if exists
  std::string m_matname;  // Matrix Name
  bool m_isvalid;
  bool m_matIsLoaded;  // Determine if the matrix has already been loaded from the file
  bool m_hasRhs;       // The right hand side exists
  bool m_hasrefX;      // A reference solution is provided
  std::string m_folder;
  DIR* m_folder_id;
  struct dirent* m_curs_id;
};
 
}  // end namespace Eigen
 
#endif