gmArmadilloSolverMatrix.h

Go to the documentation of this file.

/************************************************************************
**
** Copyright (C) 2014 by Carlos Augusto Teixera Mendes
** All rights reserved.
**
** This file is part of the "GeMA" software. It's use should respect
** the terms in the license agreement that can be found together
** with this source code.
** It is provided AS IS, with NO WARRANTY OF ANY KIND,
** INCLUDING THE WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR
** A PARTICULAR PURPOSE.
**
************************************************************************/

#ifndef _GEMA_ARMADILLO_SOLVER_MATRIX_H_
#define _GEMA_ARMADILLO_SOLVER_MATRIX_H_

#include "gmSolverMatrix.h"
#include "gmThreadManager.h"
#include "gmSparseMatrix.h"

#include <assert.h>

#include <armadillo>

template <class T> class GmAppendBuffer;

template <class T> class GmArmadilloSolverMatrixBase : public GmSolverMatrix
{
public:
  GmArmadilloSolverMatrixBase(int nlin, int ncol, GmSparseMatrixOptions options = GmSparseMatrixOptions(GM_ARMADILLO_TRIPLET_LIST))
    : _mat(nlin, ncol)
  {
    Q_UNUSED(options);
    assert(GmThreadManager::inMainThread());

    _mat.zeros();
    _sym = false;
  }

  virtual ~GmArmadilloSolverMatrixBase() {}

  // See comments on the base class
  virtual bool supportsBatchInsert() const { return false; }

  // See comments on the base class
  virtual bool supportsSparseLayouts() const { return false; }

  // See comments on the base class
  virtual bool supportsRandomSet() const { return true; }

  // See comments on the base class
  virtual ParallelAddMode supportedParallelAddMode() const { return REENTRANT_SUPPORT; }

  // See comments on the base class
  virtual bool beginBatchInsert(size_t expectedEntries = 0) { Q_UNUSED(expectedEntries); assert(GmThreadManager::inMainThread()); return true; }

  // See comments on the base class
  virtual bool endBatchInsert(bool discardData = false) { Q_UNUSED(discardData);  assert(GmThreadManager::inMainThread()); return true; }

  // See comments on the base class
  virtual GmSparseMatrixLayoutBuilder* layoutBuilder() const { assert(GmThreadManager::inMainThread()); return NULL; }

  // See comments on the base class
  virtual bool emptyLayout() const { assert(GmThreadManager::inMainThread()); return false; }

  // See comments on the base class
  virtual bool symmetric() const { return _sym.loadAcquire(); }

  // See comments on the base class
  virtual void setSymmetric(bool sym) { _sym.storeRelease(sym); }

  // See comments on the base class
  virtual int nlin() const { return _mat.n_rows; }

  // See comments on the base class
  virtual int ncol() const { return _mat.n_cols; }

  // See comments on the base class
  virtual int layoutSize() const { assert(GmThreadManager::inMainThread()); return _mat.n_elem; }

  // See comments on the base class
  virtual double at(int lin, int col) const { assert(GmThreadManager::inMainThread()); return _mat(lin, col); }

  // See comments on the base class
  virtual bool inLayout(int lin, int col) const { assert(GmThreadManager::inMainThread()); return _mat(lin, col) != 0.0; }

  // See comments on the base class
  virtual void set(int lin, int col, double val) { assert(GmThreadManager::inMainThread()); _mat(lin, col) = val; }

  // See comments on the base class
  virtual void add(int lin, int col, double val) { _mat(lin, col) += val; }

  // See comments on the base class
  virtual void set(GmSolverMatrix* A, GmSolverMatrix* B, double c, bool sameStructure)
  {
    Q_UNUSED(sameStructure);
    assert(GmThreadManager::inMainThread());
    assert(A->nlin() == _mat.n_rows && A->ncol() == _mat.n_cols);
    assert(B->nlin() == _mat.n_rows && B->ncol() == _mat.n_cols);

    GmArmadilloSolverMatrixBase<T>* newA = dynamic_cast<GmArmadilloSolverMatrixBase<T>*>(A);
    GmArmadilloSolverMatrixBase<T>* newB = dynamic_cast<GmArmadilloSolverMatrixBase<T>*>(B);
    if(newA && newB)
      _mat = newA->_mat + c * newB->_mat;
    else
    {
      for(unsigned int j = 0; j<_mat.n_cols; j++)  // In Armadillo, values are stored in column order
        for(unsigned int i = 0; i<_mat.n_rows; i++)
          _mat(i, j) = A->at(i, j) + c * B->at(i, j);
    }

    if(symmetric() && (!A->symmetric() || !B->symmetric()))
      setSymmetric(false);
  }

  // See comments on the base class
  virtual void clear(bool keepSparseLayout) 
  { 
    Q_UNUSED(keepSparseLayout);  
    assert(GmThreadManager::inMainThread()); 
    _mat.zeros(); 
    setSymmetric(false);
  }

  // See comments on the base class
  virtual void clearLineAndColumnSet(const QList<int>& indexList, bool setDiagonal, bool keepSparseLayout)
  {
    Q_UNUSED(keepSparseLayout);
    assert(GmThreadManager::inMainThread());
    foreach(int i, indexList)
    {
      _mat.row(i).zeros();
      _mat.col(i).zeros();
      if(setDiagonal)
        _mat(i, i) = 1.0;
    }
  }

  // See comments on the base class
  virtual void ensureDiagonal()
  {
    assert(GmThreadManager::inMainThread());
    for(int i = 0; i < _mat.n_rows; i++)
    {
      if(_mat(i, i) == 0.0)
        _mat(i, i) = 1.0;
    }
  }

  // See comments on the base class
  virtual void matAdd(const GmVector& a, const GmVector& b, double zeroTol = 0.0)
  {
    assert(GmThreadManager::inMainThread());
    assert(a.n_rows == _mat.n_rows);
    assert(b.n_rows == _mat.n_cols);
    for(int j = 0; j < _mat.n_cols; j++)  // In Armadillo, values are stored in column order
    {
      double bv = b[j];
      for(int i = 0; i < _mat.n_rows; i++)
      {
        double v = _mat(i, j) + a[i] * bv;
        if(fabs(v) > zeroTol)
          _mat(i, j) = v;
      }
    }
  }

  // See comments on the base class
  virtual void mul(const GmVector& a, GmVector& b) const
  {
    assert(GmThreadManager::inMainThread());
    assert(a.n_elem == _mat.n_cols);
    b = _mat * a;
  }

  // See comments on the base class
  virtual void mulAdd(const GmVector& a, GmVector& b) const
  {
    assert(GmThreadManager::inMainThread());
    assert(a.n_elem == _mat.n_cols);
    b += _mat * a;
  }

  // See comments on the base class
  virtual void mulSub(const GmVector& a, GmVector& b) const
  {
    assert(GmThreadManager::inMainThread());
    assert(a.n_elem == _mat.n_cols);
    b -= _mat * a;
  }

  // See comments on the base class
  virtual void columnMulAdd(int col, GmVector& f, double v, const bool* skipRows = NULL) const
  {
    assert(GmThreadManager::inMainThread());
    assert(f.n_elem == _mat.n_rows);
    for(int i = 0; i < _mat.n_rows; i++)
    {
      if(skipRows && skipRows[i])
        continue;
      f[i] += _mat(i, col) * v;
    }
  }

  T& armadilloMatrix() { return _mat; }

  // See comments on the base class
  virtual size_t usedMemory() const { return _mat.n_rows * _mat.n_cols * sizeof(double); }

private:
  Q_DISABLE_COPY(GmArmadilloSolverMatrixBase);

protected:
  T          _mat;  
  QAtomicInt _sym;  
};

typedef GmArmadilloSolverMatrixBase<arma::mat> GmArmadilloSolverMatrix;

class GMC_API_EXPORT GmArmadilloSparseSolverMatrix : public GmArmadilloSolverMatrixBase<arma::sp_mat>
{
public:
  GmArmadilloSparseSolverMatrix(int nlin, int ncol, GmSparseMatrixOptions options);

  virtual ~GmArmadilloSparseSolverMatrix();

  // See comments on the base class
  virtual bool supportsBatchInsert() const { return true; }

  // See comments on the base class
  virtual ParallelAddMode supportedParallelAddMode() const { return THREAD_SAFE_SUPPORT; }

  virtual bool beginBatchInsert(size_t expectedEntries = 0);
  virtual bool endBatchInsert(bool discardData = false);

  // See comments on the base class
  virtual bool supportsRandomSet() const { return true; }

  // See comments on the base class
  virtual int layoutSize() const { assert(GmThreadManager::inMainThread()); return _mat.n_nonzero; }

  virtual void add(int lin, int col, double val);

  virtual void clearLineAndColumnSet(const QList<int>& indexList, bool setDiagonal, bool keepSparseLayout);

  virtual void matAdd(const GmVector& a, const GmVector& b, double zeroTol = 0.0);

  /* This is the same implementation as in the base class.
     Repeated here to add an assert and to circumvent some weird compiler errors where a call to 
     set(int, int, double) would try to call set(GmSolverMatrix*, GmSolverMatrix*, double)
  */
  virtual void set(int lin, int col, double val) { assert(GmThreadManager::inMainThread()); assert(!_batch);  _mat(lin, col) = val; }

  virtual void set(GmSolverMatrix* A, GmSolverMatrix* B, double c, bool sameStructure);

  virtual void columnMulAdd(int col, GmVector& f, double v, const bool* skipRows = NULL) const;

  virtual size_t usedMemory() const;

private:
  bool doEnd(bool ordered);

  virtual void addFromThread(int tid, int lin, int col, double val);

protected:
  GmSparseMatrixOptions _opt;   
  bool                  _batch; 

  GmAppendBuffer<arma::uword>* _positions;  
  GmAppendBuffer<double>*      _values;     
};

#endif