mapsparsematrix.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2015:  G-CSC, Goethe University Frankfurt
 * Author: Martin Rupp
 * 
 * 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__MAP_ALGEBRA__SPARSEMATRIX__
#define __H__UG__MAP_ALGEBRA__SPARSEMATRIX__
 
#include "math.h"
#include "common/common.h"
 
#include "../common/operations_mat/operations_mat.h"
 
namespace ug{
 
 
 
template<typename TValueType> class MapSparseMatrix
{
public:
  typedef TValueType value_type;
  enum {rows_sorted=true};
 
  typedef SparseMatrix<value_type> this_type;
 
  
public:
  struct connection
  {
    size_t iIndex;    // index to
    value_type dValue; // smallmatrix value;
    connection() {}
    connection(size_t i, const value_type &v)
    : iIndex(i), dValue(v) {}
        
    void print(){std::cout << *this;}
    friend std::ostream &operator<<(std::ostream &output, const connection &c)
    {
      output << "(" << c.iIndex << "-> ";
      output << c.dValue;
      output << ")";
      return output;
    }
 
    void operator = (const connection &other)
    {
      iIndex = other.iIndex;
      dValue = other.dValue;
    }
 
    int operator < (const connection &c) const
    {
      return iIndex < c.iIndex;
    }
  };
 
public:
  // construction etc
  //----------------------
 
  MapSparseMatrix();
  virtual ~MapSparseMatrix () {}
 
 
  bool resize(size_t newRows, size_t newCols);
 
 
  bool set_as_transpose_of(const SparseMatrix<value_type> &B, double scale=1.0);
 
  bool set_as_copy_of(const MapSparseMatrix<value_type> &B, double scale=1.0);
  SparseMatrix<value_type> &operator = (const SparseMatrix<value_type> &B)
  {
    set_as_copy_of(B);
    return *this;
  }
 
 
public:
  template<typename vector_t>
  bool axpy(vector_t &dest,
      const number &alpha1, const vector_t &v1,
      const number &beta1, const vector_t &w1) const;
 
  template<typename vector_t>
  bool axpy_transposed(vector_t &dest,
      const number &alpha1, const vector_t &v1,
      const number &beta1, const vector_t &w1) const;
 
 
  // DEPRECATED!
    // apply is deprecated because of axpy(res, 0.0, res, 1.0, beta, w1)
    template<typename Vector_type>
    bool apply(Vector_type &res, const Vector_type &x) const
    {
      return axpy(res, 0.0, res, 1.0, x);
    }
 
    // apply is deprecated because of axpy(res, 0.0, res, 1.0, beta, w1)
    template<typename Vector_type>
    bool apply_transposed(Vector_type &res, const Vector_type &x) const
    {
      return axpy_transposed(res, 0.0, res, 1.0, x);
    }
 
    // matmult_minus is deprecated because of axpy(res, 1.0, res, -1.0, x);
    template<typename Vector_type>
    bool matmul_minus(Vector_type &res, const Vector_type &x) const
    {
      return axpy(res, 1.0, res, -1.0, x);
    }
 
 
 
  inline bool is_isolated(size_t i) const;
 
  bool scale(double d);
  SparseMatrix<value_type> &operator *= (double d) { scale(d); return *this; }
 
  // submatrix set/get functions
  //-------------------------------
 
  template<typename M>
  void add(const M &mat);
  template<typename M>
  void set(const M &mat);
  template<typename M>
  void get(M &mat) const;
 
  // finalizing functions
  //----------------------
 
 
 
  
  bool set(double a);
 
  const value_type &operator () (size_t r, size_t c)  const
    {
        typename std::map<size_t, TValueType>::const_iterator it = data[r].find(c);
    if(it != data[r].end(c))
    {
      static value_type v(0.0);
      return v;
    }
    return (*it).second;
    
    }
 
  value_type &operator() (size_t r, size_t c)
  {
    return data[r][c];
    }
 
public:
  // row functions
 
  void set_matrix_row(size_t row, connection *c, size_t nr);
 
  void add_matrix_row(size_t row, connection *c, size_t nr);
 
  inline size_t num_connections(size_t i) const
  {
    return data[i].size();
  }
 
  template<typename vector_t>
  inline void mat_mult_add_row(size_t row, typename vector_t::value_type &dest, double alpha, const vector_t &v) const;
public:
  // accessor functions
  //----------------------
 
  size_t num_rows() const { return data.size(); }
  
  size_t num_cols() const { return cols; }
 
  size_t total_num_connections() const
  { 
    size_t nnz=0;
    for(size_t i=0; i<data.size(); i++)
      nnz += num_connections(i);
    return nnz;
  }
 
public:
 
  // Iterators
  //---------------------------
 
  // const_row_iterator
 
 
  //typedef const connection * const_row_iterator;
  //typedef connection * const_row_iterator;
  // a row_iterator has to suppport
  // operator ++, operator +=, index() const, const value_type &value() const, value_type &value()
  // a const_row_iterator has to suppport
  // operator ++, operator +=, index() const, const value_type &value() const
 
  class row_iterator : public std::map<size_t, TValueType>::iterator
    {
    typedef typename std::map<size_t, TValueType>::iterator super;    
    public:        
    row_iterator(super it) : super(it) {}        
        value_type &value() { return super::operator*().second;   }
        size_t index() const { return super::operator*().first;   }        
    };
    class const_row_iterator : public std::map<size_t, TValueType>::const_iterator
    {
    typedef typename std::map<size_t, TValueType>::const_iterator super;
  public:        
        const_row_iterator(super it) : super(it) {}        
    const value_type &value() { return super::operator*().second;   }
        size_t index() const { return super::operator*().first;   }        
    };
  
 
  row_iterator         begin_row(size_t r)         { return row_iterator(data[r].begin());  }
    row_iterator         end_row(size_t r)           { return row_iterator(data[r].end());  }
    const_row_iterator   begin_row(size_t r) const   { return const_row_iterator(data[r].begin());  }
    const_row_iterator   end_row(size_t r)   const   { return const_row_iterator(data[r].end());  }
  
public:
  // connectivity functions
  //-------------------------
 
 
  const_row_iterator get_connection(size_t r, size_t c, bool &bFound) const
  {
    const_row_iterator it = const_row_iterator(data[r].find(c));
    bFound = (it != end_row(r));
    return it;
    }
  row_iterator get_connection(size_t r, size_t c, bool &bFound)
  {
    row_iterator it = row_iterator(data[r].find(c));
    bFound = (it != end_row(r));
    return it;
  }
 
  const_row_iterator get_connection(size_t r, size_t c) const
  {
    return const_row_iterator(data[r].find(c));
  }
  row_iterator get_connection(size_t r, size_t c)
  {
    return row_iterator(data[r].find(c));
  }
  
  
  void defragment()
    {
        
    }
 
public:
  // output functions
  //----------------------
 
  void print(const char * const name = NULL) const;
  void printtype() const;
 
  void print_to_file(const char *filename) const;
  void printrow(size_t row) const;
 
  friend std::ostream &operator<<(std::ostream &out, const MapSparseMatrix &m)
  {
    out << "MapSparseMatrix " //<< m.name
    << " [ " << m.num_rows() << " x " << m.num_cols() << " ]";
    return out;
  }
 
 
  void p() const { print(); } // for use in gdb
  void pr(size_t row) const {printrow(row); } // for use in gdb
 
private:
  // disallowed operations (not defined):
  //---------------------------------------
  MapSparseMatrix(MapSparseMatrix&); 
 
  
protected:
    std::vector<std::map<size_t, TValueType> > data;
  size_t cols;
};
 
 
template<typename T>
struct matrix_algebra_type_traits<MapSparseMatrix<T> >
{
  static const int type = MATRIX_USE_ROW_FUNCTIONS;
};
 
template<typename vector_t, typename matrix_t>
inline void MatMultTransposedAdd(vector_t &dest,
    const number &alpha1, const vector_t &v1,
    const number &beta1, const MapSparseMatrix<matrix_t> &A1, const vector_t &w1)
{
  A1.axpy_transposed(dest, alpha1, v1, beta1, w1);
}
 
} // namespace ug
 
//#include "matrixrow.h"
#include "mapsparsematrix_impl.h"
#include "mapsparsematrix_print.h"
 
#endif