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 QVector<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); }
 
  // See comments on the base class
  virtual bool hasNanInf() const { return _mat.has_inf() || _mat.has_nan(); }
 
  // This function is currently not implemented for Full matrices. If we really need
  // this, we have to remove it from here and add it to GmArmadilloSolverMatrix (that
  // will have to be transformed from a typedef into a real class). Optionally, we could 
  // use C++ 20 concepts to have an implementation of this function for T == arma::mat 
  // only
  virtual bool saveMatrixCoordinatesToFile(FILE* f) const { Q_UNUSED(f); return false; }
 
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 QVector<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;
 
  virtual bool saveMatrixCoordinatesToFile(FILE* f) 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