1
2// This file is part of Eigen, a lightweight C++ template library
3// for linear algebra.
4//
5// Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#ifndef EIGEN_BROWSE_MATRICES_H
12#define EIGEN_BROWSE_MATRICES_H
13
14namespace Eigen {
15
16enum {
17  SPD = 0x100,
18  NonSymmetric = 0x0
19};
20
21/**
22 * @brief Iterator to browse matrices from a specified folder
23 *
24 * This is used to load all the matrices from a folder.
25 * The matrices should be in Matrix Market format
26 * It is assumed that the matrices are named as matname.mtx
27 * and matname_SPD.mtx if the matrix is Symmetric and positive definite (or Hermitian)
28 * The right hand side vectors are loaded as well, if they exist.
29 * They should be named as matname_b.mtx.
30 * Note that the right hand side for a SPD matrix is named as matname_SPD_b.mtx
31 *
32 * Sometimes a reference solution is available. In this case, it should be named as matname_x.mtx
33 *
34 * Sample code
35 * \code
36 *
37 * \endcode
38 *
39 * \tparam Scalar The scalar type
40 */
41template <typename Scalar>
42class MatrixMarketIterator
43{
44    typedef typename NumTraits<Scalar>::Real RealScalar;
45  public:
46    typedef Matrix<Scalar,Dynamic,1> VectorType;
47    typedef SparseMatrix<Scalar,ColMajor> MatrixType;
48
49  public:
50    MatrixMarketIterator(const std::string &folder)
51      : m_sym(0), m_isvalid(false), m_matIsLoaded(false), m_hasRhs(false), m_hasrefX(false), m_folder(folder)
52    {
53      m_folder_id = opendir(folder.c_str());
54      if(m_folder_id)
55        Getnextvalidmatrix();
56    }
57
58    ~MatrixMarketIterator()
59    {
60      if (m_folder_id) closedir(m_folder_id);
61    }
62
63    inline MatrixMarketIterator& operator++()
64    {
65      m_matIsLoaded = false;
66      m_hasrefX = false;
67      m_hasRhs = false;
68      Getnextvalidmatrix();
69      return *this;
70    }
71    inline operator bool() const { return m_isvalid;}
72
73    /** Return the sparse matrix corresponding to the current file */
74    inline MatrixType& matrix()
75    {
76      // Read the matrix
77      if (m_matIsLoaded) return m_mat;
78
79      std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
80      if ( !loadMarket(m_mat, matrix_file))
81      {
82        std::cerr << "Warning loadMarket failed when loading \"" << matrix_file << "\"" << std::endl;
83        m_matIsLoaded = false;
84        return m_mat;
85      }
86      m_matIsLoaded = true;
87
88      if (m_sym != NonSymmetric)
89      {
90        // Check whether we need to restore a full matrix:
91        RealScalar diag_norm  = m_mat.diagonal().norm();
92        RealScalar lower_norm = m_mat.template triangularView<Lower>().norm();
93        RealScalar upper_norm = m_mat.template triangularView<Upper>().norm();
94        if(lower_norm>diag_norm && upper_norm==diag_norm)
95        {
96          // only the lower part is stored
97          MatrixType tmp(m_mat);
98          m_mat = tmp.template selfadjointView<Lower>();
99        }
100        else if(upper_norm>diag_norm && lower_norm==diag_norm)
101        {
102          // only the upper part is stored
103          MatrixType tmp(m_mat);
104          m_mat = tmp.template selfadjointView<Upper>();
105        }
106      }
107      return m_mat;
108    }
109
110    /** Return the right hand side corresponding to the current matrix.
111     * If the rhs file is not provided, a random rhs is generated
112     */
113    inline VectorType& rhs()
114    {
115       // Get the right hand side
116      if (m_hasRhs) return m_rhs;
117
118      std::string rhs_file;
119      rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx
120      m_hasRhs = Fileexists(rhs_file);
121      if (m_hasRhs)
122      {
123        m_rhs.resize(m_mat.cols());
124        m_hasRhs = loadMarketVector(m_rhs, rhs_file);
125      }
126      if (!m_hasRhs)
127      {
128        // Generate a random right hand side
129        if (!m_matIsLoaded) this->matrix();
130        m_refX.resize(m_mat.cols());
131        m_refX.setRandom();
132        m_rhs = m_mat * m_refX;
133        m_hasrefX = true;
134        m_hasRhs = true;
135      }
136      return m_rhs;
137    }
138
139    /** Return a reference solution
140     * If it is not provided and if the right hand side is not available
141     * then refX is randomly generated such that A*refX = b
142     * where A and b are the matrix and the rhs.
143     * Note that when a rhs is provided, refX is not available
144     */
145    inline VectorType& refX()
146    {
147      // Check if a reference solution is provided
148      if (m_hasrefX) return m_refX;
149
150      std::string lhs_file;
151      lhs_file = m_folder + "/" + m_matname + "_x.mtx";
152      m_hasrefX = Fileexists(lhs_file);
153      if (m_hasrefX)
154      {
155        m_refX.resize(m_mat.cols());
156        m_hasrefX = loadMarketVector(m_refX, lhs_file);
157      }
158      else
159        m_refX.resize(0);
160      return m_refX;
161    }
162
163    inline std::string& matname() { return m_matname; }
164
165    inline int sym() { return m_sym; }
166
167    bool hasRhs() {return m_hasRhs; }
168    bool hasrefX() {return m_hasrefX; }
169    bool isFolderValid() { return bool(m_folder_id); }
170
171  protected:
172
173    inline bool Fileexists(std::string file)
174    {
175      std::ifstream file_id(file.c_str());
176      if (!file_id.good() )
177      {
178        return false;
179      }
180      else
181      {
182        file_id.close();
183        return true;
184      }
185    }
186
187    void Getnextvalidmatrix( )
188    {
189      m_isvalid = false;
190      // Here, we return with the next valid matrix in the folder
191      while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
192        m_isvalid = false;
193        std::string curfile;
194        curfile = m_folder + "/" + m_curs_id->d_name;
195        // Discard if it is a folder
196        if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems
197//         struct stat st_buf;
198//         stat (curfile.c_str(), &st_buf);
199//         if (S_ISDIR(st_buf.st_mode)) continue;
200
201        // Determine from the header if it is a matrix or a right hand side
202        bool isvector,iscomplex=false;
203        if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
204        if(isvector) continue;
205        if (!iscomplex)
206        {
207          if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
208            continue;
209        }
210        if (iscomplex)
211        {
212          if(internal::is_same<Scalar, float>::value || internal::is_same<Scalar, double>::value)
213            continue;
214        }
215
216
217        // Get the matrix name
218        std::string filename = m_curs_id->d_name;
219        m_matname = filename.substr(0, filename.length()-4);
220
221        // Find if the matrix is SPD
222        size_t found = m_matname.find("SPD");
223        if( (found!=std::string::npos) && (m_sym != NonSymmetric) )
224          m_sym = SPD;
225
226        m_isvalid = true;
227        break;
228      }
229    }
230    int m_sym; // Symmetry of the matrix
231    MatrixType m_mat; // Current matrix
232    VectorType m_rhs;  // Current vector
233    VectorType m_refX; // The reference solution, if exists
234    std::string m_matname; // Matrix Name
235    bool m_isvalid;
236    bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file
237    bool m_hasRhs; // The right hand side exists
238    bool m_hasrefX; // A reference solution is provided
239    std::string m_folder;
240    DIR * m_folder_id;
241    struct dirent *m_curs_id;
242
243};
244
245} // end namespace Eigen
246
247#endif
248