sparse_matrix.h

Go to the documentation of this file.

#ifndef _SPARSE_MATRIX_H
#define _SPARSE_MATRIX_H

#include <cstdlib>
#include <map>
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
#include <map>
#include <iomanip>
#include <set>

#ifdef SKEPU_OPENCL
#ifdef USE_MAC_OPENCL
#include <OpenCL/opencl.h>
#else
#include <CL/cl.h>
#endif
#include "src/device_mem_pointer_cl.h"
#endif

#ifdef SKEPU_CUDA
#include "src/device_mem_pointer_cu.h"
#endif

#include "skepu/src/malloc_allocator.h"
#include "skepu/src/environment.h"


namespace skepu
{

enum SparseFileFormat
{
   MATRIX_MARKET_FORMAT,
   MATLAB_FORMAT,
   RUTHERFOR_BOEING_FORMAT
};


template <class T>
class SparseMatrix
{

// typedefs
public:
#ifdef SKEPU_CUDA
   typedef DeviceMemPointer_CU<T>* device_pointer_type_cu;
   typedef DeviceMemPointer_CU<size_t>* device_pointer_index_type_cu;
#endif

#ifdef SKEPU_OPENCL
   typedef DeviceMemPointer_CL<T>* device_pointer_type_cl;
   typedef DeviceMemPointer_CL<size_t>* device_pointer_index_type_cl;
#endif

private:
   T *m_values;
   size_t *m_colInd;
   size_t *m_rowPtr;
   size_t m_rows;
   size_t m_cols;
   size_t m_nnz;
   T m_zeroValue;
   bool m_dealloc;

//    bool m_valid; /*! to keep track of whether the main copy is valid or not */

   // ***** RELATED to transpose of the matrix ****/
   T *m_valuesCSC;
   size_t *m_colPtrCSC;
   size_t *m_rowIndCSC;
   bool m_transposeValid;
   SparseMatrix *m_cscMatrix;

   

public:
   class iterator;
   
   // following bool to is to control memory deallocation as we assign de-allocation responsibility to main matrix and not the transpose matrix itself.
   bool m_transMatrix;

private:

   void deleteCSCFormat() // when resizing so we need to take transpose again and also would have to allocate storgae again.
   {
      if(m_transMatrix) // do nothing then
         return;

      m_transposeValid = false;

      if(m_colPtrCSC!=NULL)
      {
         delete m_colPtrCSC;
         m_colPtrCSC = NULL;
      }
      if(m_rowIndCSC!=NULL)
      {
         delete m_rowIndCSC;
         m_rowIndCSC = NULL;
      }
      if(m_valuesCSC!=NULL)
      {
         delete m_valuesCSC;
         m_valuesCSC = NULL;
      }

      if(m_cscMatrix!=NULL)
      {
         delete m_cscMatrix;
         m_cscMatrix = NULL;
      }
   }

   void convertToCSCFormat()
   {
      if(m_transMatrix) // cannot transpose an already transposed
      {
         std::cerr<<"Cannot apply transpose operation on a matrix which is already in CSC format.\n";
         SKEPU_EXIT();
      }

      if(m_transposeValid && (m_valuesCSC!=NULL) && (m_rowIndCSC!=NULL) && (m_colPtrCSC!=NULL)) // already there
         return;

      updateHost(); // update Host for new data;

      if(m_valuesCSC==NULL)
         m_valuesCSC = new T[m_nnz];
      if(m_rowIndCSC == NULL)
         m_rowIndCSC = new size_t[m_nnz];
      if(m_colPtrCSC == NULL)
         m_colPtrCSC = new size_t[m_cols+1];

      size_t rowIdx = 0;
      size_t nxtRowIdx = 0;

      std::multimap<size_t, std::pair<size_t, T>, std::less<size_t> > cscFormat;

      for(size_t ii = 0; ii < m_rows; ii++)
      {
         rowIdx = m_rowPtr[ii];
         nxtRowIdx = m_rowPtr[ii+1];

         for(size_t jj=rowIdx; jj<nxtRowIdx; jj++)
         {
            cscFormat.insert( std::make_pair(m_colInd[jj], std::make_pair(ii, m_values[jj])) );
         }
      }

      // make multimap sorted based on key which is column index.
      typedef typename std::multimap<size_t, std::pair<size_t, T>, std::less<size_t> >::iterator mapIter;
      mapIter m_it, s_it;

      size_t col = 0, ind=0;

      m_colPtrCSC[0] = 0;

      size_t expColInd = 0; // assign it to first column index

      for (m_it = cscFormat.begin();  m_it != cscFormat.end();  m_it = s_it)
      {
         size_t colInd = (*m_it).first;

         if(expColInd < colInd) // meaning some intermediate columns are totally blank, i.e., no non-zero entries, add their record
         {
            for(size_t i=expColInd; i<colInd; i++)
            {
               m_colPtrCSC[++col] = ind;
            }
         }
         expColInd = colInd+1; // now set expected to next column.

         std::pair<mapIter, mapIter> keyRange = cscFormat.equal_range(colInd);

         // Iterate over all map elements with key == theKey
         for (s_it = keyRange.first;  s_it != keyRange.second;  ++s_it)
         {
            m_rowIndCSC[ind] = (*s_it).second.first;
            m_valuesCSC[ind] = (*s_it).second.second;
            ind++;
         }
         m_colPtrCSC[++col] = ind;
      }
      for(size_t i= (col+1); i<=m_cols; i++)
      {
         m_colPtrCSC[i] = ind;
      }

      m_transposeValid = true;
   }


   void readMTXFile(const std::string &inputfile);

public:
   
// #if SKEPU_DEBUG>0      
   std::string m_nameVerbose; // for debugging useful
// #endif

   void printMatrixInDenseFormat()
   {
      std::cerr<<"SparseMatrix ("<<m_rows <<" X "<< m_cols<<") nnz: "<<m_nnz<<"\n";

      T *temp = new T[m_cols];

      size_t rowIdx, nxtRowIdx;
      for(size_t ii = 0; ii < m_rows; ii++)
      {
         for(size_t i=0; i<m_cols; i++)
            temp[i] = T();

         rowIdx = m_rowPtr[ii];
         nxtRowIdx = m_rowPtr[ii+1];

         for(size_t jj=rowIdx; jj<nxtRowIdx; jj++)
         {
            temp[m_colInd[jj]] = m_values[jj];
         }

         for(size_t i=0; i<m_cols; i++)
            std::cerr<<std::setw(5)<<temp[i];

         std::cerr<<"\n";
      }

      delete[] temp;
   }

   // unary CSC operator
   inline SparseMatrix<T>& operator~()
   {
      convertToCSCFormat();

      if(m_cscMatrix==NULL)
      {
         m_cscMatrix = new SparseMatrix(m_cols, m_rows, m_nnz, m_valuesCSC, m_colPtrCSC, m_rowIndCSC, false, T(), true);
      }

      return *m_cscMatrix;
   }

   // ----- CONSTRUCTORS & DESTRUCTORS -------//
   SparseMatrix(size_t rows, size_t cols, size_t nnz, T *values, size_t *rowPtr, size_t *colInd, bool dealloc=true, T zeroValue=T(), bool transMatrix=false);

   SparseMatrix(size_t rows, size_t cols, size_t nnz, T min, T max, T zeroValue=T());

   SparseMatrix(const SparseMatrix &copy);

   SparseMatrix(const std::string &inputfile, enum SparseFileFormat format=MATRIX_MARKET_FORMAT, T zeroValue=T());

   ~SparseMatrix();
   
   void operator=(const SparseMatrix<T> &copy);

   // ----- FRIEND METHODS       -------//

   friend std::ostream& operator<<(std::ostream &os, SparseMatrix<T>& matrix)
   {
      matrix.updateHost();

      os<<"Matrix rows="<< matrix.total_rows() <<", cols="<<matrix.total_cols()<<", nnz="<<matrix.total_nnz()<<"\n";

      size_t rowIdx = 0;
      size_t nxtRowIdx = 0;

      for(size_t ii = 0; ii < matrix.total_rows(); ii++)
      {
         rowIdx = matrix.m_rowPtr[ii];
         nxtRowIdx = matrix.m_rowPtr[ii+1];

         for(size_t jj=rowIdx; jj<nxtRowIdx; jj++)
         {
            os<< "row: "<<std::setw(8)<<ii<<", col: "<<std::setw(8)<<matrix.m_colInd[jj]<<", value"<<std::setw(12)<<matrix.m_values[jj]<<"\n";
         }
      }

      os<<"\n";
      return os;
   }


public:

   size_t total_nnz() const
   {
      return m_nnz;
   }

   size_t total_rows() const
   {
      return m_rows;
   }

   size_t total_cols() const
   {
      return m_cols;
   }

   T* get_values() const
   {
      return m_values;
   }

   size_t* get_row_pointers() const
   {
      return m_rowPtr;
   }

   size_t* get_col_indices() const
   {
      return m_colInd;
   }

   size_t get_rowSize(size_t row) const
   {
      return (m_rowPtr[row+1]-m_rowPtr[row]);
   }

   size_t get_rowOffsetFromStart(size_t row) const
   {
      if(row<total_rows())
         return m_rowPtr[row];

      return m_nnz;
   }

public:

   iterator begin(unsigned row);

   // Element access using 'at' operator
   T at(size_t row, size_t col ) const;
   const T& at(size_t index) const;

   // Similar to 'at' operator but don't do memory check etc. must faster used internally in implementation.
   const T& at_internal(size_t row, size_t col ) const;
   T& at_internal(size_t row, size_t col );
   const T& at_internal(size_t index) const;

   // Element access using () operator
   const T& operator()(const size_t row, const size_t col) const;
   T operator()(const size_t row, const size_t col);

   // will resize the matrix, can be dangerous, used with care
   void resize(SparseMatrix<T> &copy, bool retainData);

   void updateHost() const;
   void invalidateDevice();
   void updateHostAndInvalidateDevice();

#ifdef SKEPU_OPENCL
   device_pointer_type_cl updateDevice_CL(T* start, size_t elems, Device_CL* device, bool copy);
   device_pointer_index_type_cl updateDevice_Index_CL(size_t* start, size_t elems, Device_CL* device, bool copy);
   void flush_CL();
#endif

#ifdef SKEPU_CUDA
   device_pointer_type_cu updateDevice_CU(T* start, size_t elems, unsigned int deviceID, bool copy);
   device_pointer_index_type_cu updateDevice_Index_CU(size_t* start, size_t elems, unsigned int deviceID, bool copy);
   void flush_CU();
   bool isSparseMatrixOnDevice_CU(unsigned int deviceID);
   bool isModified_CU(unsigned int deviceID);
#endif
  

private:

#ifdef SKEPU_OPENCL
   void updateHost_CL() const;
   void invalidateDeviceData_CL();
   void releaseDeviceAllocations_CL();
#endif

#ifdef SKEPU_CUDA
   void updateHost_CU() const;
   void invalidateDeviceData_CU();
   void releaseDeviceAllocations_CU();
#endif

#ifdef SKEPU_CUDA
   std::map<std::pair< unsigned int, std::pair< T*, size_t > >, device_pointer_type_cu > m_deviceMemPointers_CU;
   std::map<std::pair< unsigned int, std::pair< size_t*, size_t > >, device_pointer_index_type_cu > m_deviceMemIndexPointers_CU;
#endif

#ifdef SKEPU_OPENCL
   std::map<std::pair< cl_device_id, std::pair< T*, size_t > >, device_pointer_type_cl > m_deviceMemPointers_CL;
   std::map<std::pair< cl_device_id, std::pair< size_t*, size_t > >, device_pointer_index_type_cl > m_deviceMemIndexPointers_CL;
#endif


}; // end class SparseMatrix


} // end namespace skepu


#include "src/sparse_matrix_iterator.inl"


#include "src/sparse_matrix.inl"

#ifdef SKEPU_OPENCL
#include "src/sparse_matrix_cl.inl"
#endif

#ifdef SKEPU_CUDA
#include "src/sparse_matrix_cu.inl"
#endif

#endif  /* _MATRIX_H */