matrix_market.hpp 12.3 KB
Newer Older
1
// Software License for MTL
2
//
3
4
5
6
7
// Copyright (c) 2007 The Trustees of Indiana University.
//               2008 Dresden University of Technology and the Trustees of Indiana University.
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
8
//
9
// This file is part of the Matrix Template Library
10
//
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
// See also license.mtl.txt in the distribution.

#ifndef MTL_IO_MATRIX_MARKET_INCLUDE
#define MTL_IO_MATRIX_MARKET_INCLUDE

#include <string>
#include <fstream>
#include <iostream>
#include <limits>
#include <locale>
#include <complex>

#include <boost/utility/enable_if.hpp>
#include <boost/type_traits/is_floating_point.hpp>
#include <boost/type_traits/is_integral.hpp>
// #include <boost/algorithm/string/case_conv.hpp>
#include <boost/numeric/conversion/cast.hpp>

#include <boost/numeric/mtl/io/matrix_file.hpp>
#include <boost/numeric/mtl/io/read_filter.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_generator.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/utility/category.hpp>
#include <boost/numeric/mtl/utility/string_to_enum.hpp>
#include <boost/numeric/mtl/matrix/inserter.hpp>
#include <boost/numeric/mtl/operation/set_to_zero.hpp>

namespace mtl { namespace io {


/// Input file stream for files in matrix market format
class matrix_market_istream
{
    class pattern_type {};
    typedef matrix_market_istream        self;
    void check_stream(const std::string& file_name= std::string()) // not const to delete new_stream
    {
	if (!my_stream.good()) {
	    std::string message("matrix_market_istream: Error in input stream!\n");
	    if (file_name != std::string())
		message+= "Probably file " + file_name + " not found.\n";
	    std::cerr << message;
	    if (new_stream)
		delete new_stream, new_stream= 0; // To get valgrind quite
	    throw(file_not_found(message.c_str()));
	}
    }

  public:
    explicit matrix_market_istream(const char* p) : new_stream(new std::ifstream(p)), my_stream(*new_stream) { check_stream(p); }
    explicit matrix_market_istream(const std::string& s) : new_stream(new std::ifstream(s.c_str())), my_stream(*new_stream) { check_stream(s); }
    explicit matrix_market_istream(std::istream& s= std::cin) : new_stream(0), my_stream(s) { check_stream(); }

66
    ~matrix_market_istream()
67
68
69
    { 	if (new_stream) delete new_stream;    }

    template <typename Coll>
70
    self& operator>>(Coll& c)
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
    { 	return read(c, typename mtl::traits::category<Coll>::type());    }

    /// Close only my own file, i.e. if filename and not stream is passed in constructor
    void close() { if (new_stream) new_stream->close(); }

  protected:
    template <typename Matrix> self& read(Matrix& A, tag::matrix);

    void to_lower(std::string& s) const
    {
	using namespace std;
	for (unsigned i= 0; i < s.size(); i++)
	    s[i]= tolower(s[i]);
    }

    void set_symmetry(std::string& symmetry_text)
    {
88
	to_lower(symmetry_text);
89
90
91
92
93
94
	const char* symmetry_options[]= {"general", "symmetric", "skew-symmetric", "hermitian"};
	my_symmetry= string_to_enum(symmetry_text, symmetry_options, symmetry());
    }

    void set_sparsity(std::string& sparsity_text)
    {
95
	to_lower(sparsity_text);
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
	const char* sparsity_options[]= {"coordinate", "array"};
	my_sparsity= string_to_enum(sparsity_text, sparsity_options, sparsity());
    }

    template <typename Inserter, typename Value>
    void read_matrix(Inserter& ins, Value)
    {
	typedef typename Collection<typename Inserter::matrix_type>::size_type size_type;
	read_filter<Inserter> filter(ins);
	if (my_sparsity == coordinate) // sparse
	    // while (my_stream) { // sometimes does an extra erroneous loop
	    for (std::size_t i= 0; i < nnz; i++) {
		size_type r, c;
		my_stream >> r >> c;
		if (!my_stream) break; // in case while(my_stream) caught an empty line at the end
		insert_value(ins, r-1, c-1, filter, Value());
	    }
113
114
	else // dense
	    for (std::size_t c= 0; c < ncols; c++)
115
116
117
118
119
		for (std::size_t r= 0; r < nrows; r++)
		    insert_value(ins, r, c, filter, Value());
    }

    template <typename Inserter, typename Filter, typename Value>
120
    void insert_value(Inserter& ins, std::size_t r, std::size_t c, const Filter& filter, Value)
121
122
123
124
125
126
127
128
    {
	using mtl::conj;
	typedef typename Collection<typename Inserter::matrix_type>::value_type mvt;
	Value v;
	read_value(v);
	// std::cout << "Going to insert at [" << r << "][" << c << "] value " << which_value(v, mvt()) << "\n";
	if (filter(r, c))
	    ins[r][c] << which_value(v, mvt());
129
	if (r != c && filter(c, r))
130
131
132
133
134
135
136
137
138
139
140
	    switch(my_symmetry) {
	      case symmetric:      ins[c][r] << which_value(v, mvt()); break;
	      case skew:           ins[c][r] << -which_value(v, mvt()); break;
	      case Hermitian:      ins[c][r] << conj(which_value(v, mvt())); break;
	      default:             ; // do nothing
	    }
    }

    void read_value(pattern_type) {}
    void read_value(double& v) { my_stream >> v;}
    void read_value(long& v) { my_stream >> v;}
141
142
    void read_value(std::complex<double>& v)
    {
143
144
145
	double r, i; my_stream >> r >> i; v= std::complex<double>(r, i);
    }

146
    // Which value to be inserted? Itself if exist and 0 for pattern; complex are
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
    template <typename Value, typename MValue> MValue which_value(Value v, MValue) { return boost::numeric_cast<MValue>(v); }
    template <typename MValue> MValue which_value(pattern_type, MValue) { return boost::numeric_cast<MValue>(0.0); }
    template <typename MValue> MValue which_value(std::complex<double>, MValue) { MTL_THROW(runtime_error("Cannot convert complex value in real\n")); return 1; }
    std::complex<long double> which_value(std::complex<double> v, std::complex<long double>) { return boost::numeric_cast<std::complex<long double> >(v); }
    std::complex<double> which_value(std::complex<double> v, std::complex<double>) { return v; }
    std::complex<float> which_value(std::complex<double> v, std::complex<float>) { return std::complex<float>(float(real(v)), float(imag(v))); }

    std::ifstream      *new_stream;
    std::istream       &my_stream;
    enum symmetry {general, symmetric, skew, Hermitian} my_symmetry;
    enum sparsity {coordinate, array} my_sparsity;
    std::size_t nrows, ncols, nnz;
};



// Matrix version
template <typename Matrix>
matrix_market_istream& matrix_market_istream::read(Matrix& A, tag::matrix)
{
    std::string marker, type, sparsity_text, value_format, symmetry_text;
    my_stream >> marker >> type >> sparsity_text >> value_format >> symmetry_text;
169
170
#if 0
    std::cout << marker << ", " << type << ", " << sparsity_text << ", "
171
172
	      << value_format << ", " << symmetry_text << "\n";
#endif
173
    MTL_THROW_IF(marker != std::string("%%MatrixMarket"),
174
		 runtime_error("File not in Matrix Market format"));
175
    MTL_THROW_IF(type != std::string("matrix"),
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
		 runtime_error("Try to read matrix from non-matrix file"));

    set_symmetry(symmetry_text);
    set_sparsity(sparsity_text);

    char first, comment[80];
    do {
	my_stream >> first;
	if (first == '%') // comments start with % -> ignore them
	    my_stream.getline(comment, 80, '\n'); // read rest of line
	else
	    my_stream.putback(first); // if not commment we still need it
    } while (first == '%');

    my_stream >> nrows >> ncols;
191
    // std::cout << nrows << "x" << ncols << ", " << nnz << " non-zeros\n";
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
    A.change_dim(nrows, ncols);
    set_to_zero(A);

    std::size_t slot_size;
    if (sparsity_text == std::string("coordinate")) {
	my_stream >> nnz; slot_size= std::max(std::size_t(double(nnz) / double(A.dim1()) * 1.25), std::size_t(1));
    } else
	slot_size= A.dim2(); // maximal value (if A is dense it does not matter anyway)

    // Create enough space in sparse matrices
    matrix::inserter<Matrix> ins(A, slot_size);

    if (value_format == std::string("real"))
	read_matrix(ins, double());
    else if (value_format == std::string("integer"))
	read_matrix(ins, long());
    else if (value_format == std::string("complex"))
	read_matrix(ins, std::complex<double>());
    else if (value_format == std::string("pattern"))
	read_matrix(ins, pattern_type());
    else
	MTL_THROW(runtime_error("Unknown tag for matrix value type in file"));

    return *this;
}


219
class matrix_market_ostream
220
221
222
223
224
225
226
227
228
229
{
    typedef matrix_market_ostream        self;
public:
    explicit matrix_market_ostream(const char* p) : new_stream(new std::ofstream(p)), my_stream(*new_stream) {}
    explicit matrix_market_ostream(const std::string& s) : new_stream(new std::ofstream(s.c_str())), my_stream(*new_stream) {}
    explicit matrix_market_ostream(std::ostream& s= std::cout) : new_stream(0), my_stream(s) {}

    ~matrix_market_ostream() { if (new_stream) delete new_stream; }

    template <typename Coll>
230
231
    self& operator<<(const Coll& c)
    {
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
	return write(c, typename mtl::traits::category<Coll>::type());
    }

    /// Close only my own file, i.e. if filename and not stream is passed in constructor
    void close() { if (new_stream) new_stream->close(); }

private:
    template <typename Matrix> self& write(const Matrix& A, tag::matrix)
    {
	matrix_status_line(A);
	if (sparsity(A) == std::string("coordinate "))
	    return write_sparse_matrix(A);
	else
	    return write_dense_matrix(A);
    }

    template <typename Matrix> self& write_sparse_matrix(const Matrix& A)
    {
	my_stream << num_rows(A) << " " << num_cols(A) << " " << A.nnz() << "\n";
251
252
253
254

	typename mtl::traits::row<Matrix>::type             row(A);
	typename mtl::traits::col<Matrix>::type             col(A);
	typename mtl::traits::const_value<Matrix>::type     value(A);
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
	typedef typename mtl::traits::range_generator<tag::major, Matrix>::type  cursor_type;
	typedef typename mtl::traits::range_generator<tag::nz, cursor_type>::type icursor_type;

	for (cursor_type cursor = begin<tag::major>(A), cend = end<tag::major>(A); cursor != cend; ++cursor)
	    for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
		my_stream << row(*icursor)+1 << " " << col(*icursor)+1 << " ", write_value(value(*icursor)), my_stream << "\n";
	return *this;
    }

    template <typename Matrix> self& write_dense_matrix(const Matrix& A)
    {
	my_stream << num_rows(A) << " " << num_cols(A) << "\n";
	for (std::size_t c = 0; c < num_cols(A); ++c)
	    for (std::size_t r = 0; r < num_rows(A); ++r) {
		write_value(A[r][c]), my_stream << " ";
	    my_stream << "\n";
	}
	return *this;
    }

    template <typename Value>
    typename boost::enable_if<boost::is_integral<Value> >::type
    write_value(const Value& v) { my_stream << v; }

    template <typename Value>
    typename boost::enable_if<boost::is_floating_point<Value> >::type
281
282
    write_value(const Value& v)
    {
283
	my_stream.precision(std::numeric_limits<Value>::digits10 + 1);
284
285
286
	my_stream.setf(std::ios::scientific);
	my_stream << v;
	my_stream.unsetf(std::ios::scientific);
287
288
289
    }

    template <typename Value>
290
291
    void write_value(const std::complex<Value>& v)
    {
292
	my_stream.precision(std::numeric_limits<Value>::digits10 + 1);
293
294
295
	my_stream.setf(std::ios::scientific);
	my_stream << real(v) << " " << imag(v);
	my_stream.unsetf(std::ios::scientific);
296
297
298
299
300
301
302
    }

    // Will be generalized via traits::is_symmetric and alike
    template <typename Matrix>
    std::string symmetry(const Matrix&) const { return std::string("general\n"); }

    template <typename Matrix>
303
    std::string sparsity(const Matrix&) const
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
    {
	return std::string( mtl::traits::is_sparse<Matrix>::value ? "coordinate " : "array " );
    }

    template <typename Value>
    typename boost::enable_if<boost::is_integral<Value>, std::string>::type
    value(const Value&) const { return std::string("integer "); }

    template <typename Value>
    typename boost::enable_if<boost::is_floating_point<Value>, std::string>::type
    value(const Value&) const { return std::string("real "); }

    template <typename Value>
    std::string value(const std::complex<Value>&) const { return std::string("complex "); }

    template <typename Matrix>
    void matrix_status_line(const Matrix& A) const
    {
	typedef typename Collection<Matrix>::value_type value_type;
	std::string st(std::string("%%MatrixMarket  matrix ") + sparsity(A) + value(value_type()) + symmetry(A));
	my_stream << st;
    }

protected:
    std::ofstream      *new_stream;
    std::ostream       &my_stream;
};


}} // namespace mtl::io

#endif // MTL_IO_MATRIX_MARKET_INCLUDE