gmSparseMatrixLayout.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_SPARSE_MATRIX_LAYOUT_H_
#define _GEMA_SPARSE_MATRIX_LAYOUT_H_

#include "gmCoreConfig.h"

#include <armadillo>
#include <QObject>

#if 0
#define ENABLE_MATRIX_LAYOUT_INDEX_CACHE 1
#endif

#ifdef ENABLE_MATRIX_LAYOUT_INDEX_CACHE
#include "gmThreadLocalStorage.h"
#endif


class GmLogCategory;

enum GmSparseMatrixLayoutTypes
{
  GM_CSC_SPARSE_FORMAT,  
  GM_CSR_SPARSE_FORMAT,  
};


struct GMC_API_EXPORT GmSparseMatrixLayout : public QObject
{
  Q_OBJECT

public:
  GmSparseMatrixLayout(GmSparseMatrixLayoutTypes type) : _type(type), _shared(false) {}

  virtual ~GmSparseMatrixLayout() {}

  virtual bool empty() const = 0;

  virtual void clear() = 0;

  virtual int index(int row, int col) const = 0;

  virtual int dindex(int n) const = 0;

  virtual void printStatistics(const GmLogCategory& logger) const = 0;

  void emitLayoutCompleted() { emit layoutCompleted(); }

  GmSparseMatrixLayoutTypes _type;   
  mutable bool              _shared; 

#ifdef ENABLE_TESTS
  virtual void checkLayoutData() const = 0;
#endif

signals:
  void layoutCompleted();

private:
  Q_DISABLE_COPY(GmSparseMatrixLayout)
};


template <class IndexType, GmSparseMatrixLayoutTypes T> 
struct GMC_API_EXPORT GmCSxSparseMatrixLayout: public GmSparseMatrixLayout
{
public:
  GmCSxSparseMatrixLayout(bool arma) : GmSparseMatrixLayout(T)
  {
    _ptr  = _index = NULL;
    _n    = _nnz   = 0;
    _arma = arma;

#ifdef ENABLE_MATRIX_LAYOUT_INDEX_CACHE
    _lastIndex.init(-1);
#endif
  }

  virtual ~GmCSxSparseMatrixLayout()
  {
    delete[] _ptr;
    delete[] _index;
  }

  // See comments on the base class
  virtual bool empty() const { return _n == 0; }

  // See comments on the base class
  virtual void clear()
  {
    delete[] _ptr;
    delete[] _index;
    _ptr = _index = NULL;
    _n   = _nnz   = 0;

#ifdef ENABLE_MATRIX_LAYOUT_INDEX_CACHE
    _lastIndex.init(-1);
#endif
  }

  // Index implementation is specialized explicitly in the cpp file to avoid the need of an 'if' checking 
  // the layout type.  Although this is static and the compiler should be able to optimize the if, looking
  // at the generated optimized code does not clarify if the compiler is really doing that :(.
  virtual int index(int row, int col) const;

  // See comments on the base class.  TODO: Should we try to optimize searching first in the center values?
  virtual int dindex(int n) const { return indexWorker(_ptr[n], _ptr[n+1], n); }

  virtual void printStatistics(const GmLogCategory& logger) const;

#ifdef ENABLE_TESTS
  virtual void checkLayoutData() const;
#endif

  IndexType* _ptr;   
  IndexType* _index; 
  int        _n;     
  int        _nnz;   
  int        _arma;  

#ifdef ENABLE_MATRIX_LAYOUT_INDEX_CACHE
  mutable GmTLS<int, true> _lastIndex;  
#endif

protected:
  int indexWorker(int start, int end, int ind) const
  {
    // Simple linear search. Number of non zeros should be small!
    // Should we change for a binary search if the array length gets past a limit?
    // Should we unroll part of the loop?  Should we try to use vector instructions?
    // Should we remember the last search and begin from there?
    // TODO: Optimize.
    // Interesting links:
    //    - https://dirtyhandscoding.wordpress.com/2017/08/25/performance-comparison-linear-search-vs-binary-search/
    //    - https://schani.wordpress.com/2010/04/30/linear-vs-binary-search/
    //    - https://blog.demofox.org/2017/06/20/simd-gpu-friendly-branchless-binary-search/
    //    - https://www.agner.org/optimize/
#if 1
    // Standard
    int i = start;
    for(; i < end; ++i)
    {
      if(_index[i] >= ind)
        return (_index[i] == ind) ? i : -1;
    }
    return -1;
#endif
#if 0
    // Cache
    int& li = _lastIndex.localData();
    int i = (li > start && li < end && _index[li] <= ind) ? li : start;
    for(; i < end; ++i)
    {
      if(_index[i] >= ind)
      {
        li = i;
        return (_index[i] == ind) ? i : -1;
      }
    }
    li = -1;
    return -1;
#endif
#if 0
    // Cache unroll
    int& li = _lastIndex.localData();
    int i = (li > start && li < end && _index[li] <= ind) ? li : start;
    while(i+3 < end)
    {
      if(_index[i + 0] >= ind) { li = i + 0; return (_index[i + 0] == ind) ? i + 0 : -1; }
      if(_index[i + 1] >= ind) { li = i + 1; return (_index[i + 1] == ind) ? i + 1 : -1; }
      if(_index[i + 2] >= ind) { li = i + 2; return (_index[i + 2] == ind) ? i + 2 : -1; }
      if(_index[i + 3] >= ind) { li = i + 3; return (_index[i + 3] == ind) ? i + 3 : -1; }
      i += 4;
    }
    for(; i < end; ++i)
    {
      if(_index[i] >= ind)
      {
        li = i;
        return (_index[i] == ind) ? i : -1;
      }
    }
    li = -1;
    return -1;
#endif
#if 0
    // Bin search
    IndexType* pos = std::lower_bound(_index+start, _index+end, ind);
    return (pos < _index+end && *pos == ind) ? pos - _index : -1;
#endif
#if 0
    // Unroll
    int i = start;
    while(i+3 < end)
    {
      if(_index[i + 0] >= ind) return (_index[i + 0] == ind) ? i + 0 : -1;
      if(_index[i + 1] >= ind) return (_index[i + 1] == ind) ? i + 1 : -1;
      if(_index[i + 2] >= ind) return (_index[i + 2] == ind) ? i + 2 : -1;
      if(_index[i + 3] >= ind) return (_index[i + 3] == ind) ? i + 3 : -1;
      i += 4;
    }
    for(; i < end; ++i)
    {
      if(_index[i] >= ind)
        return (_index[i] == ind) ? i : -1;
    }
    return -1;
#endif
#if 0
    // Count
    int cnt = start;
    for(int i = start; i<end; i++)
      cnt += (_index[i] < ind);
    return (cnt < end && _index[cnt] == ind) ? cnt : -1;
#endif
  }
};



typedef GmCSxSparseMatrixLayout<int, GM_CSR_SPARSE_FORMAT> GmCSRSparseMatrixLayout;

typedef GmCSxSparseMatrixLayout<arma::uword, GM_CSC_SPARSE_FORMAT> GmCSCArmadilloSparseMatrixLayout;

// Currently, the <int, CSC> and <uword, CSR> configurations are used only during unit tests
#ifdef ENABLE_TESTS
typedef GmCSxSparseMatrixLayout<arma::uword, GM_CSR_SPARSE_FORMAT> GmCSRArmadilloSparseMatrixLayout;
typedef GmCSxSparseMatrixLayout<int,         GM_CSC_SPARSE_FORMAT> GmCSCSparseMatrixLayout;
#endif


#endif