sparsematrix.h

Go to the documentation of this file.

/*
 * Copyright (c) 2010-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__CPU_ALGEBRA__SPARSEMATRIX__
#define __H__UG__CPU_ALGEBRA__SPARSEMATRIX__
 
#include <math.h>
#include "common/common.h"
#include "../algebra_common/sparsematrix_util.h"
#include <iostream>
#include <algorithm>
#include "common/util/ostream_util.h"
 
#include "../algebra_common/connection.h"
#include "../algebra_common/matrixrow.h"
#include "../common/operations_mat/operations_mat.h"
 
#define PROFILE_SPMATRIX(name) PROFILE_BEGIN_GROUP(name, "SparseMatrix algebra")
 
#ifndef NDEBUG
#define CHECK_ROW_ITERATORS
#endif
 
namespace ug{
 
 
 
// example for the variable CRS storage structure:
// say we have:
// rowStart = 0 3 8
// rowEnd = 3 6 11
// rowMax = 3 8 11
// cols ( | marking end of row): 2 5 6 | 2 6 7 x x| 8 9 10
 
// now insert (0 3): row 0 is full (rowEnd[0]==rowMax[0]), copy it to the end, and insert index
// rowStart = 11 3 8
// rowEnd = 15 6 11
// rowMax = 17 8 11
// cols ( | marking end of row): x x x | 2 6 7 x x| 8 9 10 | 2 3 5 6 x x |
 
// now insert (1 3): row 1 not full, we can add it
// rowStart 11 3 8
// rowEnd 15 7 11
// rowMax = 17 8 11
// cols : x x x | 2 3 6 7 x | 8 9 10 | 2 3 5 6 x x |
 
// defragment:
// rowStart 0 4 8
// rowEnd 4 8 11
// rowMax = 4 8 11
// cols : 2 3 5 6 | 2 3 6 7 | 8 9 10
 
 
template<typename TValueType> class SparseMatrix
{
public:
  typedef TValueType value_type;
  enum {rows_sorted=true};
 
  typedef SparseMatrix<value_type> this_type;
 
public:
  typedef AlgebraicConnection<TValueType> connection;
  typedef MatrixRow<this_type> row_type;
  typedef ConstMatrixRow<this_type> const_row_type;
 
public:
  // construction etc
  //----------------------
 
  SparseMatrix();
  virtual ~SparseMatrix () {}
 
  void clear_and_free();
 
  void resize_and_clear(size_t newRows, size_t newCols);
  void resize_and_keep_values(size_t newRows, size_t newCols);
  void clear_retain_structure();
 
  void set_as_transpose_of(const SparseMatrix<value_type> &B, double scale=1.0);
  void set_as_transpose_of2(const SparseMatrix<value_type> &B, double scale=1.0);
 
  void set_as_copy_of(const SparseMatrix<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>
  void axpy(vector_t &dest,
      const number &alpha1, const vector_t &v1,
      const number &beta1, const vector_t &w1) const;
 
  template<typename vector_t>
  void axpy_transposed(vector_t &dest,
      const number &alpha1, const vector_t &v1,
      const number &beta1, const vector_t &w1) const;
 
  template<typename vector_t>
  void apply_ignore_zero_rows(vector_t &dest,
      const number &beta1, const vector_t &w1) const;
 
  template<typename vector_t>
  void apply_transposed_ignore_zero_rows(vector_t &dest,
      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
    {
      axpy(res, 0.0, res, 1.0, x);
      return true;
    }
 
    // 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
    {
      axpy_transposed(res, 0.0, res, 1.0, x);
      return true;
    }
 
    // 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
    {
      axpy(res, 1.0, res, -1.0, x);
      return true;
    }
 
 
 
  inline bool is_isolated(size_t i) const;
 
  void 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
  //----------------------
 
 
 
  inline void check_rc(size_t r, size_t c) const
  {
    UG_ASSERT(r < num_rows() && c < num_cols(), "tried to access element (" << r << ", " << c << ") of " << num_rows() << " x " << num_cols() << " matrix.");
  }
 
  inline void check_row_modifiable(size_t r) const
  {
#ifdef CHECK_ROW_ITERATORS
    UG_ASSERT(nrOfRowIterators[r] == 0, "row " << r << " is used by an iterator and should not be modified.")
#endif
  }
 
 
  void set(double a);
 
  const value_type &operator () (size_t r, size_t c)  const
    {
    check_rc(r, c);
        int j=get_index_const(r, c);
    if(j == -1)
    {
      static value_type v(0.0);
      return v;
    }
        UG_ASSERT(cols[j]==(int)c && j >= rowStart[r] && j < rowEnd[r], "");
        return values[j];
    }
 
  value_type &operator() (size_t r, size_t c)
  {
    check_rc(r, c);
    int j=get_index(r, c);
        UG_ASSERT(j != -1 && cols[j]==(int)c && j >= rowStart[r] && j < rowEnd[r], "");
        return values[j];
    }
 
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);
 
 
  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
  //----------------------
 
  inline size_t num_connections(size_t i) const
  {
    if(rowStart[i] == -1) return 0;
    else return rowEnd[i]-rowStart[i];
  }
 
  size_t num_rows() const { return rowEnd.size(); }
 
  size_t num_cols() const { return m_numCols; }
 
  size_t total_num_connections() const { return nnz; }
 
public:
 
  // Iterators
  //---------------------------
 
  // 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
    {
        SparseMatrix &A;
#ifdef CHECK_ROW_ITERATORS
        size_t _row;
#endif
        size_t i;
    public:
        inline void check() const {
#ifdef CHECK_ROW_ITERATORS
        A.check_row(_row, i);
#endif
      }
        row_iterator(SparseMatrix &_A, size_t row, size_t _i) : A(_A)
#ifdef CHECK_ROW_ITERATORS
         , _row(row)
#endif
         , i(_i) { A.add_iterator(row); }
        row_iterator(row_iterator &&other) : A(other.A),
#ifdef CHECK_ROW_ITERATORS
      _row(other._row),
#endif
      i(other.i) {
#ifdef CHECK_ROW_ITERATORS
        A.add_iterator(other._row);
#else
        A.add_iterator(-1);
#endif
      }
        row_iterator(const row_iterator &other) : A(other.A),
#ifdef CHECK_ROW_ITERATORS
      _row(other._row),
#endif
      i(other.i) {
#ifdef CHECK_ROW_ITERATORS
        A.add_iterator(other._row);
#else
        A.add_iterator(-1);
#endif
      }
        ~row_iterator() {
#ifdef CHECK_ROW_ITERATORS
        A.remove_iterator(_row);
#else
        A.remove_iterator(-1);
#endif
      }
        row_iterator *operator ->() { return this; }
        value_type &value() { check(); return A.values[i];   }
        size_t index() const { check(); return A.cols[i]; }
        bool operator != (const row_iterator &o) const { return i != o.i;  }
        void operator ++ () { ++i; }
    void operator += (int nr) { i+=nr; }
    bool operator == (const row_iterator &other) const { return other.i == i; check(); }
    };
    class const_row_iterator
    {
        const SparseMatrix &A;
#ifdef CHECK_ROW_ITERATORS
        size_t _row; // int?
#endif
        size_t i;
    public:
        inline void check() const {
#ifdef CHECK_ROW_ITERATORS
        A.check_row(_row, i);
#endif
      }
        const_row_iterator(const SparseMatrix &_A, size_t row, size_t _i) : A(_A),
#ifdef CHECK_ROW_ITERATORS
      _row(row),
#endif
      i(_i) {A.add_iterator(row);}
        const_row_iterator(const const_row_iterator &other) : A(other.A),
#ifdef CHECK_ROW_ITERATORS
      _row(other._row),
#endif
      i(other.i) {
#ifdef CHECK_ROW_ITERATORS
        A.add_iterator(_row);
#else
        A.add_iterator(-1);
#endif
      }
        const_row_iterator(const_row_iterator&& other) : A(other.A),
#ifdef CHECK_ROW_ITERATORS
         _row(other._row),
#endif
        i(other.i) {
#ifdef CHECK_ROW_ITERATORS
        A.add_iterator(_row);
#else
        A.add_iterator(-1);
#endif
      }
      const_row_iterator& operator=(const_row_iterator const&) = delete;
      const_row_iterator& operator=(const_row_iterator&&) = delete;
        ~const_row_iterator() {
#ifdef CHECK_ROW_ITERATORS
        A.remove_iterator(_row);
#else
        A.remove_iterator(-1);
#endif
      }
        const const_row_iterator *operator ->() const { return this; }
        const value_type &value() const { check(); return A.values[i];   }
        size_t index() const { check(); return A.cols[i];     }
        bool operator != (const const_row_iterator &o) const { return i != o.i; }
        void operator ++ () { ++i; }
        void operator += (int nr) { i+=nr; }
    bool operator == (const const_row_iterator &other) const { return other.i == i; }
    public: // BUG
     // int row() const{return _row;}
     size_t idx() const{return i;}
    };
 
 
 
 
  row_iterator         begin_row(size_t r)         { return row_iterator(*this, r, rowStart[r]);  }
    row_iterator         end_row(size_t r)           { return row_iterator(*this, r, rowEnd[r]);  }
    const_row_iterator   begin_row(size_t r) const   { return const_row_iterator(*this, r, rowStart[r]);  }
    const_row_iterator   end_row(size_t r)   const   { return const_row_iterator(*this, r, rowEnd[r]);  }
 
    row_type    get_row(size_t r)     { return row_type(*this, r); }
    const_row_type  get_row(size_t r) const { return const_row_type(*this, r); }
 
public:
  // connectivity functions
  //-------------------------
 
    bool has_connection(size_t r, size_t c) const
    {
      check_rc(r, c);
      bool bFound;
      get_connection(r, c, bFound);
      return bFound;
    }
 
  row_iterator get_iterator_or_next(size_t r, size_t c)
  {
    check_rc(r, c);
    if(rowStart[r] == -1 || rowStart[r] == rowEnd[r])
          return end_row(r);
        else
        {
          int j=get_index_internal(r, c);
          if(j > maxValues) return end_row(r);
          else return row_iterator(*this, r, j);
        }
    }
 
  const_row_iterator get_connection(size_t r, size_t c, bool &bFound) const
  {
    check_rc(r, c);
    int j=get_index_const(r, c);
    if(j != -1)
    {
      bFound = true;
      return const_row_iterator(*this, r, j);
    }
    else
    {
      bFound = false;
      return end_row(r);
    }
    }
  row_iterator get_connection(size_t r, size_t c, bool &bFound)
  {
    check_rc(r, c);
    int j=get_index_const(r, c);
    if(j != -1)
    {
      bFound = true;
      return row_iterator(*this, r, j);
    }
    else
    {
      bFound = false;
      return end_row(r);
    }
  }
 
  const_row_iterator get_connection(size_t r, size_t c) const
  {
    bool b;
    return get_connection(r, c, b);
  }
  row_iterator get_connection(size_t r, size_t c)
  {
    check_rc(r, c);
    assert(bNeedsValues);
    int j=get_index(r, c);
    return row_iterator(*this, r, j);
  }
 
 
  void defragment()
    {
    if(num_rows() != 0 && num_cols() != 0)
      copyToNewSize(nnz);
    }
 
  void defragment() const
  {
    (const_cast<this_type*>(this))->defragment();
  }
 
  void copy_crs(size_t &numRows, size_t &numCols,
      std::vector<value_type> &argValues, std::vector<int> &argRowStart,
      std::vector<int> &argColInd) const
  {
    numRows = num_rows();
    numCols = num_cols();
    defragment();
    argValues = values;
    argRowStart = rowStart;
    argColInd = cols;
  }
 
  void get_crs(size_t &numRows, size_t &numCols,
      value_type *&pValues, size_t *pRowStart, size_t *pColInd, size_t &nnz) const
  {
    defragment();
    pValues = &values[0];
    // FIXME:
    //pRowStart = &rowStart[0]; // assigning int * to size_t *
    //pColInd = &cols[0];       // assigning int * to size_t *
    UG_THROW("SparseMatrix::get_crs() needs to be fixed.");
    numRows = num_rows();
    numCols = num_cols();
    nnz = total_num_connections();
  }
 
  void get_values(std::vector<value_type> &argValues) const
  {
    defragment();
    argValues = values;
  }
 
 
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 SparseMatrix &m)
  {
    out << "SparseMatrix " //<< 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:
  // private functions
 
  void add_iterator(size_t row) const
  {
#ifdef CHECK_ROW_ITERATORS
    nrOfRowIterators[row]++;
#endif
    iIterators++;
  }
  void remove_iterator(size_t row) const
  {
#ifdef CHECK_ROW_ITERATORS
    nrOfRowIterators[row]--;
    UG_ASSERT(nrOfRowIterators[row] >= 0, row);
#endif
    iIterators--;
    UG_ASSERT(iIterators >= 0, row);
 
  }
  inline void check_row(size_t row, int i) const
  {
    UG_ASSERT(i < rowEnd[row] && i >= rowStart[row], "row iterator row " << row << " pos " << i << " out of bounds [" << rowStart[row] << ", " << rowEnd[row] << "]");
  }
    void assureValuesSize(size_t s);
    size_t get_nnz() const { return nnz; }
 
private:
  // disallowed operations (not defined):
  //---------------------------------------
  SparseMatrix(SparseMatrix&); 
 
 
protected:
  int get_index_internal(size_t row, int col) const;
    int get_index_const(int r, int c) const;
    int get_index(int r, int c);
    void copyToNewSize(size_t newSize)
    {
      copyToNewSize(newSize, num_cols());
    }
    void copyToNewSize(size_t newSize, size_t maxCols);
  void check_fragmentation() const;
  int get_nnz_max_cols(size_t maxCols);
 
public: // bug
  int col(size_t i) const{
    assert(i<cols.size());
    return cols[i];
  }
 
#ifndef NDEBUG
  int iters() const{
    return iIterators;
  }
#endif
 
protected:
    std::vector<int> rowStart;
    std::vector<int> rowEnd;
    std::vector<int> rowMax;
    std::vector<int> cols;
    size_t fragmented;
    size_t nnz;
    bool bNeedsValues;
 
    std::vector<value_type> values;
    int maxValues;
    int m_numCols;
    mutable int iIterators;
 
#ifdef CHECK_ROW_ITERATORS
public:
    mutable std::vector<int> nrOfRowIterators;
#endif
};
 
 
template<typename T>
struct matrix_algebra_type_traits<SparseMatrix<T> >
{
  enum{
    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 SparseMatrix<matrix_t> &A1, const vector_t &w1)
{
  A1.axpy_transposed(dest, alpha1, v1, beta1, w1);
}
 
 
 
 
 
 
 
 
// end group cpu_algebra
 
} // namespace ug
 
//#include "matrixrow.h"
#include "sparsematrix_impl.h"
#include "sparsematrix_print.h"
 
#endif