matrix_io_mtx.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2013:  G-CSC, Goethe University Frankfurt
 * Author: Torbjörn Klatt
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
#ifndef __H__UG__LIB_ALGEBRA__MATRIX_IO_MTX_H
#define __H__UG__LIB_ALGEBRA__MATRIX_IO_MTX_H
 
#include <sstream>
#include <vector>
#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
 
#include "lib_algebra/common/matrixio/matrix_io.h"
#include "lib_algebra/common/matrixio/mm_type_code.h"
 
namespace ug
{
 
 
class MatrixIOMtx : private MatrixIO
{
  private:
    // Full path name of the matrix exchange file
    // (docu in matrix_io.h)
    std::string *m_pMatFileName;
    // Internal file stream for reading from and writing into the matrix
    // exchange file (docu in matrix_io.h)
    std::fstream m_matFileStream;
    size_t m_firstDataLine;
    // Matrix exchange file type (set to constant MatrixFileType::MatrixMarket)
    // (docu in matrix_io.h)
    MatrixFileType m_matFileType;
    MMTypeCode m_mmTypeCode;
    // Number of rows as specified in the matrix exchange file
    // (docu in matrix_io.h)
    size_t m_rows;
    // Number of columns as specified in the matrix exchange file
    // (docu in matrix_io.h)
    size_t m_cols;
    // Number of data lines as specified in the matrix exchange file
    // (docu in matrix_io.h)
    size_t m_lines;
 
  public:
    MatrixIOMtx();
    MatrixIOMtx( std::string mFile, int openMode=MatrixIO::EXISTING );
    
    ~MatrixIOMtx();
 
    void set_mat_file_name( std::string mFile, int openMode=MatrixIO::EXISTING );
    
    std::string get_mat_file_name() const;
 
    void set_mat_dims( size_t rows, size_t cols, size_t lines );
 
    size_t get_num_rows() const;
 
    size_t get_num_cols() const;
 
    size_t get_num_lines() const;
 
    bool is_sparse() const;
 
    template<typename matrix_type>
    void read_into( matrix_type &matrix )
    {
      PROFILE_FUNC();
      UG_ASSERT( m_rows > 0 || m_cols > 0 || m_lines > 0,
                 "MatrixMarket matrix dimensions seem not be reasonable: (m, n, nnz)=("
                 << m_rows << "," << m_cols << "," << m_lines << ")" );
 
      // open the file and go one before first data line
      open_file();
 
      std::string dummy;
      for( size_t i = 0; i < m_firstDataLine; i++ ) {
        std::getline( m_matFileStream, dummy );
      }
 
      matrix.resize_and_clear( m_rows, m_cols );
      size_t x, y;
      double val;
      if ( m_mmTypeCode.is_sparse() ) {
        for( size_t i = 0; i < m_lines; i++ ) {
          read_entry( &x, &y, &val );
          // we need to substract 1 from row and column indices as MM is 1-based
          matrix( x - 1, y - 1 ) = val;
          if ( m_mmTypeCode.is_symmetric() && x != y ) {
            matrix( y - 1, x - 1 ) = val;
          } else if ( m_mmTypeCode.is_skew_symmetric() && x != y ) {
            matrix( y - 1, x - 1 ) = -val;
          }
        }
      } else {
        UG_THROW( "Other than sparse MatrixMarket matrices are not yet implemented." );
      }
 
      // close the file
      close_file();
    }
 
 
    template<typename matrix_type>
    void write_from( matrix_type &matrix,
                     std::string comment="%Generated with ug4." )
    {
    PROFILE_FUNC();
 
    // analyze the given matrix
    std::vector< std::vector<size_t> > rowIndexPerCol = determine_matrix_characteristics( matrix );
 
    // write the MatrixMarket banner
    write_banner();
 
    // open the file for appending
    open_file( std::ios_base::out | std::ios_base::app );
 
    // add a comment if it's not empty
    if ( !comment.empty() ) {
    if ( comment.find_first_of( '%' ) != 0 ) {
      UG_WARNING( "Given comment did not start with '%'. Prepending it to make it valid." );
      comment.insert( 0, "%" );
    }
    m_matFileStream << comment << "\n";
    }
 
    // write characteristics
    m_matFileStream << m_rows << " " << m_cols << " " << m_lines << "\n";
 
 
    // write entries to the file
    for ( size_t col = 0; col < m_cols; col++ ) {
    for (size_t row = 0; row < rowIndexPerCol.at(col).size(); row++ ) {
      // we need to add 1 to row and column index, as MM indices are 1-based
      write_entry( rowIndexPerCol.at(col).at(row) + 1, col + 1,
             matrix(rowIndexPerCol.at(col).at(row), col) );
    }
    }
    close_file();
    }
 
 
  private:
    // Opens file as fstream. (docu in matrix_io.h)
    void open_file( std::ios_base::openmode mode=std::ios_base::in );
    
    // Closes fstream. (docu in matrix_io.h)
    void close_file();
 
    void query_matrix_type();
 
    void query_matrix_characteristics();
 
    template<typename matrix_type>
    std::vector< std::vector<size_t> > determine_matrix_characteristics(const matrix_type &matrix )
    {
      PROFILE_FUNC();
 
      // read matrix dimensions
      size_t rows = matrix.num_rows();
      size_t cols = matrix.num_cols();
      size_t offDiagEntries = 0;
      size_t diagEntries = 0;
 
      // As MatrixMarket specifies a column-first-ordering, we need to get to know
      // in which rows of each column non-zero entries are.
      // During this, we also get to know how many non-zero entries there are in
      // total and can detect symmetries of the matrix.
      std::vector< std::vector<size_t> > rowIndexPerCol;
      rowIndexPerCol.resize(cols);
      bool isSymmetric = true;
      bool isSkew = true;
      bool changed = false;
      for ( size_t r = 0; r < rows; r++ ) { // iterate rows
        for (typename matrix_type::const_row_iterator conn = matrix.begin_row(r);
              conn != matrix.end_row(r); ++conn ) {
          if ( conn.value() != 0.0 ) {
            // first add index to list
            if ( isSymmetric || isSkew ) {
              if ( conn.index() <= r ) {
                rowIndexPerCol.at( conn.index() ).push_back(r);
              }
            } else {
              rowIndexPerCol.at( conn.index() ).push_back(r);
            }
 
            // increment counters
            (conn.index() == r) ? diagEntries++ : offDiagEntries++ ;
 
            // check, whether it's still a symmetric or skew-symmetric matrix
            if ( r != conn.index() ) {
              if ( matrix(r, conn.index() ) != matrix( conn.index(), r ) ) {
                if ( isSymmetric ) {
                  changed = true;
                }
                isSymmetric = false;
              }
              if ( matrix(r, conn.index() ) != -1.0 * matrix( conn.index(), r ) ) {
                if ( isSkew ) {
                  changed = true;
                }
                isSkew = false;
              }
            }
 
            // We assumed the matrix to be symmetric or skew-symmetric, but it's
            // not. Thus we need to redo everything done before.
            if ( changed ) {
              offDiagEntries = 0;
              diagEntries = 0;
              rowIndexPerCol.clear();
              rowIndexPerCol.resize( cols );
              r = -1; // at the end of the current row-loop this is incremented again
              changed = false;
              break;
            }
          }
        }
      }
 
      // make sure the matrix is not both, symmetric and skew-symmetric
      UG_ASSERT( ( ( isSymmetric && !isSkew ) ||
                   ( isSkew && !isSymmetric ) ||
                   ( !isSkew && !isSymmetric ) ),
                 "Error on detecting symmetry of matrix. (skew:" << isSkew
                 << " sym:" << isSymmetric << ")" );
 
      // set MMTypeCode
      size_t entries = 0;
      m_mmTypeCode.set_class_type( MMTypeCode::COORDINATE );
      m_mmTypeCode.set_numeric_type( MMTypeCode::REAL );
      if ( isSymmetric ) {
        entries = ( offDiagEntries / 2 )+ diagEntries;
        m_mmTypeCode.set_algebraic_type( MMTypeCode::SYMMETRIC );
      } else if ( isSkew ) {
        entries = ( offDiagEntries / 2 ) + diagEntries;
        m_mmTypeCode.set_algebraic_type( MMTypeCode::SKEW );
      } else {
        entries = offDiagEntries + diagEntries;
        m_mmTypeCode.set_algebraic_type( MMTypeCode::GENERAL );
      }
 
      // now we can set the dimensions
      set_mat_dims( rows, cols, entries );
 
      return rowIndexPerCol;
    }
 
    void read_entry( size_t *m, size_t *n, CPUAlgebra::matrix_type::value_type *val );
 
    void write_banner();
    
    void write_entry( size_t m, size_t n, CPUAlgebra::matrix_type::value_type val );
};
 
// end group matrixio
 
} // namespace ug
 
#endif // __H__UG__LIB_ALGEBRA__MATRIX_IO_MTX_H
 
// EOF