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 |
|
---|
14 | namespace Eigen {
|
---|
15 |
|
---|
16 | enum {
|
---|
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 | */
|
---|
41 | template <typename Scalar>
|
---|
42 | class MatrixMarketIterator
|
---|
43 | {
|
---|
44 | public:
|
---|
45 | typedef Matrix<Scalar,Dynamic,1> VectorType;
|
---|
46 | typedef SparseMatrix<Scalar,ColMajor> MatrixType;
|
---|
47 |
|
---|
48 | public:
|
---|
49 | MatrixMarketIterator(const std::string folder):m_sym(0),m_isvalid(false),m_matIsLoaded(false),m_hasRhs(false),m_hasrefX(false),m_folder(folder)
|
---|
50 | {
|
---|
51 | m_folder_id = opendir(folder.c_str());
|
---|
52 | if (!m_folder_id){
|
---|
53 | m_isvalid = false;
|
---|
54 | std::cerr << "The provided Matrix folder could not be opened \n\n";
|
---|
55 | abort();
|
---|
56 | }
|
---|
57 | Getnextvalidmatrix();
|
---|
58 | }
|
---|
59 |
|
---|
60 | ~MatrixMarketIterator()
|
---|
61 | {
|
---|
62 | if (m_folder_id) closedir(m_folder_id);
|
---|
63 | }
|
---|
64 |
|
---|
65 | inline MatrixMarketIterator& operator++()
|
---|
66 | {
|
---|
67 | m_matIsLoaded = false;
|
---|
68 | m_hasrefX = false;
|
---|
69 | m_hasRhs = false;
|
---|
70 | Getnextvalidmatrix();
|
---|
71 | return *this;
|
---|
72 | }
|
---|
73 | inline operator bool() const { return m_isvalid;}
|
---|
74 |
|
---|
75 | /** Return the sparse matrix corresponding to the current file */
|
---|
76 | inline MatrixType& matrix()
|
---|
77 | {
|
---|
78 | // Read the matrix
|
---|
79 | if (m_matIsLoaded) return m_mat;
|
---|
80 |
|
---|
81 | std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
|
---|
82 | if ( !loadMarket(m_mat, matrix_file))
|
---|
83 | {
|
---|
84 | m_matIsLoaded = false;
|
---|
85 | return m_mat;
|
---|
86 | }
|
---|
87 | m_matIsLoaded = true;
|
---|
88 |
|
---|
89 | if (m_sym != NonSymmetric)
|
---|
90 | { // Store the upper part of the matrix. It is needed by the solvers dealing with nonsymmetric matrices ??
|
---|
91 | MatrixType B;
|
---|
92 | B = m_mat;
|
---|
93 | m_mat = B.template selfadjointView<Lower>();
|
---|
94 | }
|
---|
95 | return m_mat;
|
---|
96 | }
|
---|
97 |
|
---|
98 | /** Return the right hand side corresponding to the current matrix.
|
---|
99 | * If the rhs file is not provided, a random rhs is generated
|
---|
100 | */
|
---|
101 | inline VectorType& rhs()
|
---|
102 | {
|
---|
103 | // Get the right hand side
|
---|
104 | if (m_hasRhs) return m_rhs;
|
---|
105 |
|
---|
106 | std::string rhs_file;
|
---|
107 | rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx
|
---|
108 | m_hasRhs = Fileexists(rhs_file);
|
---|
109 | if (m_hasRhs)
|
---|
110 | {
|
---|
111 | m_rhs.resize(m_mat.cols());
|
---|
112 | m_hasRhs = loadMarketVector(m_rhs, rhs_file);
|
---|
113 | }
|
---|
114 | if (!m_hasRhs)
|
---|
115 | {
|
---|
116 | // Generate a random right hand side
|
---|
117 | if (!m_matIsLoaded) this->matrix();
|
---|
118 | m_refX.resize(m_mat.cols());
|
---|
119 | m_refX.setRandom();
|
---|
120 | m_rhs = m_mat * m_refX;
|
---|
121 | m_hasrefX = true;
|
---|
122 | m_hasRhs = true;
|
---|
123 | }
|
---|
124 | return m_rhs;
|
---|
125 | }
|
---|
126 |
|
---|
127 | /** Return a reference solution
|
---|
128 | * If it is not provided and if the right hand side is not available
|
---|
129 | * then refX is randomly generated such that A*refX = b
|
---|
130 | * where A and b are the matrix and the rhs.
|
---|
131 | * Note that when a rhs is provided, refX is not available
|
---|
132 | */
|
---|
133 | inline VectorType& refX()
|
---|
134 | {
|
---|
135 | // Check if a reference solution is provided
|
---|
136 | if (m_hasrefX) return m_refX;
|
---|
137 |
|
---|
138 | std::string lhs_file;
|
---|
139 | lhs_file = m_folder + "/" + m_matname + "_x.mtx";
|
---|
140 | m_hasrefX = Fileexists(lhs_file);
|
---|
141 | if (m_hasrefX)
|
---|
142 | {
|
---|
143 | m_refX.resize(m_mat.cols());
|
---|
144 | m_hasrefX = loadMarketVector(m_refX, lhs_file);
|
---|
145 | }
|
---|
146 | return m_refX;
|
---|
147 | }
|
---|
148 |
|
---|
149 | inline std::string& matname() { return m_matname; }
|
---|
150 |
|
---|
151 | inline int sym() { return m_sym; }
|
---|
152 |
|
---|
153 | inline bool hasRhs() {return m_hasRhs; }
|
---|
154 | inline bool hasrefX() {return m_hasrefX; }
|
---|
155 |
|
---|
156 | protected:
|
---|
157 |
|
---|
158 | inline bool Fileexists(std::string file)
|
---|
159 | {
|
---|
160 | std::ifstream file_id(file.c_str());
|
---|
161 | if (!file_id.good() )
|
---|
162 | {
|
---|
163 | return false;
|
---|
164 | }
|
---|
165 | else
|
---|
166 | {
|
---|
167 | file_id.close();
|
---|
168 | return true;
|
---|
169 | }
|
---|
170 | }
|
---|
171 |
|
---|
172 | void Getnextvalidmatrix( )
|
---|
173 | {
|
---|
174 | m_isvalid = false;
|
---|
175 | // Here, we return with the next valid matrix in the folder
|
---|
176 | while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
|
---|
177 | m_isvalid = false;
|
---|
178 | std::string curfile;
|
---|
179 | curfile = m_folder + "/" + m_curs_id->d_name;
|
---|
180 | // Discard if it is a folder
|
---|
181 | if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems
|
---|
182 | // struct stat st_buf;
|
---|
183 | // stat (curfile.c_str(), &st_buf);
|
---|
184 | // if (S_ISDIR(st_buf.st_mode)) continue;
|
---|
185 |
|
---|
186 | // Determine from the header if it is a matrix or a right hand side
|
---|
187 | bool isvector,iscomplex=false;
|
---|
188 | if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
|
---|
189 | if(isvector) continue;
|
---|
190 | if (!iscomplex)
|
---|
191 | {
|
---|
192 | if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
|
---|
193 | continue;
|
---|
194 | }
|
---|
195 | if (iscomplex)
|
---|
196 | {
|
---|
197 | if(internal::is_same<Scalar, float>::value || internal::is_same<Scalar, double>::value)
|
---|
198 | continue;
|
---|
199 | }
|
---|
200 |
|
---|
201 |
|
---|
202 | // Get the matrix name
|
---|
203 | std::string filename = m_curs_id->d_name;
|
---|
204 | m_matname = filename.substr(0, filename.length()-4);
|
---|
205 |
|
---|
206 | // Find if the matrix is SPD
|
---|
207 | size_t found = m_matname.find("SPD");
|
---|
208 | if( (found!=std::string::npos) && (m_sym != NonSymmetric) )
|
---|
209 | m_sym = SPD;
|
---|
210 |
|
---|
211 | m_isvalid = true;
|
---|
212 | break;
|
---|
213 | }
|
---|
214 | }
|
---|
215 | int m_sym; // Symmetry of the matrix
|
---|
216 | MatrixType m_mat; // Current matrix
|
---|
217 | VectorType m_rhs; // Current vector
|
---|
218 | VectorType m_refX; // The reference solution, if exists
|
---|
219 | std::string m_matname; // Matrix Name
|
---|
220 | bool m_isvalid;
|
---|
221 | bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file
|
---|
222 | bool m_hasRhs; // The right hand side exists
|
---|
223 | bool m_hasrefX; // A reference solution is provided
|
---|
224 | std::string m_folder;
|
---|
225 | DIR * m_folder_id;
|
---|
226 | struct dirent *m_curs_id;
|
---|
227 |
|
---|
228 | };
|
---|
229 |
|
---|
230 | } // end namespace Eigen
|
---|
231 |
|
---|
232 | #endif
|
---|