Line data Source code
1 : /*
2 : * Copyright (c) 2012-2013: G-CSC, Goethe University Frankfurt
3 : * Author: Torbjörn Klatt
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : /**
34 : * \file matrix_io_mtx.h
35 : * \author Torbjoern Klatt
36 : * \date 2012-05-06
37 : *
38 : * \details For some parts the implementation of a pure C library for the
39 : * MatrixMarket exchange format was used as a source of inspiration.
40 : * These can be found at http://math.nist.gov/MatrixMarket/mmio-c.html
41 : */
42 :
43 : #ifndef __H__UG__LIB_ALGEBRA__MATRIX_IO_MTX_H
44 : #define __H__UG__LIB_ALGEBRA__MATRIX_IO_MTX_H
45 :
46 : #include <sstream>
47 : #include <vector>
48 : #include <boost/algorithm/string.hpp>
49 : #include <boost/lexical_cast.hpp>
50 :
51 : #include "lib_algebra/common/matrixio/matrix_io.h"
52 : #include "lib_algebra/common/matrixio/mm_type_code.h"
53 :
54 : namespace ug
55 : {
56 :
57 : /// \addtogroup matrixio
58 : /// \{
59 :
60 : /**
61 : * \brief Provides I/O functionality for MatrixMarket exchange file format
62 : *
63 : * Given the path to a \c .mtx file it queries the therein stored matrix type
64 : * and characteristics.
65 : * With <tt>readInto(matrixType &matrix)</tt> an easy way of reading the
66 : * represented matrix into a <tt>CPUAlgebra::matrix_type</tt> object is
67 : * provided.
68 : *
69 : * \note So far, only real sparse matrices either symmetric or skew-symmetric or
70 : * none of both are supported.
71 : * Support for dense matrices (either really dense or just in array format) will
72 : * come soon.
73 : * Non-real (i.e. complex, integer or pattern) matrices and hermitian matrices
74 : * might be supported in the future.
75 : *
76 : * The specification of the MatrixMarket matrix exchange format can be found
77 : * here: http://math.nist.gov/MatrixMarket/formats.html#MMformat
78 : *
79 : * \note Please include the generic MatrixIO header file
80 : * (lib_algebra/common/matrixio/matrix_io.h) for usage and not the header file
81 : * of this specific exchange file format.
82 : */
83 : class MatrixIOMtx : private MatrixIO
84 : {
85 : private:
86 : // Full path name of the matrix exchange file
87 : // (docu in matrix_io.h)
88 : std::string *m_pMatFileName;
89 : // Internal file stream for reading from and writing into the matrix
90 : // exchange file (docu in matrix_io.h)
91 : std::fstream m_matFileStream;
92 : /// Line number of the first data line (0 based)
93 : size_t m_firstDataLine;
94 : // Matrix exchange file type (set to constant MatrixFileType::MatrixMarket)
95 : // (docu in matrix_io.h)
96 : MatrixFileType m_matFileType;
97 : /// Characteristics of MatrixMarket matrix
98 : MMTypeCode m_mmTypeCode;
99 : // Number of rows as specified in the matrix exchange file
100 : // (docu in matrix_io.h)
101 : size_t m_rows;
102 : // Number of columns as specified in the matrix exchange file
103 : // (docu in matrix_io.h)
104 : size_t m_cols;
105 : // Number of data lines as specified in the matrix exchange file
106 : // (docu in matrix_io.h)
107 : size_t m_lines;
108 :
109 : public:
110 : /**
111 : * \brief Default constructor
112 : */
113 : MatrixIOMtx();
114 : /**
115 : * \brief Constructor specifying matrix exchange file path
116 : *
117 : * After successful call to MatrixIOMtx::set_mat_file_name it calls
118 : * MatrixIOMtx::query_matrix_type and
119 : * MatrixIOMtx::query_matrix_characteristics as long as openMode is
120 : * MatrixIO::OpenMode::EXISTING.
121 : *
122 : * \param[in] mFile Full path and name of matrix exchange file.
123 : * \param[in] openMode how to deal with non-existing files (a value of
124 : * MatrixIO::OpenMode)
125 : */
126 : MatrixIOMtx( std::string mFile, int openMode=MatrixIO::EXISTING );
127 :
128 : /**
129 : * \brief Destructor
130 : *
131 : * The destructor calls MatrixIOMtx::close_file to make sure, that write
132 : * operations to the associated exchange file are flushed and terminated
133 : * correctly.
134 : */
135 : ~MatrixIOMtx();
136 :
137 : /**
138 : * \brief Sets associated file to another one
139 : *
140 : * \note It does not check, whether the given file is a MatrixMarket
141 : * exchange file.
142 : * Only it's existance is verified.
143 : *
144 : * \param[in] mFile Full path and name of matrix exchange file.
145 : * \param[in] openMode how to deal with non-existing files
146 : *
147 : * \throws std::runtime_error if exchange file is not found and
148 : * <tt>openMode=MatrixIO::OpenMode::Existing</tt>.
149 : * \throws UGError if \c openMode is non of MatrixIO::OpenMode
150 : */
151 : void set_mat_file_name( std::string mFile, int openMode=MatrixIO::EXISTING );
152 :
153 : /**
154 : * \brief Retreive associated exchange file path and name
155 : *
156 : * \return String representation of exchange file with full path as set by
157 : * constructor or MatrixIOMtx::set_mat_file_name.
158 : */
159 : std::string get_mat_file_name() const;
160 :
161 : /**
162 : * \brief Sets the size of the matrix and number of data lines
163 : *
164 : * In case of reading from a matrix exchange file, this function is called
165 : * by MatrixIOMtx::query_matrix_characteristics.
166 : *
167 : * \param[in] rows Number of rows of the matrix
168 : * \param[in] cols Number of columns of the matrix
169 : * \param[in] lines Number of data lines in the exchange file
170 : *
171 : * \throws std::runtime_error if \c rows is not positive
172 : * \throws std::runtime_error if \c cols is not positive
173 : * \throws std::runtime_error if \c lines is not positive
174 : */
175 : void set_mat_dims( size_t rows, size_t cols, size_t lines );
176 :
177 : /**
178 : * \brief Retrieves number of rows as specified in the exchange file
179 : *
180 : * \return Number of rows as \c size_t
181 : */
182 : size_t get_num_rows() const;
183 :
184 : /**
185 : * \brief Retrieves number of columns as specified in the exchange file
186 : *
187 : * \return Number of columns as \c size_t
188 : */
189 : size_t get_num_cols() const;
190 :
191 : /**
192 : * \brief Retrieves number of data lines as specified in the exchange file
193 : *
194 : * \return Number of data lines as \c size_t
195 : */
196 : size_t get_num_lines() const;
197 :
198 : /**
199 : * \brief Tells whether the associated exchange file holds a sparse matrix
200 : *
201 : * \return \c true if it is in coordinate format, \c false otherwise
202 : */
203 : bool is_sparse() const;
204 :
205 : /**
206 : * \brief Reads data from associated exchange file into ug matrix
207 : *
208 : * The given ug-matrix is first resized.
209 : * Each data line of the underlying exchange file is read separately and
210 : * independently.
211 : * If the matrix is (skew-)symmetric as defined in the exchange file, this
212 : * is considered during writing the values to the correct positions.
213 : *
214 : * \note This function does not take care about memory alignement of the
215 : * given ug-matrix (e.g. SparseMatrix::defragment ).
216 : *
217 : * \param[in,out] matrix Matrix of CPUAlgebra::matrix_type to read into
218 : *
219 : * \throws std::runtime_error if one of \c m_rows, \c m_cols or \c m_lines
220 : * is not positive.
221 : * \throws UGError if exchange file's matrix is not in coordinate
222 : * format (as only sparse matrices are supported yet)
223 : */
224 : template<typename matrix_type>
225 6 : void read_into( matrix_type &matrix )
226 : {
227 : PROFILE_FUNC();
228 : UG_ASSERT( m_rows > 0 || m_cols > 0 || m_lines > 0,
229 : "MatrixMarket matrix dimensions seem not be reasonable: (m, n, nnz)=("
230 : << m_rows << "," << m_cols << "," << m_lines << ")" );
231 :
232 : // open the file and go one before first data line
233 6 : open_file();
234 :
235 : std::string dummy;
236 18 : for( size_t i = 0; i < m_firstDataLine; i++ ) {
237 12 : std::getline( m_matFileStream, dummy );
238 : }
239 :
240 6 : matrix.resize_and_clear( m_rows, m_cols );
241 : size_t x, y;
242 : double val;
243 6 : if ( m_mmTypeCode.is_sparse() ) {
244 15728 : for( size_t i = 0; i < m_lines; i++ ) {
245 15722 : read_entry( &x, &y, &val );
246 : // we need to substract 1 from row and column indices as MM is 1-based
247 15722 : matrix( x - 1, y - 1 ) = val;
248 15722 : if ( m_mmTypeCode.is_symmetric() && x != y ) {
249 6844 : matrix( y - 1, x - 1 ) = val;
250 8878 : } else if ( m_mmTypeCode.is_skew_symmetric() && x != y ) {
251 1760 : matrix( y - 1, x - 1 ) = -val;
252 : }
253 : }
254 : } else {
255 0 : UG_THROW( "Other than sparse MatrixMarket matrices are not yet implemented." );
256 : }
257 :
258 : // close the file
259 6 : close_file();
260 6 : }
261 :
262 :
263 : /**
264 : * \brief Writes a ug matrix into the associated MatrixMarket exchange file
265 : *
266 : * From this analysis the MatrixMarket exchange file banner is constructed
267 : * and the characteristics and data lines written to that file consecutively.
268 : *
269 : * A comment is added right after the banner.
270 : * The default comment is <tt>Generated with ug4.</tt>.
271 : * To suppress the comment line completely, pass an empty string as the
272 : * \c comment parameter to this function.
273 : *
274 : * \note Please make sure, the lines of the given comment are not longer
275 : * than 1025 characters and all liness start with a \c %.
276 : *
277 : * \todo Implement validity check for the comment.
278 : *
279 : * The resulting MatrixMarket exchange file is in line with the
280 : * specifications.
281 : *
282 : * \param[in] matrix Matrix of CPUAlgebra::matrix_type to read from
283 : * \param[in] comment Optional comment to be put atop of the MatrixMarket
284 : * matrix exchange file.
285 : */
286 : template<typename matrix_type>
287 3 : void write_from( matrix_type &matrix,
288 : std::string comment="%Generated with ug4." )
289 : {
290 : PROFILE_FUNC();
291 :
292 : // analyze the given matrix
293 3 : std::vector< std::vector<size_t> > rowIndexPerCol = determine_matrix_characteristics( matrix );
294 :
295 : // write the MatrixMarket banner
296 3 : write_banner();
297 :
298 : // open the file for appending
299 3 : open_file( std::ios_base::out | std::ios_base::app );
300 :
301 : // add a comment if it's not empty
302 3 : if ( !comment.empty() ) {
303 0 : if ( comment.find_first_of( '%' ) != 0 ) {
304 : UG_WARNING( "Given comment did not start with '%'. Prepending it to make it valid." );
305 0 : comment.insert( 0, "%" );
306 : }
307 0 : m_matFileStream << comment << "\n";
308 : }
309 :
310 : // write characteristics
311 9 : m_matFileStream << m_rows << " " << m_cols << " " << m_lines << "\n";
312 :
313 :
314 : // write entries to the file
315 1445 : for ( size_t col = 0; col < m_cols; col++ ) {
316 17164 : for (size_t row = 0; row < rowIndexPerCol.at(col).size(); row++ ) {
317 : // we need to add 1 to row and column index, as MM indices are 1-based
318 15722 : write_entry( rowIndexPerCol.at(col).at(row) + 1, col + 1,
319 7861 : matrix(rowIndexPerCol.at(col).at(row), col) );
320 : }
321 : }
322 3 : close_file();
323 3 : }
324 :
325 :
326 : private:
327 : // Opens file as fstream. (docu in matrix_io.h)
328 : void open_file( std::ios_base::openmode mode=std::ios_base::in );
329 :
330 : // Closes fstream. (docu in matrix_io.h)
331 : void close_file();
332 :
333 : /**
334 : * \brief Parses banner of MatrixMarket file
335 : *
336 : * The banner is the first line of a MatrixMarket file with the syntax:
337 : *
338 : * <tt>%%MatrixMarket matrix [format] [numeric] [type]</tt>
339 : *
340 : * <dl>
341 : * <dt><tt>format</tt></dt>
342 : * <dd>Either \c coordinate or \c array, representing a sparse or dense
343 : * matrix respectively</dd>
344 : * <dt><tt>numeric</tt></dt>
345 : * <dd>Either \c real, \c complex, \c integer or \c pattern,
346 : * representing the numerical type of the values</dd>
347 : * <dt><tt>type</tt></dt>
348 : * <dd>Either \c general, \c symmetric, \c skew-symmetric or
349 : * \c hermitian, representing the algebraic type of the matrix</dd>
350 : * </dl>
351 : *
352 : * \throws std::runtime_error if the banner starts not with
353 : * <tt>%%MatrixMarket</tt>
354 : * \throws std::runtime_error if the second word of the banner is not
355 : * \c matrix
356 : * \throws UGError if the \c format identifier is invalid
357 : * \throws UGError if the \c numeric identifier is invalid
358 : * \throws UGError if the \c type identifier is invalid
359 : */
360 : void query_matrix_type();
361 :
362 : /**
363 : * \brief Retreives matrix size and number non-zeros from exchange file
364 : *
365 : * The first non-empty line after the comments (lines starting with \c %)
366 : * consists of three values for the coordinate format:
367 : *
368 : * <tt>nRows nCols nLines</tt>
369 : *
370 : * <dl>
371 : * <dt><tt>nRows</tt></dt><dd>Number of rows of the described matrix</dd>
372 : * <dt><tt>nCols</tt></dt><dd>Number of columns of the described matrix</dd>
373 : * <dt><tt>nLines</tt></dt><dd>Number of data lines. For types \c symmetric,
374 : * \c skew-symmetric and \c hermitian, this is not equal to the number of
375 : * non-zero elements.</dd>
376 : * </dl>
377 : *
378 : * \throws std::runtime_error if the file stream ended unexpectetly (e.g.
379 : * the exchange file does not contain data lines)
380 : * \throws UGError if the exchange file does not describe a
381 : * coordinate/sparse matrix (as only sparse
382 : * matrices are supported yet)
383 : */
384 : void query_matrix_characteristics();
385 :
386 : /**
387 : * \brief Determines characteristics of given matrix
388 : *
389 : * Symmetries of the given matrix are determined and saved.
390 : *
391 : * \param[in] matrix Matrix of CPUAlgebra::matrix_type to analyse
392 : * \return 2D vector with coordinates of non-zero entries of the given matrix
393 : * with column indices as first and row indices as second dimension.
394 : *
395 : * \throws std::runtime_error if the matrix was determined to be symmetric
396 : * and skew-symmetric at the same time
397 : */
398 :
399 :
400 : template<typename matrix_type>
401 3 : std::vector< std::vector<size_t> > determine_matrix_characteristics(const matrix_type &matrix )
402 : {
403 : PROFILE_FUNC();
404 :
405 : // read matrix dimensions
406 : size_t rows = matrix.num_rows();
407 : size_t cols = matrix.num_cols();
408 : size_t offDiagEntries = 0;
409 : size_t diagEntries = 0;
410 :
411 : // As MatrixMarket specifies a column-first-ordering, we need to get to know
412 : // in which rows of each column non-zero entries are.
413 : // During this, we also get to know how many non-zero entries there are in
414 : // total and can detect symmetries of the matrix.
415 : std::vector< std::vector<size_t> > rowIndexPerCol;
416 3 : rowIndexPerCol.resize(cols);
417 : bool isSymmetric = true;
418 : bool isSkew = true;
419 : bool changed = false;
420 1448 : for ( size_t r = 0; r < rows; r++ ) { // iterate rows
421 : for (typename matrix_type::const_row_iterator conn = matrix.begin_row(r);
422 13610 : conn != matrix.end_row(r); ++conn ) {
423 12168 : if ( conn.value() != 0.0 ) {
424 : // first add index to list
425 12168 : if ( isSymmetric || isSkew ) {
426 9509 : if ( conn.index() <= r ) {
427 5204 : rowIndexPerCol.at( conn.index() ).push_back(r);
428 : }
429 : } else {
430 2659 : rowIndexPerCol.at( conn.index() ).push_back(r);
431 : }
432 :
433 : // increment counters
434 12168 : (conn.index() == r) ? diagEntries++ : offDiagEntries++ ;
435 :
436 : // check, whether it's still a symmetric or skew-symmetric matrix
437 12168 : if ( r != conn.index() ) {
438 21291 : if ( matrix(r, conn.index() ) != matrix( conn.index(), r ) ) {
439 4241 : if ( isSymmetric ) {
440 : changed = true;
441 : }
442 : isSymmetric = false;
443 : }
444 11086 : if ( matrix(r, conn.index() ) != -1.0 * matrix( conn.index(), r ) ) {
445 9325 : if ( isSkew ) {
446 : changed = true;
447 : }
448 : isSkew = false;
449 : }
450 : }
451 :
452 : // We assumed the matrix to be symmetric or skew-symmetric, but it's
453 : // not. Thus we need to redo everything done before.
454 12168 : if ( changed ) {
455 : offDiagEntries = 0;
456 : diagEntries = 0;
457 : rowIndexPerCol.clear();
458 3 : rowIndexPerCol.resize( cols );
459 3 : r = -1; // at the end of the current row-loop this is incremented again
460 : changed = false;
461 3 : break;
462 : }
463 : }
464 : }
465 : }
466 :
467 : // make sure the matrix is not both, symmetric and skew-symmetric
468 : UG_ASSERT( ( ( isSymmetric && !isSkew ) ||
469 : ( isSkew && !isSymmetric ) ||
470 : ( !isSkew && !isSymmetric ) ),
471 : "Error on detecting symmetry of matrix. (skew:" << isSkew
472 : << " sym:" << isSymmetric << ")" );
473 :
474 : // set MMTypeCode
475 : size_t entries = 0;
476 3 : m_mmTypeCode.set_class_type( MMTypeCode::COORDINATE );
477 3 : m_mmTypeCode.set_numeric_type( MMTypeCode::REAL );
478 3 : if ( isSymmetric ) {
479 1 : entries = ( offDiagEntries / 2 )+ diagEntries;
480 1 : m_mmTypeCode.set_algebraic_type( MMTypeCode::SYMMETRIC );
481 2 : } else if ( isSkew ) {
482 1 : entries = ( offDiagEntries / 2 ) + diagEntries;
483 1 : m_mmTypeCode.set_algebraic_type( MMTypeCode::SKEW );
484 : } else {
485 1 : entries = offDiagEntries + diagEntries;
486 1 : m_mmTypeCode.set_algebraic_type( MMTypeCode::GENERAL );
487 : }
488 :
489 : // now we can set the dimensions
490 3 : set_mat_dims( rows, cols, entries );
491 :
492 3 : return rowIndexPerCol;
493 0 : }
494 :
495 : /**
496 : * \brief Reads and parses the next data line
497 : *
498 : * \note This function does not open the filestream by itself, but requires
499 : * it to be already open.
500 : *
501 : * \param[out] m pointer to store the row index to (1 based)
502 : * \param[out] n pointer to store the column index to (1 based)
503 : * \param[out] val pointer to store the value to
504 : *
505 : * \throws std::runtime_error if the filestream is not open
506 : * \throws std::runtime_error if the matrix is sparse but there were less or
507 : * more than 3 elements in the line
508 : * \throws ug::UGError if the exchange file does not describe a
509 : * coordinate/sparse matrix (as only sparse
510 : * matrices are supported yet)
511 : */
512 : void read_entry( size_t *m, size_t *n, CPUAlgebra::matrix_type::value_type *val );
513 :
514 : /**
515 : * \brief Writes the banner to the associated MatrixMarket file
516 : *
517 : * \note Any existing content of the associated file will be delted.
518 : */
519 : void write_banner();
520 :
521 : /**
522 : * \brief Appends a data line to the associated file
523 : *
524 : * The value is formatted to scientific notation with 13 digits of precision.
525 : *
526 : * \param[out] m row index (1 based)
527 : * \param[out] n column index (1 based)
528 : * \param[out] val value at [row,col]
529 : *
530 : * \throws std::runtime_error if the filestream is not open
531 : * \throws std::runtime_error if either the specified row or column index
532 : * is not positive
533 : */
534 : void write_entry( size_t m, size_t n, CPUAlgebra::matrix_type::value_type val );
535 : };
536 :
537 : // end group matrixio
538 : /// \}
539 :
540 : } // namespace ug
541 :
542 : #endif // __H__UG__LIB_ALGEBRA__MATRIX_IO_MTX_H
543 :
544 : // EOF
|
