ug4
Loading...
Searching...
No Matches
matrix_io_mtx.h
Go to the documentation of this file.
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
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
53
54namespace ug
55{
56
59
83class 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;
94 // Matrix exchange file type (set to constant MatrixFileType::MatrixMarket)
95 // (docu in matrix_io.h)
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:
113 MatrixIOMtx();
126 MatrixIOMtx( std::string mFile, int openMode=MatrixIO::EXISTING );
127
135 ~MatrixIOMtx();
136
151 void set_mat_file_name( std::string mFile, int openMode=MatrixIO::EXISTING );
152
159 std::string get_mat_file_name() const;
160
175 void set_mat_dims( size_t rows, size_t cols, size_t lines );
176
182 size_t get_num_rows() const;
183
189 size_t get_num_cols() const;
190
196 size_t get_num_lines() const;
197
203 bool is_sparse() const;
204
224 template<typename matrix_type>
225 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 open_file();
234
235 std::string dummy;
236 for( size_t i = 0; i < m_firstDataLine; i++ ) {
237 std::getline( m_matFileStream, dummy );
238 }
239
240 matrix.resize_and_clear( m_rows, m_cols );
241 size_t x, y;
242 double val;
243 if ( m_mmTypeCode.is_sparse() ) {
244 for( size_t i = 0; i < m_lines; i++ ) {
245 read_entry( &x, &y, &val );
246 // we need to substract 1 from row and column indices as MM is 1-based
247 matrix( x - 1, y - 1 ) = val;
248 if ( m_mmTypeCode.is_symmetric() && x != y ) {
249 matrix( y - 1, x - 1 ) = val;
250 } else if ( m_mmTypeCode.is_skew_symmetric() && x != y ) {
251 matrix( y - 1, x - 1 ) = -val;
252 }
253 }
254 } else {
255 UG_THROW( "Other than sparse MatrixMarket matrices are not yet implemented." );
256 }
257
258 // close the file
259 close_file();
260 }
261
262
286 template<typename matrix_type>
287 void write_from( matrix_type &matrix,
288 std::string comment="%Generated with ug4." )
289 {
290 PROFILE_FUNC();
291
292 // analyze the given matrix
293 std::vector< std::vector<size_t> > rowIndexPerCol = determine_matrix_characteristics( matrix );
294
295 // write the MatrixMarket banner
296 write_banner();
297
298 // open the file for appending
299 open_file( std::ios_base::out | std::ios_base::app );
300
301 // add a comment if it's not empty
302 if ( !comment.empty() ) {
303 if ( comment.find_first_of( '%' ) != 0 ) {
304 UG_WARNING( "Given comment did not start with '%'. Prepending it to make it valid." );
305 comment.insert( 0, "%" );
306 }
307 m_matFileStream << comment << "\n";
308 }
309
310 // write characteristics
311 m_matFileStream << m_rows << " " << m_cols << " " << m_lines << "\n";
312
313
314 // write entries to the file
315 for ( size_t col = 0; col < m_cols; col++ ) {
316 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 write_entry( rowIndexPerCol.at(col).at(row) + 1, col + 1,
319 matrix(rowIndexPerCol.at(col).at(row), col) );
320 }
321 }
322 close_file();
323 }
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
360 void query_matrix_type();
361
385
400 template<typename matrix_type>
401 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 rowIndexPerCol.resize(cols);
417 bool isSymmetric = true;
418 bool isSkew = true;
419 bool changed = false;
420 for ( size_t r = 0; r < rows; r++ ) { // iterate rows
421 for (typename matrix_type::const_row_iterator conn = matrix.begin_row(r);
422 conn != matrix.end_row(r); ++conn ) {
423 if ( conn.value() != 0.0 ) {
424 // first add index to list
425 if ( isSymmetric || isSkew ) {
426 if ( conn.index() <= r ) {
427 rowIndexPerCol.at( conn.index() ).push_back(r);
428 }
429 } else {
430 rowIndexPerCol.at( conn.index() ).push_back(r);
431 }
432
433 // increment counters
434 (conn.index() == r) ? diagEntries++ : offDiagEntries++ ;
435
436 // check, whether it's still a symmetric or skew-symmetric matrix
437 if ( r != conn.index() ) {
438 if ( matrix(r, conn.index() ) != matrix( conn.index(), r ) ) {
439 if ( isSymmetric ) {
440 changed = true;
441 }
442 isSymmetric = false;
443 }
444 if ( matrix(r, conn.index() ) != -1.0 * matrix( conn.index(), r ) ) {
445 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 if ( changed ) {
455 offDiagEntries = 0;
456 diagEntries = 0;
457 rowIndexPerCol.clear();
458 rowIndexPerCol.resize( cols );
459 r = -1; // at the end of the current row-loop this is incremented again
460 changed = false;
461 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;
478 if ( isSymmetric ) {
479 entries = ( offDiagEntries / 2 )+ diagEntries;
481 } else if ( isSkew ) {
482 entries = ( offDiagEntries / 2 ) + diagEntries;
484 } else {
485 entries = offDiagEntries + diagEntries;
487 }
488
489 // now we can set the dimensions
490 set_mat_dims( rows, cols, entries );
491
492 return rowIndexPerCol;
493 }
494
512 void read_entry( size_t *m, size_t *n, CPUAlgebra::matrix_type::value_type *val );
513
519 void write_banner();
520
534 void write_entry( size_t m, size_t n, CPUAlgebra::matrix_type::value_type val );
535};
536
537// end group matrixio
539
540} // namespace ug
541
542#endif // __H__UG__LIB_ALGEBRA__MATRIX_IO_MTX_H
543
544// EOF
Type representation for MatrixMarket matrix exchange files.
Definition mm_type_code.h:60
void set_class_type(int type)
Sets a new class type from an enum value.
Definition mm_type_code.h:202
bool is_skew_symmetric() const
Tells whether MMTypeCode is skew-symmetric.
Definition mm_type_code.h:192
void set_numeric_type(int type)
Sets a new numeric type from an enum value.
Definition mm_type_code.h:236
void set_algebraic_type(int type)
Sets a new algebraic type from an enum value.
Definition mm_type_code.h:281
bool is_symmetric() const
Tells whether MMTypeCode is symmetric.
Definition mm_type_code.h:184
bool is_sparse() const
Tells whether MMTypeCode is sparse.
Definition mm_type_code.h:158
@ REAL
Definition mm_type_code.h:78
@ GENERAL
Definition mm_type_code.h:90
@ SYMMETRIC
Definition mm_type_code.h:91
@ SKEW
Definition mm_type_code.h:92
@ COORDINATE
Definition mm_type_code.h:68
Representation of a matrix exchange file format.
Definition matrix_io.h:66
Generic matrix I/O functionality.
Definition matrix_io.h:126
@ EXISTING
Only existing files are valid files (i.e. non-existing file result in error)
Definition matrix_io.h:131
Provides I/O functionality for MatrixMarket exchange file format.
Definition matrix_io_mtx.h:84
void read_into(matrix_type &matrix)
Reads data from associated exchange file into ug matrix.
Definition matrix_io_mtx.h:225
size_t m_lines
Definition matrix_io_mtx.h:107
~MatrixIOMtx()
Destructor.
Definition matrix_io_mtx.cpp:65
size_t m_firstDataLine
Line number of the first data line (0 based)
Definition matrix_io_mtx.h:93
size_t m_rows
Definition matrix_io_mtx.h:101
void query_matrix_type()
Parses banner of MatrixMarket file.
Definition matrix_io_mtx.cpp:152
std::fstream m_matFileStream
Definition matrix_io_mtx.h:91
MatrixFileType m_matFileType
Definition matrix_io_mtx.h:96
MatrixIOMtx()
Default constructor.
Definition matrix_io_mtx.cpp:40
size_t get_num_lines() const
Retrieves number of data lines as specified in the exchange file.
Definition matrix_io_mtx.cpp:119
void open_file(std::ios_base::openmode mode=std::ios_base::in)
Definition matrix_io_mtx.cpp:132
void query_matrix_characteristics()
Retreives matrix size and number non-zeros from exchange file.
Definition matrix_io_mtx.cpp:199
void write_banner()
Writes the banner to the associated MatrixMarket file.
Definition matrix_io_mtx.cpp:267
std::vector< std::vector< size_t > > determine_matrix_characteristics(const matrix_type &matrix)
Determines characteristics of given matrix.
Definition matrix_io_mtx.h:401
MMTypeCode m_mmTypeCode
Characteristics of MatrixMarket matrix.
Definition matrix_io_mtx.h:98
void write_from(matrix_type &matrix, std::string comment="%Generated with ug4.")
Writes a ug matrix into the associated MatrixMarket exchange file.
Definition matrix_io_mtx.h:287
void set_mat_dims(size_t rows, size_t cols, size_t lines)
Sets the size of the matrix and number of data lines.
Definition matrix_io_mtx.cpp:97
void close_file()
Definition matrix_io_mtx.cpp:144
std::string get_mat_file_name() const
Retreive associated exchange file path and name.
Definition matrix_io_mtx.cpp:92
void write_entry(size_t m, size_t n, CPUAlgebra::matrix_type::value_type val)
Appends a data line to the associated file.
Definition matrix_io_mtx.cpp:287
void set_mat_file_name(std::string mFile, int openMode=MatrixIO::EXISTING)
Sets associated file to another one.
Definition matrix_io_mtx.cpp:73
size_t get_num_rows() const
Retrieves number of rows as specified in the exchange file.
Definition matrix_io_mtx.cpp:109
void read_entry(size_t *m, size_t *n, CPUAlgebra::matrix_type::value_type *val)
Reads and parses the next data line.
Definition matrix_io_mtx.cpp:241
std::string * m_pMatFileName
Definition matrix_io_mtx.h:88
size_t get_num_cols() const
Retrieves number of columns as specified in the exchange file.
Definition matrix_io_mtx.cpp:114
size_t m_cols
Definition matrix_io_mtx.h:104
bool is_sparse() const
Tells whether the associated exchange file holds a sparse matrix.
Definition matrix_io_mtx.cpp:124
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define UG_THROW(msg)
Definition error.h:57
#define UG_WARNING(msg)
Definition log.h:328
the ug namespace
#define PROFILE_FUNC()
Definition profiler.h:257