SpDCCols.cpp

Go to the documentation of this file.

/****************************************************************/
/* Parallel Combinatorial BLAS Library (for Graph Computations) */
/* version 1.6 -------------------------------------------------*/
/* date: 6/15/2017 ---------------------------------------------*/
/* authors: Ariful Azad, Aydin Buluc  --------------------------*/
/****************************************************************/
/*
 Copyright (c) 2010-2017, The Regents of the University of California
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
 in the Software without restriction, including without limitation the rights
 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 copies of the Software, and to permit persons to whom the Software is
 furnished to do so, subject to the following conditions:
 
 The above copyright notice and this permission notice shall be included in
 all copies or substantial portions of the Software.
 
 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 */


#include "SpDCCols.h"
#include "Deleter.h"
#include <algorithm>
#include <functional>
#include <vector>
#include <climits>
#include <iomanip>
#include <cassert>

namespace combblas {

/****************************************************************************/
/********************* PUBLIC CONSTRUCTORS/DESTRUCTORS **********************/
/****************************************************************************/

template <class IT, class NT>
const IT SpDCCols<IT,NT>::esscount = static_cast<IT>(4);


template <class IT, class NT>
SpDCCols<IT,NT>::SpDCCols():dcsc(NULL), m(0), n(0), nnz(0), splits(0){
}

// Allocate all the space necessary
template <class IT, class NT>
SpDCCols<IT,NT>::SpDCCols(IT size, IT nRow, IT nCol, IT nzc)
:m(nRow), n(nCol), nnz(size), splits(0)
{
    if(nnz > 0)
        dcsc = new Dcsc<IT,NT>(nnz, nzc);
    else
        dcsc = NULL; 
}

template <class IT, class NT>
SpDCCols<IT,NT>::~SpDCCols()
{
    if(nnz > 0)
    {
        if(dcsc != NULL) 
        {   
            if(splits > 0)
            {
                for(int i=0; i<splits; ++i)
                    delete dcscarr[i];
                delete [] dcscarr;
            }
            else
            {
                delete dcsc;    
            }
        }
    }
}

// Copy constructor (constructs a new object. i.e. this is NEVER called on an existing object)
// Derived's copy constructor can safely call Base's default constructor as base has no data members 
template <class IT, class NT>
SpDCCols<IT,NT>::SpDCCols(const SpDCCols<IT,NT> & rhs)
: m(rhs.m), n(rhs.n), nnz(rhs.nnz), splits(rhs.splits)
{
    if(splits > 0)
    {
        for(int i=0; i<splits; ++i)
            CopyDcsc(rhs.dcscarr[i]);
    }
    else
    {
        CopyDcsc(rhs.dcsc);
    }
}

template <class IT, class NT>
SpDCCols<IT,NT>::SpDCCols(const SpTuples<IT, NT> & rhs, bool transpose)
: m(rhs.m), n(rhs.n), nnz(rhs.nnz), splits(0)
{    
    
    if(nnz == 0)    // m by n matrix of complete zeros
    {
        if(transpose) std::swap(m,n);
        dcsc = NULL;    
    } 
    else
    {
        if(transpose)
        {
      std::swap(m,n);
            IT localnzc = 1;
            for(IT i=1; i< rhs.nnz; ++i)
            {
                if(rhs.rowindex(i) != rhs.rowindex(i-1))
                {
                    ++localnzc;
                }
            }
            dcsc = new Dcsc<IT,NT>(rhs.nnz,localnzc);   
            dcsc->jc[0]  = rhs.rowindex(0); 
            dcsc->cp[0] = 0;

            for(IT i=0; i<rhs.nnz; ++i)
            {
                dcsc->ir[i]  = rhs.colindex(i);     // copy rhs.jc to ir since this transpose=true
                dcsc->numx[i] = rhs.numvalue(i);
            }

            IT jspos  = 1;      
            for(IT i=1; i<rhs.nnz; ++i)
            {
                if(rhs.rowindex(i) != dcsc->jc[jspos-1])
                {
                    dcsc->jc[jspos] = rhs.rowindex(i);  // copy rhs.ir to jc since this transpose=true
                    dcsc->cp[jspos++] = i;
                }
            }       
            dcsc->cp[jspos] = rhs.nnz;
        }
        else
        {
            IT localnzc = 1;
            for(IT i=1; i<rhs.nnz; ++i)
            {
                if(rhs.colindex(i) != rhs.colindex(i-1))
                {
                    ++localnzc;
                }
            }
            dcsc = new Dcsc<IT,NT>(rhs.nnz,localnzc);   
            dcsc->jc[0]  = rhs.colindex(0); 
            dcsc->cp[0] = 0;

            for(IT i=0; i<rhs.nnz; ++i)
            {
                dcsc->ir[i]  = rhs.rowindex(i);     // copy rhs.ir to ir since this transpose=false
                dcsc->numx[i] = rhs.numvalue(i);
            }

            IT jspos = 1;       
            for(IT i=1; i<rhs.nnz; ++i)
            {
                if(rhs.colindex(i) != dcsc->jc[jspos-1])
                {
                    dcsc->jc[jspos] = rhs.colindex(i);  // copy rhs.jc to jc since this transpose=true
                    dcsc->cp[jspos++] = i;
                }
            }       
            dcsc->cp[jspos] = rhs.nnz;
        }
    } 
}




template <class IT, class NT>
SpDCCols<IT,NT>::SpDCCols(IT nRow, IT nCol, IT nTuples, const std::tuple<IT, IT, NT>*  tuples, bool transpose)
: m(nRow), n(nCol), nnz(nTuples), splits(0)
{
    
    if(nnz == 0)    // m by n matrix of complete zeros
    {
        dcsc = NULL;
    }
    else
    {
        int totThreads=1;
#ifdef _OPENMP
#pragma omp parallel
        {
            totThreads = omp_get_num_threads();
        }
#endif
        
        std::vector <IT> tstart(totThreads);
        std::vector <IT> tend(totThreads);
        std::vector <IT> tdisp(totThreads+1);
        
        // extra memory, but replaces an O(nnz) loop by an O(nzc) loop
        IT* temp_jc = new IT[nTuples];
        IT* temp_cp = new IT[nTuples];

#ifdef _OPENMP
#pragma omp parallel
#endif
        {
            int threadID = 0;
#ifdef _OPENMP
            threadID = omp_get_thread_num();
#endif
            IT start = threadID * (nTuples / totThreads);
            IT end = (threadID + 1) * (nTuples / totThreads);
            if(threadID == (totThreads-1)) end = nTuples;
            IT curpos=start;
            if(end>start) // no work for the current thread
            {
                temp_jc[start] = std::get<1>(tuples[start]);
                temp_cp[start] = start;
                for (IT i = start+1; i < end; ++i)
                {
                    if(std::get<1>(tuples[i]) != temp_jc[curpos] )
                    {
                        temp_jc[++curpos] = std::get<1>(tuples[i]);
                        temp_cp[curpos] = i;
                    }
                }
            }
           
            tstart[threadID] = start;
            if(end>start) tend[threadID] = curpos+1;
            else tend[threadID] = end; // start=end
        }

        
        // serial part
        for(int t=totThreads-1; t>0; --t)
        {
            if(tend[t] > tstart[t] && tend[t-1] > tstart[t-1])
            {
                if(temp_jc[tstart[t]] == temp_jc[tend[t-1]-1])
                {
                    tstart[t] ++;
                }
            }
        }
        
        tdisp[0] = 0;
        for(int t=0; t<totThreads; ++t)
        {
            tdisp[t+1] = tdisp[t] + tend[t] - tstart[t];
        }

        IT localnzc = tdisp[totThreads];
        dcsc = new Dcsc<IT,NT>(nTuples,localnzc);
    
#ifdef _OPENMP
#pragma omp parallel
#endif
        {
            int threadID = 0;
#ifdef _OPENMP
            threadID = omp_get_thread_num();
#endif
            std::copy(temp_jc + tstart[threadID],  temp_jc + tend[threadID], dcsc->jc + tdisp[threadID]);
            std::copy(temp_cp + tstart[threadID],  temp_cp + tend[threadID], dcsc->cp + tdisp[threadID]);
        }
        dcsc->cp[localnzc] = nTuples;

        delete [] temp_jc;
        delete [] temp_cp;
        
#ifdef _OPENMP
#pragma omp parallel for schedule (static)
#endif
        for(IT i=0; i<nTuples; ++i)
        {
            dcsc->ir[i]  = std::get<0>(tuples[i]);
            dcsc->numx[i] = std::get<2>(tuples[i]);
        }
     }
    
    if(transpose) Transpose(); // this is not efficient, think to improve later. We included this parameter anyway to make this constructor different from another constracttor when the fourth argument is passed as 0.
}



/*
template <class IT, class NT>
SpDCCols<IT,NT>::SpDCCols(IT nRow, IT nCol, IT nTuples, const tuple<IT, IT, NT>*  tuples)
: m(nRow), n(nCol), nnz(nTuples), splits(0)
{
    
    if(nnz == 0)    // m by n matrix of complete zeros
    {
        dcsc = NULL;
    }
    else
    {
        IT localnzc = 1;
#pragma omp parallel for schedule (static) default(shared) reduction(+:localnzc)
        for(IT i=1; i<nTuples; ++i) // not scaling well, try my own version
        {
            if(std::get<1>(tuples[i]) != std::get<1>(tuples[i-1]))
            {
                ++localnzc;
            }
        }
        
        dcsc = new Dcsc<IT,NT>(nTuples,localnzc);
        dcsc->jc[0]  = std::get<1>(tuples[0]);
        dcsc->cp[0] = 0;
        
#pragma omp parallel for schedule (static)
        for(IT i=0; i<nTuples; ++i)
        {
            dcsc->ir[i]  = std::get<0>(tuples[i]);
            dcsc->numx[i] = std::get<2>(tuples[i]);
        }
        
        IT jspos = 1;
        for(IT i=1; i<nTuples; ++i) // now this loop
        {
            if(std::get<1>(tuples[i]) != dcsc->jc[jspos-1])
            {
                dcsc->jc[jspos] = std::get<1>(tuples[i]);
                dcsc->cp[jspos++] = i;
            }
        }
        dcsc->cp[jspos] = nTuples;
    }
}
*/

/****************************************************************************/
/************************** PUBLIC OPERATORS ********************************/
/****************************************************************************/

template <class IT, class NT>
SpDCCols<IT,NT> & SpDCCols<IT,NT>::operator=(const SpDCCols<IT,NT> & rhs)
{
    // this pointer stores the address of the class instance
    // check for self assignment using address comparison
    if(this != &rhs)        
    {
        if(dcsc != NULL && nnz > 0)
        {
            delete dcsc;
        }
        if(rhs.dcsc != NULL)    
        {
            dcsc = new Dcsc<IT,NT>(*(rhs.dcsc));
            nnz = rhs.nnz;
        }
        else
        {
            dcsc = NULL;
            nnz = 0;
        }
        
        m = rhs.m; 
        n = rhs.n;
        splits = rhs.splits;
    }
    return *this;
}

template <class IT, class NT>
SpDCCols<IT,NT> & SpDCCols<IT,NT>::operator+= (const SpDCCols<IT,NT> & rhs)
{
    // this pointer stores the address of the class instance
    // check for self assignment using address comparison
    if(this != &rhs)        
    {
        if(m == rhs.m && n == rhs.n)
        {
            if(rhs.nnz == 0)
            {
                return *this;
            }
            else if(nnz == 0)
            {
                dcsc = new Dcsc<IT,NT>(*(rhs.dcsc));
                nnz = dcsc->nz;
            }
            else
            {
                (*dcsc) += (*(rhs.dcsc));
                nnz = dcsc->nz;
            }       
        }
        else
        {
            std::cout<< "Not addable: " << m  << "!=" << rhs.m << " or " << n << "!=" << rhs.n <<std::endl;     
        }
    }
    else
    {
        std::cout<< "Missing feature (A+A): Use multiply with 2 instead !"<<std::endl;  
    }
    return *this;
}

template <class IT, class NT>
template <typename _UnaryOperation, typename GlobalIT>
SpDCCols<IT,NT>* SpDCCols<IT,NT>::PruneI(_UnaryOperation __unary_op, bool inPlace, GlobalIT rowOffset, GlobalIT colOffset)
{
    if(nnz > 0)
    {
        Dcsc<IT,NT>* ret = dcsc->PruneI (__unary_op, inPlace, rowOffset, colOffset);
        if (inPlace)
        {
            nnz = dcsc->nz;
    
            if(nnz == 0) 
            {   
                delete dcsc;
                dcsc = NULL;
            }
            return NULL;
        }
        else
        {
            // wrap the new pruned Dcsc into a new SpDCCols
            SpDCCols<IT,NT>* retcols = new SpDCCols<IT, NT>();
            retcols->dcsc = ret;
            retcols->nnz = retcols->dcsc->nz;
            retcols->n = n;
            retcols->m = m;
            return retcols;
        }
    }
    else
    {
        if (inPlace)
        {
            return NULL;
        }
        else
        {
            SpDCCols<IT,NT>* retcols = new SpDCCols<IT, NT>();
            retcols->dcsc = NULL;
            retcols->nnz = 0;
            retcols->n = n;
            retcols->m = m;
            return retcols;
        }
    }
}

template <class IT, class NT>
template <typename _UnaryOperation>
SpDCCols<IT,NT>* SpDCCols<IT,NT>::Prune(_UnaryOperation __unary_op, bool inPlace)
{
    if(nnz > 0)
    {
        Dcsc<IT,NT>* ret = dcsc->Prune (__unary_op, inPlace);
        if (inPlace)
        {
            nnz = dcsc->nz;
    
            if(nnz == 0) 
            {   
                delete dcsc;
                dcsc = NULL;
            }
            return NULL;
        }
        else
        {
            // wrap the new pruned Dcsc into a new SpDCCols
            SpDCCols<IT,NT>* retcols = new SpDCCols<IT, NT>();
            retcols->dcsc = ret;
            retcols->nnz = retcols->dcsc->nz;
            retcols->n = n;
            retcols->m = m;
            return retcols;
        }
    }
    else
    {
        if (inPlace)
        {
            return NULL;
        }
        else
        {
            SpDCCols<IT,NT>* retcols = new SpDCCols<IT, NT>();
            retcols->dcsc = NULL;
            retcols->nnz = 0;
            retcols->n = n;
            retcols->m = m;
            return retcols;
        }
    }
}



template <class IT, class NT>
template <typename _BinaryOperation>
SpDCCols<IT,NT>* SpDCCols<IT,NT>::PruneColumn(NT* pvals, _BinaryOperation __binary_op, bool inPlace)
{
    if(nnz > 0)
    {
        Dcsc<IT,NT>* ret = dcsc->PruneColumn (pvals, __binary_op, inPlace);
        if (inPlace)
        {
            nnz = dcsc->nz;
            
            if(nnz == 0)
            {
                delete dcsc;
                dcsc = NULL;
            }
            return NULL;
        }
        else
        {
            // wrap the new pruned Dcsc into a new SpDCCols
            SpDCCols<IT,NT>* retcols = new SpDCCols<IT, NT>();
            retcols->dcsc = ret;
            retcols->nnz = retcols->dcsc->nz;
            retcols->n = n;
            retcols->m = m;
            return retcols;
        }
    }
    else
    {
        if (inPlace)
        {
            return NULL;
        }
        else
        {
            SpDCCols<IT,NT>* retcols = new SpDCCols<IT, NT>();
            retcols->dcsc = NULL;
            retcols->nnz = 0;
            retcols->n = n;
            retcols->m = m;
            return retcols;
        }
    }
}


template <class IT, class NT>
template <typename _BinaryOperation>
SpDCCols<IT,NT>* SpDCCols<IT,NT>::PruneColumn(IT* pinds, NT* pvals, _BinaryOperation __binary_op, bool inPlace)
{
    if(nnz > 0)
    {
        Dcsc<IT,NT>* ret = dcsc->PruneColumn (pinds, pvals, __binary_op, inPlace);
        if (inPlace)
        {
            nnz = dcsc->nz;
            
            if(nnz == 0)
            {
                delete dcsc;
                dcsc = NULL;
            }
            return NULL;
        }
        else
        {
            // wrap the new pruned Dcsc into a new SpDCCols
            SpDCCols<IT,NT>* retcols = new SpDCCols<IT, NT>();
            retcols->dcsc = ret;
            retcols->nnz = retcols->dcsc->nz;
            retcols->n = n;
            retcols->m = m;
            return retcols;
        }
    }
    else
    {
        if (inPlace)
        {
            return NULL;
        }
        else
        {
            SpDCCols<IT,NT>* retcols = new SpDCCols<IT, NT>();
            retcols->dcsc = NULL;
            retcols->nnz = 0;
            retcols->n = n;
            retcols->m = m;
            return retcols;
        }
    }
}



template <class IT, class NT>
void SpDCCols<IT,NT>::EWiseMult (const SpDCCols<IT,NT> & rhs, bool exclude)
{
    if(this != &rhs)        
    {
        if(m == rhs.m && n == rhs.n)
        {
            if(rhs.nnz == 0)
            {
                if(exclude) // then we don't exclude anything
                {
                    return;
                }
                else    // A .* zeros() is zeros()
                {
                    *this = SpDCCols<IT,NT>(0,m,n,0);   // completely reset the matrix
                }
            }
            else if (rhs.nnz != 0 && nnz != 0)
            {
                dcsc->EWiseMult (*(rhs.dcsc), exclude);
                nnz = dcsc->nz;
                if(nnz == 0 )
                    dcsc = NULL;
            }       
        }
        else
        {
            std::cout<< "Matrices do not conform for A .* op(B) !"<<std::endl;      
        }
    }
    else
    {
        std::cout<< "Missing feature (A .* A): Use Square_EWise() instead !"<<std::endl;    
    }
}

template <class IT, class NT>
void SpDCCols<IT,NT>::EWiseScale(NT ** scaler, IT m_scaler, IT n_scaler)
{
    if(m == m_scaler && n == n_scaler)
    {
        if(nnz > 0)
            dcsc->EWiseScale (scaler);
    }
    else
    {
        std::cout<< "Matrices do not conform for EWiseScale !"<<std::endl;      
    }
}   


/****************************************************************************/
/********************* PUBLIC MEMBER FUNCTIONS ******************************/
/****************************************************************************/

template <class IT, class NT>
void SpDCCols<IT,NT>::CreateImpl(IT * _cp, IT * _jc, IT * _ir, NT * _numx, IT _nz, IT _nzc, IT _m, IT _n)
{
    m = _m;
    n = _n;
    nnz =  _nz;
    
    if(nnz > 0)
        dcsc = new Dcsc<IT,NT>(_cp, _jc, _ir, _numx, _nz, _nzc, false); // memory not owned by DCSC
    else
        dcsc = NULL; 
}

template <class IT, class NT>
void SpDCCols<IT,NT>::CreateImpl(const std::vector<IT> & essentials)
{
    assert(essentials.size() == esscount);
    nnz = essentials[0];
    m = essentials[1];
    n = essentials[2];

    if(nnz > 0)
        dcsc = new Dcsc<IT,NT>(nnz,essentials[3]);
    else
        dcsc = NULL; 
}

template <class IT, class NT>
void SpDCCols<IT,NT>::CreateImpl(IT size, IT nRow, IT nCol, std::tuple<IT, IT, NT> * mytuples)
{
    SpTuples<IT,NT> tuples(size, nRow, nCol, mytuples);        
    tuples.SortColBased();
    
#ifdef DEBUG
  std::pair<IT,IT> rlim = tuples.RowLimits(); 
  std::pair<IT,IT> clim = tuples.ColLimits();

  std::ofstream oput;
  std::stringstream ss;
  std::string rank;
    int myrank;
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    ss << myrank;
    ss >> rank;
  std::string ofilename = "Read";
    ofilename += rank;
    oput.open(ofilename.c_str(), std::ios_base::app );
    oput << "Creating of dimensions " << nRow << "-by-" << nCol << " of size: " << size << 
            " with row range (" << rlim.first  << "," << rlim.second << ") and column range (" << clim.first  << "," << clim.second << ")" << std::endl;
    if(tuples.getnnz() > 0)
    { 
        IT minfr = joker::get<0>(tuples.front());
        IT minto = joker::get<1>(tuples.front());
        IT maxfr = joker::get<0>(tuples.back());
        IT maxto = joker::get<1>(tuples.back());

        oput << "Min: " << minfr << ", " << minto << "; Max: " << maxfr << ", " << maxto << std::endl;
    }
    oput.close();
#endif

    SpDCCols<IT,NT> object(tuples, false);  
    *this = object;
}


template <class IT, class NT>
std::vector<IT> SpDCCols<IT,NT>::GetEssentials() const
{
    std::vector<IT> essentials(esscount);
    essentials[0] = nnz;
    essentials[1] = m;
    essentials[2] = n;
    essentials[3] = (nnz > 0) ? dcsc->nzc : 0;
    return essentials;
}

template <class IT, class NT>
template <typename NNT>
SpDCCols<IT,NT>::operator SpDCCols<IT,NNT> () const
{
    Dcsc<IT,NNT> * convert;
    if(nnz > 0)
        convert = new Dcsc<IT,NNT>(*dcsc);
    else
        convert = NULL;

    return SpDCCols<IT,NNT>(m, n, convert);
}


template <class IT, class NT>
template <typename NIT, typename NNT>
SpDCCols<IT,NT>::operator SpDCCols<NIT,NNT> () const
{
    Dcsc<NIT,NNT> * convert;
    if(nnz > 0)
        convert = new Dcsc<NIT,NNT>(*dcsc);
    else
        convert = NULL;

    return SpDCCols<NIT,NNT>(m, n, convert);
}


template <class IT, class NT>
Arr<IT,NT> SpDCCols<IT,NT>::GetArrays() const
{
    Arr<IT,NT> arr(3,1);

    if(nnz > 0)
    {
        arr.indarrs[0] = LocArr<IT,IT>(dcsc->cp, dcsc->nzc+1);
        arr.indarrs[1] = LocArr<IT,IT>(dcsc->jc, dcsc->nzc);
        arr.indarrs[2] = LocArr<IT,IT>(dcsc->ir, dcsc->nz);
        arr.numarrs[0] = LocArr<NT,IT>(dcsc->numx, dcsc->nz);
    }
    else
    {
        arr.indarrs[0] = LocArr<IT,IT>(NULL, 0);
        arr.indarrs[1] = LocArr<IT,IT>(NULL, 0);
        arr.indarrs[2] = LocArr<IT,IT>(NULL, 0);
        arr.numarrs[0] = LocArr<NT,IT>(NULL, 0);
    
    }
    return arr;
}

template <class IT, class NT>
void SpDCCols<IT,NT>::Transpose()
{
    if(nnz > 0)
    {
        SpTuples<IT,NT> Atuples(*this);
        Atuples.SortRowBased();

        // destruction of (*this) is handled by the assignment operator
        *this = SpDCCols<IT,NT>(Atuples,true);
    }
    else
    {
        *this = SpDCCols<IT,NT>(0, n, m, 0);
    }
}


template <class IT, class NT>
SpDCCols<IT,NT> SpDCCols<IT,NT>::TransposeConst() const
{
    SpTuples<IT,NT> Atuples(*this);
    Atuples.SortRowBased();

    return SpDCCols<IT,NT>(Atuples,true);
}

template <class IT, class NT>
SpDCCols<IT,NT> * SpDCCols<IT,NT>::TransposeConstPtr() const
{
    SpTuples<IT,NT> Atuples(*this);
    Atuples.SortRowBased();
    
    return new SpDCCols<IT,NT>(Atuples,true);
}

template <class IT, class NT>
void SpDCCols<IT,NT>::Split(SpDCCols<IT,NT> & partA, SpDCCols<IT,NT> & partB) 
{
    IT cut = n/2;
    if(cut == 0)
    {
        std::cout<< "Matrix is too small to be splitted" << std::endl;
        return;
    }

    Dcsc<IT,NT> *Adcsc = NULL;
    Dcsc<IT,NT> *Bdcsc = NULL;

    if(nnz != 0)
    {
        dcsc->Split(Adcsc, Bdcsc, cut);
    }

    partA = SpDCCols<IT,NT> (m, cut, Adcsc);
    partB = SpDCCols<IT,NT> (m, n-cut, Bdcsc);
    
    // handle destruction through assignment operator
    *this = SpDCCols<IT, NT>();     
}

template <class IT, class NT>
void SpDCCols<IT,NT>::ColSplit(int parts, std::vector< SpDCCols<IT,NT> > & matrices)
{
    if(parts < 2)
    {
        matrices.emplace_back(*this);
    }
    else
    {
        std::vector<IT> cuts(parts-1);
        for(int i=0; i< (parts-1); ++i)
        {
            cuts[i] = (i+1) * (n/parts);
        }
        if(n < parts)
        {
            std::cout<< "Matrix is too small to be splitted" << std::endl;
            return;
        }
        std::vector< Dcsc<IT,NT> * > dcscs(parts, NULL);
        
        if(nnz != 0)
        {
            dcsc->ColSplit(dcscs, cuts);
        }
        
        for(int i=0; i< (parts-1); ++i)
        {
            SpDCCols<IT,NT> matrix = SpDCCols<IT,NT>(m, (n/parts), dcscs[i]);
            matrices.emplace_back(matrix);
        }
        SpDCCols<IT,NT> matrix = SpDCCols<IT,NT>(m, n-cuts[parts-2], dcscs[parts-1]);
        matrices.emplace_back(matrix);
    }
    *this = SpDCCols<IT, NT>();         // handle destruction through assignment operator
}


template <class IT, class NT>
void SpDCCols<IT,NT>::ColConcatenate(std::vector< SpDCCols<IT,NT> > & matrices)
{
    std::vector< SpDCCols<IT,NT> * > nonempties;
    std::vector< Dcsc<IT,NT> * > dcscs;
    std::vector< IT > offsets;
    IT runningoffset = 0;

    for(size_t i=0; i< matrices.size(); ++i)
    {
        if(matrices[i].nnz != 0)
        {
            nonempties.push_back(&(matrices[i]));
            dcscs.push_back(matrices[i].dcsc);
            offsets.push_back(runningoffset);
        }
        runningoffset += matrices[i].n;
    }
    
    if(nonempties.size() < 1)
    {
#ifdef DEBUG
        std::cout << "Nothing to ColConcatenate" << std::endl;
#endif
        n = runningoffset;
    }
    else if(nonempties.size() < 2)
    {
        *this =  *(nonempties[0]);
        n = runningoffset; 
    }
    else // nonempties.size() > 1
    {
        Dcsc<IT,NT> * Cdcsc = new Dcsc<IT,NT>();
        Cdcsc->ColConcatenate(dcscs, offsets);
        *this = SpDCCols<IT,NT> (nonempties[0]->m, runningoffset, Cdcsc);
    }
    
    // destruct parameters
    for(size_t i=0; i< matrices.size(); ++i)
    {
        matrices[i] = SpDCCols<IT,NT>();
    }
}



template <class IT, class NT>
void SpDCCols<IT,NT>::Merge(SpDCCols<IT,NT> & partA, SpDCCols<IT,NT> & partB) 
{
    assert( partA.m == partB.m );

    Dcsc<IT,NT> * Cdcsc = new Dcsc<IT,NT>();

    if(partA.nnz == 0 && partB.nnz == 0)
    {
        Cdcsc = NULL;
    }
    else if(partA.nnz == 0)
    {
        Cdcsc = new Dcsc<IT,NT>(*(partB.dcsc));
    std::transform(Cdcsc->jc, Cdcsc->jc + Cdcsc->nzc, Cdcsc->jc, std::bind2nd(std::plus<IT>(), partA.n));
    }
    else if(partB.nnz == 0)
    {
        Cdcsc = new Dcsc<IT,NT>(*(partA.dcsc));
    }
    else
    {
        Cdcsc->Merge(partA.dcsc, partB.dcsc, partA.n);
    }
    *this = SpDCCols<IT,NT> (partA.m, partA.n + partB.n, Cdcsc);

    partA = SpDCCols<IT, NT>(); 
    partB = SpDCCols<IT, NT>();
}

template <class IT, class NT>
template <class SR>
int SpDCCols<IT,NT>::PlusEq_AnXBt(const SpDCCols<IT,NT> & A, const SpDCCols<IT,NT> & B)
{
    if(A.isZero() || B.isZero())
    {
        return -1;  // no need to do anything
    }
    Isect<IT> *isect1, *isect2, *itr1, *itr2, *cols, *rows;
    SpHelper::SpIntersect(*(A.dcsc), *(B.dcsc), cols, rows, isect1, isect2, itr1, itr2);
    
    IT kisect = static_cast<IT>(itr1-isect1);       // size of the intersection ((itr1-isect1) == (itr2-isect2))
    if(kisect == 0)
    {
        DeleteAll(isect1, isect2, cols, rows);
        return -1;
    }
    
    StackEntry< NT, std::pair<IT,IT> > * multstack;
    IT cnz = SpHelper::SpCartesian< SR > (*(A.dcsc), *(B.dcsc), kisect, isect1, isect2, multstack);  
    DeleteAll(isect1, isect2, cols, rows);

    IT mdim = A.m;  
    IT ndim = B.m;      // since B has already been transposed
    if(isZero())
    {
        dcsc = new Dcsc<IT,NT>(multstack, mdim, ndim, cnz);
    }
    else
    {
        dcsc->AddAndAssign(multstack, mdim, ndim, cnz);
    }
    nnz = dcsc->nz;

    delete [] multstack;
    return 1;   
}

template <class IT, class NT>
template <typename SR>
int SpDCCols<IT,NT>::PlusEq_AnXBn(const SpDCCols<IT,NT> & A, const SpDCCols<IT,NT> & B)
{
    if(A.isZero() || B.isZero())
    {
        return -1;  // no need to do anything
    }
    StackEntry< NT, std::pair<IT,IT> > * multstack;
    int cnz = SpHelper::SpColByCol< SR > (*(A.dcsc), *(B.dcsc), A.n, multstack);  
    
    IT mdim = A.m;  
    IT ndim = B.n;
    if(isZero())
    {
        dcsc = new Dcsc<IT,NT>(multstack, mdim, ndim, cnz);
    }
    else
    {
        dcsc->AddAndAssign(multstack, mdim, ndim, cnz);
    }
    nnz = dcsc->nz;
    
    delete [] multstack;
    return 1;   
}


template <class IT, class NT>
template <typename SR>
int SpDCCols<IT,NT>::PlusEq_AtXBn(const SpDCCols<IT,NT> & A, const SpDCCols<IT,NT> & B)
{
    std::cout << "PlusEq_AtXBn function has not been implemented yet !" << std::endl;
    return 0;
}

template <class IT, class NT>
template <typename SR>
int SpDCCols<IT,NT>::PlusEq_AtXBt(const SpDCCols<IT,NT> & A, const SpDCCols<IT,NT> & B)
{
    std::cout << "PlusEq_AtXBt function has not been implemented yet !" << std::endl;
    return 0;
}


template <class IT, class NT>
SpDCCols<IT,NT> SpDCCols<IT,NT>::operator() (IT ri, IT ci) const
{
    IT * itr = std::find(dcsc->jc, dcsc->jc + dcsc->nzc, ci);
    if(itr != dcsc->jc + dcsc->nzc)
    {
        IT irbeg = dcsc->cp[itr - dcsc->jc];
        IT irend = dcsc->cp[itr - dcsc->jc + 1];

        IT * ele = std::find(dcsc->ir + irbeg, dcsc->ir + irend, ri);
        if(ele != dcsc->ir + irend) 
        {   
            SpDCCols<IT,NT> SingEleMat(1, 1, 1, 1); // 1-by-1 matrix with 1 nonzero 
            *(SingEleMat.dcsc->numx) = dcsc->numx[ele - dcsc->ir];
            *(SingEleMat.dcsc->ir) = *ele; 
            *(SingEleMat.dcsc->jc) = *itr;
            (SingEleMat.dcsc->cp)[0] = 0;
            (SingEleMat.dcsc->cp)[1] = 1;
            return SingEleMat;
        }
        else
        {
            return SpDCCols<IT,NT>();  // 0-by-0 empty matrix
        }
    }
    else
    {
        return SpDCCols<IT,NT>();    // 0-by-0 empty matrix     
    }
}

template <class IT, class NT>
SpDCCols<IT,NT> SpDCCols<IT,NT>::operator() (const std::vector<IT> & ri, const std::vector<IT> & ci) const
{
    typedef PlusTimesSRing<NT,NT> PT;   

    IT rsize = ri.size();
    IT csize = ci.size();

    if(rsize == 0 && csize == 0)
    {
        // return an m x n matrix of complete zeros
        // since we don't know whether columns or rows are indexed
        return SpDCCols<IT,NT> (0, m, n, 0);        
    }
    else if(rsize == 0)
    {
        return ColIndex(ci);
    }
    else if(csize == 0)
    {
        SpDCCols<IT,NT> LeftMatrix(rsize, rsize, this->m, ri, true);
        return LeftMatrix.OrdColByCol< PT >(*this);
    }
    else    // this handles the (rsize=1 && csize=1) case as well
    {
        SpDCCols<IT,NT> LeftMatrix(rsize, rsize, this->m, ri, true);
        SpDCCols<IT,NT> RightMatrix(csize, this->n, csize, ci, false);
        return LeftMatrix.OrdColByCol< PT >( OrdColByCol< PT >(RightMatrix) );
    }
}

template <class IT, class NT>
std::ofstream & SpDCCols<IT,NT>::put(std::ofstream & outfile) const 
{
    if(nnz == 0)
    {
        outfile << "Matrix doesn't have any nonzeros" <<std::endl;
        return outfile;
    }
    SpTuples<IT,NT> tuples(*this); 
    outfile << tuples << std::endl;
    return outfile;
}


template <class IT, class NT>
std::ifstream & SpDCCols<IT,NT>::get(std::ifstream & infile)
{
    std::cout << "Getting... SpDCCols" << std::endl;
    IT m, n, nnz;
    infile >> m >> n >> nnz;
    SpTuples<IT,NT> tuples(nnz, m, n);        
    infile >> tuples;
    tuples.SortColBased();
        
    SpDCCols<IT,NT> object(tuples, false);
    *this = object;
    return infile;
}


template<class IT, class NT>
void SpDCCols<IT,NT>::PrintInfo(std::ofstream &  out) const
{
    out << "m: " << m ;
    out << ", n: " << n ;
    out << ", nnz: "<< nnz ;

    if(splits > 0)
    {
        out << ", local splits: " << splits << std::endl;
    }
    else
    {
        if(dcsc != NULL)
        {
            out << ", nzc: "<< dcsc->nzc << std::endl;
        }
        else
        {
            out <<", nzc: "<< 0 << std::endl;
        }
    }
}

template<class IT, class NT>
void SpDCCols<IT,NT>::PrintInfo() const
{
    std::cout << "m: " << m ;
    std::cout << ", n: " << n ;
    std::cout << ", nnz: "<< nnz ;

    if(splits > 0)
    {
        std::cout << ", local splits: " << splits << std::endl;
    }
    else
    {
        if(dcsc != NULL)
        {
            std::cout << ", nzc: "<< dcsc->nzc << std::endl;
        }
        else
        {
            std::cout <<", nzc: "<< 0 << std::endl;
        }

        if(m < PRINT_LIMIT && n < PRINT_LIMIT)  // small enough to print
        {
            NT ** A = SpHelper::allocate2D<NT>(m,n);
            for(IT i=0; i< m; ++i)
                for(IT j=0; j<n; ++j)
                    A[i][j] = NT();
            if(dcsc != NULL)
            {
                for(IT i=0; i< dcsc->nzc; ++i)
                {
                    for(IT j = dcsc->cp[i]; j<dcsc->cp[i+1]; ++j)
                    {
                        IT colid = dcsc->jc[i];
                        IT rowid = dcsc->ir[j];
                        A[rowid][colid] = dcsc->numx[j];
                    }
                }
            } 
            for(IT i=0; i< m; ++i)
            {
                for(IT j=0; j<n; ++j)
                {
                    std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(2) << A[i][j];
                    std::cout << " ";
                }
                std::cout << std::endl;
            }
            SpHelper::deallocate2D(A,m);
        }
    }
}


/****************************************************************************/
/********************* PRIVATE CONSTRUCTORS/DESTRUCTORS *********************/
/****************************************************************************/

template <class IT, class NT>
SpDCCols<IT,NT>::SpDCCols(IT nRow, IT nCol, Dcsc<IT,NT> * mydcsc)
:dcsc(mydcsc), m(nRow), n(nCol), splits(0)
{
    if (mydcsc == NULL) 
        nnz = 0;
    else
        nnz = mydcsc->nz;
}

template <class IT, class NT>
SpDCCols<IT,NT>::SpDCCols (IT size, IT nRow, IT nCol, const std::vector<IT> & indices, bool isRow)
:m(nRow), n(nCol), nnz(size), splits(0)
{
    if(size > 0)
        dcsc = new Dcsc<IT,NT>(size,indices,isRow);
    else
        dcsc = NULL; 
}


/****************************************************************************/
/************************* PRIVATE MEMBER FUNCTIONS *************************/
/****************************************************************************/

template <class IT, class NT>
inline void SpDCCols<IT,NT>::CopyDcsc(Dcsc<IT,NT> * source)
{
    // source dcsc will be NULL if number of nonzeros is zero 
    if(source != NULL)  
        dcsc = new Dcsc<IT,NT>(*source);
    else
        dcsc = NULL;
}

template <class IT, class NT>
SpDCCols<IT,NT> SpDCCols<IT,NT>::ColIndex(const std::vector<IT> & ci) const
{
    IT csize = ci.size();
    if(nnz == 0)    // nothing to index
    {
        return SpDCCols<IT,NT>(0, m, csize, 0); 
    }
    else if(ci.empty())
    {
        return SpDCCols<IT,NT>(0, m,0, 0);
    }

    // First pass for estimation
    IT estsize = 0;
    IT estnzc = 0;
    for(IT i=0, j=0;  i< dcsc->nzc && j < csize;)
    {
        if((dcsc->jc)[i] < ci[j])
        {
            ++i;
        }
        else if ((dcsc->jc)[i] > ci[j])
        {
            ++j;
        }
        else
        {
            estsize +=  (dcsc->cp)[i+1] - (dcsc->cp)[i];
            ++estnzc;
            ++i;
            ++j;
        }
    }
    
    SpDCCols<IT,NT> SubA(estsize, m, csize, estnzc);    
    if(estnzc == 0)
    {
        return SubA;        // no need to run the second pass
    }
    SubA.dcsc->cp[0] = 0;
    IT cnzc = 0;
    IT cnz = 0;
    for(IT i=0, j=0;  i < dcsc->nzc && j < csize;)
    {
        if((dcsc->jc)[i] < ci[j])
        {
            ++i;
        }
        else if ((dcsc->jc)[i] > ci[j])     // an empty column for the output
        {
            ++j;
        }
        else
        {
            IT columncount = (dcsc->cp)[i+1] - (dcsc->cp)[i];
            SubA.dcsc->jc[cnzc++] = j;
            SubA.dcsc->cp[cnzc] = SubA.dcsc->cp[cnzc-1] + columncount;
            std::copy(dcsc->ir + dcsc->cp[i], dcsc->ir + dcsc->cp[i+1], SubA.dcsc->ir + cnz);
            std::copy(dcsc->numx + dcsc->cp[i], dcsc->numx + dcsc->cp[i+1], SubA.dcsc->numx + cnz);
            cnz += columncount;
            ++i;
            ++j;
        }
    }
    return SubA;
}

template <class IT, class NT>
template <typename SR, typename NTR>
SpDCCols< IT, typename promote_trait<NT,NTR>::T_promote > SpDCCols<IT,NT>::OrdOutProdMult(const SpDCCols<IT,NTR> & rhs) const
{
    typedef typename promote_trait<NT,NTR>::T_promote T_promote;  

    if(isZero() || rhs.isZero())
    {
        return SpDCCols< IT, T_promote > (0, m, rhs.n, 0);      // return an empty matrix   
    }
    SpDCCols<IT,NTR> Btrans = rhs.TransposeConst();

    Isect<IT> *isect1, *isect2, *itr1, *itr2, *cols, *rows;
    SpHelper::SpIntersect(*dcsc, *(Btrans.dcsc), cols, rows, isect1, isect2, itr1, itr2);
    
    IT kisect = static_cast<IT>(itr1-isect1);       // size of the intersection ((itr1-isect1) == (itr2-isect2))
    if(kisect == 0)
    {
        DeleteAll(isect1, isect2, cols, rows);
        return SpDCCols< IT, T_promote > (0, m, rhs.n, 0);  
    }
    StackEntry< T_promote, std::pair<IT,IT> > * multstack;
    IT cnz = SpHelper::SpCartesian< SR > (*dcsc, *(Btrans.dcsc), kisect, isect1, isect2, multstack);  
    DeleteAll(isect1, isect2, cols, rows);

    Dcsc<IT, T_promote> * mydcsc = NULL;
    if(cnz > 0)
    {
        mydcsc = new Dcsc< IT,T_promote > (multstack, m, rhs.n, cnz);
        delete [] multstack;
    }
    return SpDCCols< IT,T_promote > (m, rhs.n, mydcsc); 
}


template <class IT, class NT>
template <typename SR, typename NTR>
SpDCCols< IT, typename promote_trait<NT,NTR>::T_promote > SpDCCols<IT,NT>::OrdColByCol(const SpDCCols<IT,NTR> & rhs) const
{
    typedef typename promote_trait<NT,NTR>::T_promote T_promote;  

    if(isZero() || rhs.isZero())
    {
        return SpDCCols<IT, T_promote> (0, m, rhs.n, 0);        // return an empty matrix   
    }
    StackEntry< T_promote, std::pair<IT,IT> > * multstack;
    IT cnz = SpHelper::SpColByCol< SR > (*dcsc, *(rhs.dcsc), n, multstack);  
    
    Dcsc<IT,T_promote > * mydcsc = NULL;
    if(cnz > 0)
    {
        mydcsc = new Dcsc< IT,T_promote > (multstack, m, rhs.n, cnz);
        delete [] multstack;
    }
    return SpDCCols< IT,T_promote > (m, rhs.n, mydcsc); 
}

}