sparse_matrix.h

Go to the documentation of this file.

#ifndef SPARSE_MATRIX_H
#define SPARSE_MATRIX_H

#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdlib>

#include <vector>

#include <map>

#include <iomanip>

#include <starpu.h>

#include "skepu/src/debug.h"

#include "skepu/src/malloc_allocator.h"





namespace skepu
{
  
 enum SparseFileFormat
 {
  MATRIX_MARKET_FORMAT,
  MATLAB_FORMAT,
  RUTHERFOR_BOEING_FORMAT
 }; 
 

template<typename T>
class SparseMatrix
{

public:
    T *m_values;
    unsigned int *m_colInd;
    unsigned int *m_rowPtr; 
    unsigned int m_rows;
    unsigned int m_cols;
    unsigned int m_nnz;
    T m_zeroValue;
    bool m_dealloc;



// ***** RELATED to transpose of the matrix ****/   
    T *m_valuesCSC;
    unsigned int *m_colPtrCSC;
    unsigned int *m_rowIndCSC; 
    bool m_transposeValid;
    SparseMatrix *m_cscMatrix;
    
    // 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;


public:
    class iterator; 
    
private: //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_ABORT();
        }
            
        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 unsigned int[m_nnz];
        if( m_colPtrCSC == NULL)
            m_colPtrCSC = new unsigned int[m_cols+1]; 
            
        int rowIdx = 0;
        int nxtRowIdx = 0;
        
        std::multimap<unsigned int, std::pair<unsigned int, T>, std::less<unsigned int> > cscFormat;
           
        for(int ii = 0; ii < m_rows; ii++)
        {
            rowIdx = m_rowPtr[ii];
            nxtRowIdx = m_rowPtr[ii+1];
    
            for(int 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<unsigned int, std::pair<unsigned int, T>, std::less<unsigned int> >::iterator mapIter;
        mapIter m_it, s_it;
        
        int col = 0, ind=0;
    
        m_colPtrCSC[0] = 0;
        
        unsigned int expColInd = 0; // assign it to first column index
        
        for (m_it = cscFormat.begin();  m_it != cscFormat.end();  m_it = s_it)
        {
            unsigned int colInd = (*m_it).first;
            
            if(expColInd < colInd) // meaning some intermediate columns are totally blank, i.e., no non-zero entries, add their record
            {
                for(int 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(int i= (col+1); i<=m_cols; i++)
        {
            m_colPtrCSC[i] = ind;
        }
        
        m_transposeValid = true;
    }
    
    
    void readMTXFile(const std::string &inputfile);
    
public:

    void printMatrixInDenseFormat()
    {
        std::cerr<<"SparseMatrix ("<<m_rows <<" X "<< m_cols<<") nnz: "<<m_nnz<<"\n";
        
        T *temp = new T[m_cols];
        
        int rowIdx, nxtRowIdx;
        for(int ii = 0; ii < m_rows; ii++)
        {
            for(int i=0; i<m_cols; i++)
                temp[i] = T();
            
            rowIdx = m_rowPtr[ii];
            nxtRowIdx = m_rowPtr[ii+1];
    
            for(int jj=rowIdx; jj<nxtRowIdx; jj++)
            {
                temp[m_colInd[jj]] = m_values[jj];
            }
            for(int 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(unsigned int rows, unsigned int cols, unsigned int nnz, T *values, unsigned int *rowPtr, unsigned int *colInd, bool dealloc=true, T zeroValue=T(), bool transMatrix=false);
    
    SparseMatrix(unsigned int rows, unsigned int cols, unsigned int 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();
    

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

    friend std::ostream& operator<<(std::ostream &os, SparseMatrix<T>& matrix)
    {
        matrix.acquireRead();
        
        os<<"Matrix rows="<< matrix.total_rows() <<", cols="<<matrix.total_cols()<<", nnz="<<matrix.total_nnz()<<"\n";
        
        int rowIdx = 0;
        int nxtRowIdx = 0;
        
        for(int ii = 0; ii < matrix.total_rows(); ii++)
        {
            rowIdx = matrix.m_rowPtr[ii];
            nxtRowIdx = matrix.m_rowPtr[ii+1];

            for(int 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: 

    unsigned int total_nnz() const { return m_nnz;}
    
    unsigned int total_rows() const {return m_rows; }
    
    unsigned int total_cols() const {return m_cols; }
    
    T* get_values() { /* unregisterSparseMatrix(); */ acquireReadWrite(); return m_values; }
    
    unsigned int* get_row_pointers() const {return m_rowPtr; }
    
    unsigned int* get_col_indices() const {return m_colInd; }
    
    int get_rowSize(unsigned int row) const 
    {
        return (m_rowPtr[row+1]-m_rowPtr[row]);
    }
    
    int get_rowOffsetFromStart(unsigned int 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(unsigned int row, unsigned int col ) const;
    const T& at(unsigned int index) const;
    
    // Similar to 'at' operator but don't do memory check etc. must faster used internally in implementation.
    const T& at_internal(unsigned int row, unsigned int col ) const;
    T& at_internal(unsigned int row, unsigned int col );
    const T& at_internal(unsigned int index) const;
    
    // Element access using () operator
    const T& operator()(const unsigned int row, const unsigned int col) const;
    T operator()(const unsigned int row, const unsigned int 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();


    void acquireRead() const;
    void acquireReadWrite();
        
private:

    void unpartitionMatrix() const
    {
        DEBUG_TEXT_LEVEL1("*****  SPARSE MATRIX UNPARTIOTIINING parts: "<< parts <<" *****\n")
        
        if(parts>1)
            starpu_data_unpartition(csr_matrix_handle, 0); // need to check as it might only unpartition only one filter.

        parts=1;
    }
    
    void release_acquire()
    {
        if(isAcquire)
        {
            if(parts>1) // partitioning
            {
                DEBUG_TEXT_LEVEL1("*****  SPARSE MATRIX RELEASE **** parts: "<< parts <<"\n")
                
                // parts should be equal to what parts variable points to
                assert (parts == starpu_data_get_nb_children(csr_matrix_handle));
                for (int x = 0; x < parts; x++) 
                {
                    starpu_data_release(starpu_data_get_sub_data(csr_matrix_handle, 1, x));
                }
            }
            else // no partitioning
            {
                DEBUG_TEXT_LEVEL1("*****  SPARSE MATRIX RELEASE ****\n")
                starpu_data_release(csr_matrix_handle);
            }
            isAcquire = false;
        }
    }
    

public:

    starpu_data_handle_t csr_matrix_handle;
    mutable bool isOnStarPU;
    mutable bool isReadBack;
    mutable bool isWriteBack;
    mutable int parts; // partitioning
    
    mutable bool isAcquire;
    mutable enum starpu_data_access_mode mode;
    
    struct starpu_data_filter csr_matrix_filter;
    
    
    starpu_data_handle_t registerSparseMatrix()
    {
        release_acquire();
        
        if(!isOnStarPU)
        {
            starpu_csr_data_register(&csr_matrix_handle,0, m_nnz, m_rows, (uintptr_t)m_values, (uint32_t *)m_colInd, (uint32_t *)m_rowPtr, 0, sizeof(m_values[0]));
            isOnStarPU=true;
            DEBUG_TEXT_LEVEL1("*****  SPARSE MATRIX REGISTERING **** rows: "<< m_rows <<" cols: "<< m_cols <<" \n")
        }
        else if(parts>1)
        {
            unpartitionMatrix();
        }
        return csr_matrix_handle;
    }
    
    void unregisterSparseMatrix(bool update=true)
    {
        if(isOnStarPU)
        {
            release_acquire(); // Should we release before unregister? Yes!
                
            if(parts>1) // do we need to unpartition before unregistering? Yes!
            {
                unpartitionMatrix();
            }
            
            if(update)  
                starpu_data_unregister(csr_matrix_handle);
                
            else // dont need to get updated data back
                starpu_data_unregister_no_coherency(csr_matrix_handle);
                
            isOnStarPU = false;
            isAcquire = false;      
        }
    }
    
    starpu_data_handle_t registerPartitions(int _parts=1)
    {
        release_acquire();
        
        if(!isOnStarPU) // if not registered, register it first
        {
            starpu_csr_data_register(&csr_matrix_handle,0, m_nnz, m_rows, (uintptr_t)m_values, (uint32_t *)m_colInd, (uint32_t *)m_rowPtr, 0, sizeof(m_values[0]));
            isOnStarPU=true;
            DEBUG_TEXT_LEVEL1("*****  SPARSE MATRIX REGISTERING **** "<< m_nnz <<"\n")
            parts=1;
        }
        
        if(_parts == parts)
            return csr_matrix_handle;
                
        if(parts>1)
        {
            starpu_data_unpartition(csr_matrix_handle, 0); // first unpartition previous partitioning.
            DEBUG_TEXT_LEVEL1("*****  SPARSE MATRIX UNPARTITIONING **** parts: "<< parts <<"\n")
            parts = 1;
        }
        if(_parts>1)
        {   
            csr_matrix_filter.nchildren = _parts;
            starpu_data_partition(csr_matrix_handle, &csr_matrix_filter);
            
            parts = _parts;
            
            DEBUG_TEXT_LEVEL1("*****  SPARSE MATRIX PARTITIONING **** parts: "<< _parts <<"\n")
        }
    
        return csr_matrix_handle;
    }
};

}

#include "src/sparse_matrix_iterator.inl"

#include "src/sparse_matrix.inl"

#endif