dcsc.cpp

Go to the documentation of this file.

/****************************************************************/
/* Parallel Combinatorial BLAS Library (for Graph Computations) */
/* version 1.6 -------------------------------------------------*/
/* date: 05/15/2016 --------------------------------------------*/
/* authors: Ariful Azad, Aydin Buluc, Adam Lugowski ------------*/
/****************************************************************/
/*
 Copyright (c) 2010-2016, 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 "dcsc.h"
#include <algorithm>
#include <functional>
#include <iostream>
#include "Friends.h"
#include "SpHelper.h"

namespace combblas {

template <class IT, class NT>
Dcsc<IT,NT>::Dcsc ():cp(NULL), jc(NULL), ir(NULL), numx(NULL),nz(0), nzc(0), memowned(true){}

template <class IT, class NT>
Dcsc<IT,NT>::Dcsc (IT nnz, IT nzcol): nz(nnz),nzc(nzcol),memowned(true)
{
    assert (nz != 0);
    assert (nzc != 0);
    cp = new IT[nzc+1];
    jc = new IT[nzc];
    ir = new IT[nz];
    numx = new NT[nz];
}

template <class IT, class NT>
inline void Dcsc<IT,NT>::getindices (StackEntry<NT, std::pair<IT,IT> > * multstack, IT & rindex, IT & cindex, IT & j, IT nnz)
{
    if(j<nnz)
    {
        cindex = multstack[j].key.first;
        rindex = multstack[j].key.second;
    }
    else
    {
        rindex = std::numeric_limits<IT>::max();
        cindex = std::numeric_limits<IT>::max();
    }
    ++j;
}
    
template <class IT, class NT>
Dcsc<IT,NT> & Dcsc<IT,NT>::AddAndAssign (StackEntry<NT, std::pair<IT,IT> > * multstack, IT mdim, IT ndim, IT nnz)
{
    if(nnz == 0)    return *this;
        
    IT estnzc = nzc + nnz;
    IT estnz  = nz + nnz;
    Dcsc<IT,NT> temp(estnz, estnzc);

    IT curnzc = 0;      // number of nonzero columns constructed so far
    IT curnz = 0;
    IT i = 0;
    IT j = 0;
    IT rindex, cindex;
    getindices(multstack, rindex, cindex,j,nnz);

    temp.cp[0] = 0;
    while(i< nzc && cindex < std::numeric_limits<IT>::max())    // i runs over columns of "this",  j runs over all the nonzeros of "multstack"
    {
        if(jc[i] > cindex)
        {
            IT columncount = 0;
            temp.jc[curnzc++] = cindex;
            do
            {
                temp.ir[curnz]      = rindex;
                temp.numx[curnz++]  = multstack[j-1].value;

                getindices(multstack, rindex, cindex,j,nnz);
                ++columncount;
            }
            while(temp.jc[curnzc-1] == cindex); // loop until cindex changes

            temp.cp[curnzc] = temp.cp[curnzc-1] + columncount;
        }
        else if(jc[i] < cindex)
        {
            temp.jc[curnzc++] = jc[i++];
            for(IT k = cp[i-1]; k< cp[i]; ++k)
            {
                temp.ir[curnz]      = ir[k];
                temp.numx[curnz++]  = numx[k];
            }
            temp.cp[curnzc] = temp.cp[curnzc-1] + (cp[i] - cp[i-1]);
        }
        else    // they are equal, merge the column
        {
            temp.jc[curnzc++] = jc[i];
            IT ii = cp[i];
            IT prevnz = curnz;      
            while (ii < cp[i+1] && cindex == jc[i]) // cindex would be MAX if multstack is deplated
            {
                if (ir[ii] < rindex)
                {
                    temp.ir[curnz] = ir[ii];
                    temp.numx[curnz++] = numx[ii++];
                }
                else if (ir[ii] > rindex)
                {
                    temp.ir[curnz] = rindex;
                    temp.numx[curnz++] = multstack[j-1].value;

                    getindices(multstack, rindex, cindex,j,nnz);
                }
                else
                {
                    temp.ir[curnz] = ir[ii];
                    temp.numx[curnz++] = numx[ii++] + multstack[j-1].value;

                    getindices(multstack, rindex, cindex,j,nnz);
                }
            }
            while (ii < cp[i+1])
            {
                temp.ir[curnz] = ir[ii];
                temp.numx[curnz++] = numx[ii++];
            }
            while (cindex == jc[i])
            {
                temp.ir[curnz] = rindex;
                temp.numx[curnz++] = multstack[j-1].value;

                getindices(multstack, rindex, cindex,j,nnz);
            }
            temp.cp[curnzc] = temp.cp[curnzc-1] + curnz-prevnz;
            ++i;
        }
    }
    while(i< nzc)
    {
        temp.jc[curnzc++] = jc[i++];
        for(IT k = cp[i-1]; k< cp[i]; ++k)
        {
            temp.ir[curnz]      = ir[k];
            temp.numx[curnz++]  = numx[k];
        }
        temp.cp[curnzc] = temp.cp[curnzc-1] + (cp[i] - cp[i-1]);
    }
    while(cindex < std::numeric_limits<IT>::max())
    {
        IT columncount = 0;
        temp.jc[curnzc++] = cindex;
        do
        {
            temp.ir[curnz]      = rindex;
            temp.numx[curnz++]  = multstack[j-1].value;

            getindices(multstack, rindex, cindex,j,nnz);
            ++columncount;
        }
        while(temp.jc[curnzc-1] == cindex); // loop until cindex changes

        temp.cp[curnzc] = temp.cp[curnzc-1] + columncount;
    }
    temp.Resize(curnzc, curnz);
    *this = temp;
    return *this;
}


template <class IT, class NT>
Dcsc<IT,NT>::Dcsc (StackEntry<NT, std::pair<IT,IT> > * multstack, IT mdim, IT ndim, IT nnz): nz(nnz),memowned(true)
{
    nzc = std::min(ndim, nnz);  // nzc can't exceed any of those

    assert(nz != 0 );
    cp = new IT[nzc+1]; // to be shrinked
    jc = new IT[nzc];   // to be shrinked
    ir = new IT[nz];
    numx = new NT[nz];
    
    IT curnzc = 0;              // number of nonzero columns constructed so far
    IT cindex = multstack[0].key.first;
    IT rindex = multstack[0].key.second;

    ir[0]   = rindex;
    numx[0] = multstack[0].value;
    jc[curnzc] = cindex;
    cp[curnzc] = 0; 
    ++curnzc;

    for(IT i=1; i<nz; ++i)
    {
        cindex = multstack[i].key.first;
        rindex = multstack[i].key.second;

        ir[i]   = rindex;
        numx[i] = multstack[i].value;
        if(cindex != jc[curnzc-1])
        {
            jc[curnzc] = cindex;
            cp[curnzc++] = i;
        }
    }
    cp[curnzc] = nz;
    Resize(curnzc, nz); // only shrink cp & jc arrays
}

template <class IT, class NT>
Dcsc<IT,NT>::Dcsc (IT nnz, const std::vector<IT> & indices, bool isRow): nz(nnz),nzc(nnz),memowned(true)
{
    assert((nnz != 0) && (indices.size() == nnz));
    cp = new IT[nnz+1]; 
    jc = new IT[nnz];
    ir = new IT[nnz];
    numx = new NT[nnz];

    SpHelper::iota(cp, cp+nnz+1, 0);  // insert sequential values {0,1,2,..}
    std::fill_n(numx, nnz, static_cast<NT>(1));
    
    if(isRow)
    {
        SpHelper::iota(ir, ir+nnz, 0);  
        std::copy (indices.begin(), indices.end(), jc);
    }
    else
    {
        SpHelper::iota(jc, jc+nnz, 0);
        std::copy (indices.begin(), indices.end(), ir); 
    }
}


template <class IT, class NT>
template <typename NNT>
Dcsc<IT,NT>::operator Dcsc<IT,NNT>() const
{
    Dcsc<IT,NNT> convert(nz, nzc);  
    
    for(IT i=0; i< nz; ++i)
    {       
        convert.numx[i] = static_cast<NNT>(numx[i]);        
    }
    std::copy(ir, ir+nz, convert.ir);   // copy(first, last, result)
    std::copy(jc, jc+nzc, convert.jc);
    std::copy(cp, cp+nzc+1, convert.cp);
    return convert;
}

template <class IT, class NT>
template <typename NIT, typename NNT>
Dcsc<IT,NT>::operator Dcsc<NIT,NNT>() const
{
    Dcsc<NIT,NNT> convert(nz, nzc); 
    
    for(IT i=0; i< nz; ++i)
        convert.numx[i] = static_cast<NNT>(numx[i]);        
    for(IT i=0; i< nz; ++i)
        convert.ir[i] = static_cast<NIT>(ir[i]);
    for(IT i=0; i< nzc; ++i)
        convert.jc[i] = static_cast<NIT>(jc[i]);
    for(IT i=0; i<= nzc; ++i)
        convert.cp[i] = static_cast<NIT>(cp[i]);
    return convert;
}

template <class IT, class NT>
Dcsc<IT,NT>::Dcsc (const Dcsc<IT,NT> & rhs): nz(rhs.nz), nzc(rhs.nzc),memowned(true)
{
    if(nz > 0)
    {
        numx = new NT[nz];
        ir = new IT[nz];
        std::copy(rhs.numx, rhs.numx + nz, numx);   // numx can be a non-POD type
        std::copy(rhs.ir, rhs.ir + nz, ir);
    }
    else
    {
        numx = NULL;
        ir = NULL;
    }
    if(nzc > 0)
    {
        jc = new IT[nzc];
        cp = new IT[nzc+1];
        std::copy(rhs.jc, rhs.jc + nzc, jc);
        std::copy(rhs.cp, rhs.cp + nzc + 1, cp);
    }
    else
    {
        jc = NULL;
        cp = NULL;
    }
}

template <class IT, class NT>
Dcsc<IT,NT> & Dcsc<IT,NT>::operator=(const Dcsc<IT,NT> & rhs)
{
    if(this != &rhs)        
    {
        // make empty first !
        if(nz > 0)
        {
            delete[] numx;
            delete[] ir;    
        }
        if(nzc > 0)
        {
            delete[] jc;
            delete[] cp;
        }
        nz = rhs.nz;
        nzc = rhs.nzc;
        if(nz > 0)
        {
            numx = new NT[nz];
            ir = new IT[nz];
            std::copy(rhs.numx, rhs.numx + nz, numx);   // numx can be a non-POD type
            std::copy(rhs.ir, rhs.ir + nz, ir);
        }
        else
        {
            numx = NULL;
            ir = NULL;
        }
        if(nzc > 0)
        {
            jc = new IT[nzc];
                    cp = new IT[nzc+1];
                    std::copy(rhs.jc, rhs.jc + nzc, jc);
                    std::copy(rhs.cp, rhs.cp + nzc + 1, cp);
        }
        else
        {
            jc = NULL;
            cp = NULL;
        }
    }
    return *this;
}

template <class IT, class NT>
Dcsc<IT, NT> & Dcsc<IT,NT>::operator+=(const Dcsc<IT,NT> & rhs) // add and assign operator
{
    IT estnzc = nzc + rhs.nzc;
    IT estnz  = nz + rhs.nz;
    Dcsc<IT,NT> temp(estnz, estnzc);

    IT curnzc = 0;
    IT curnz = 0;
    IT i = 0;
    IT j = 0;
    temp.cp[0] = 0;
    while(i< nzc && j<rhs.nzc)
    {
        if(jc[i] > rhs.jc[j])
        {
            temp.jc[curnzc++] = rhs.jc[j++];
            for(IT k = rhs.cp[j-1]; k< rhs.cp[j]; ++k)
            {
                temp.ir[curnz]      = rhs.ir[k];
                temp.numx[curnz++]  = rhs.numx[k];
            }
            temp.cp[curnzc] = temp.cp[curnzc-1] + (rhs.cp[j] - rhs.cp[j-1]);
        }
        else if(jc[i] < rhs.jc[j])
        {
            temp.jc[curnzc++] = jc[i++];
            for(IT k = cp[i-1]; k< cp[i]; k++)
            {
                temp.ir[curnz]      = ir[k];
                temp.numx[curnz++]  = numx[k];
            }
            temp.cp[curnzc] = temp.cp[curnzc-1] + (cp[i] - cp[i-1]);
        }
        else
        {
            temp.jc[curnzc++] = jc[i];
            IT ii = cp[i];
            IT jj = rhs.cp[j];
            IT prevnz = curnz;      
            while (ii < cp[i+1] && jj < rhs.cp[j+1])
            {
                if (ir[ii] < rhs.ir[jj])
                {
                    temp.ir[curnz] = ir[ii];
                    temp.numx[curnz++] = numx[ii++];
                }
                else if (ir[ii] > rhs.ir[jj])
                {
                    temp.ir[curnz] = rhs.ir[jj];
                    temp.numx[curnz++] = rhs.numx[jj++];
                }
                else
                {
                    temp.ir[curnz] = ir[ii];
                    temp.numx[curnz++] = numx[ii++] + rhs.numx[jj++];       // might include zeros
                }
            }
            while (ii < cp[i+1])
            {
                temp.ir[curnz] = ir[ii];
                temp.numx[curnz++] = numx[ii++];
            }
            while (jj < rhs.cp[j+1])
            {
                temp.ir[curnz] = rhs.ir[jj];
                temp.numx[curnz++] = rhs.numx[jj++];
            }
            temp.cp[curnzc] = temp.cp[curnzc-1] + curnz-prevnz;
            ++i;
            ++j;
        }
    }
    while(i< nzc)
    {
        temp.jc[curnzc++] = jc[i++];
        for(IT k = cp[i-1]; k< cp[i]; ++k)
        {
            temp.ir[curnz]  = ir[k];
            temp.numx[curnz++] = numx[k];
        }
        temp.cp[curnzc] = temp.cp[curnzc-1] + (cp[i] - cp[i-1]);
    }
    while(j < rhs.nzc)
    {
        temp.jc[curnzc++] = rhs.jc[j++];
        for(IT k = rhs.cp[j-1]; k< rhs.cp[j]; ++k)
        {
            temp.ir[curnz]  = rhs.ir[k];
            temp.numx[curnz++]  = rhs.numx[k];
        }
        temp.cp[curnzc] = temp.cp[curnzc-1] + (rhs.cp[j] - rhs.cp[j-1]);
    }
    temp.Resize(curnzc, curnz);
    *this = temp;
    return *this;
}

template <class IT, class NT>
bool Dcsc<IT,NT>::operator==(const Dcsc<IT,NT> & rhs)
{
    if(nzc != rhs.nzc) return false;
    bool same = std::equal(cp, cp+nzc+1, rhs.cp);
    same = same && std::equal(jc, jc+nzc, rhs.jc);
    same = same && std::equal(ir, ir+nz, rhs.ir);
    
#ifdef DEBUG
  std::vector<NT> error(nz);
  std::transform(numx, numx+nz, rhs.numx, error.begin(), absdiff<NT>());
  std::vector< std::pair<NT, NT> > error_original_pair(nz);
    for(IT i=0; i < nz; ++i)
        error_original_pair[i] = std::make_pair(error[i], numx[i]);
    if(error_original_pair.size() > 10) // otherwise would crush for small data
    {
        partial_sort(error_original_pair.begin(), error_original_pair.begin()+10, error_original_pair.end(), std::greater< std::pair<NT,NT> >());
    std::cout << "Highest 10 different entries are: " << std::endl;
        for(IT i=0; i < 10; ++i)
            std::cout << "Diff: " << error_original_pair[i].first << " on " << error_original_pair[i].second << std::endl;
    }
    else
    {
        sort(error_original_pair.begin(), error_original_pair.end(), std::greater< std::pair<NT,NT> >());
        std::cout << "Highest different entries are: " << std::endl;
        for(typename std::vector< std::pair<NT, NT> >::iterator it=error_original_pair.begin(); it != error_original_pair.end(); ++it)
            std::cout << "Diff: " << it->first << " on " << it->second << std::endl;
    }
    std::cout << "Same before num: " << same << std::endl;
#endif

    ErrorTolerantEqual<NT> epsilonequal;
    same = same && std::equal(numx, numx+nz, rhs.numx, epsilonequal );
    return same;
}

template <class IT, class NT>
void Dcsc<IT,NT>::EWiseMult(const Dcsc<IT,NT> & rhs, bool exclude)  
{
    // We have a class with a friend function and a member function with the same name. Calling the friend function from the member function 
    // might (if the signature is the same) give compilation errors if not preceded by :: that denotes the global scope.
    *this = combblas::EWiseMult((*this), &rhs, exclude);    // call the binary version
}


template <class IT, class NT>
template <typename _UnaryOperation, typename GlobalIT>
Dcsc<IT,NT>* Dcsc<IT,NT>::PruneI(_UnaryOperation __unary_op, bool inPlace, GlobalIT rowOffset, GlobalIT colOffset)
{
    // Two-pass algorithm
    IT prunednnz = 0;
    IT prunednzc = 0;
    for(IT i=0; i<nzc; ++i)
    {
        bool colexists = false;
        for(IT j=cp[i]; j < cp[i+1]; ++j)
        {
            if(!(__unary_op(make_tuple(rowOffset+ir[j], colOffset+jc[i], numx[j]))))    // keep this nonzero
            {
                ++prunednnz;
                colexists = true;
            }
        }
        if(colexists)   ++prunednzc;
    }
    IT * oldcp = cp; 
    IT * oldjc = jc;
    IT * oldir = ir;    
    NT * oldnumx = numx;    

    cp = new IT[prunednzc+1];
    jc = new IT[prunednzc];
    ir = new IT[prunednnz];
    numx = new NT[prunednnz];

    IT cnzc = 0;
    IT cnnz = 0;
    cp[cnzc] = 0;
    for(IT i=0; i<nzc; ++i)
    {
        for(IT j = oldcp[i]; j < oldcp[i+1]; ++j)
        {
            if(!(__unary_op(make_tuple(rowOffset+oldir[j], colOffset+oldjc[i], oldnumx[j])))) // keep this nonzero
            {
                ir[cnnz] = oldir[j];    
                numx[cnnz++] =  oldnumx[j];
            }
        }
        if(cnnz > cp[cnzc])
        {
            jc[cnzc] = oldjc[i];
            cp[cnzc+1] = cnnz;
            ++cnzc;
        }
    }
    assert(cnzc == prunednzc);
    assert(cnnz == prunednnz);
    if (inPlace)
    {
        // delete the memory pointed by previous pointers
        DeleteAll(oldnumx, oldir, oldjc, oldcp);
        nz = cnnz;
        nzc = cnzc;
        return NULL;
    }
    else
    {
        // create a new object to store the data
        Dcsc<IT,NT>* ret = new Dcsc<IT,NT>();
        ret->cp = cp; 
        ret->jc = jc;
        ret->ir = ir;   
        ret->numx = numx;
        ret->nz = cnnz;
        ret->nzc = cnzc;

        // put the previous pointers back       
        cp = oldcp;
        jc = oldjc;
        ir = oldir;
        numx = oldnumx;
        
        return ret;
    }
}

template <class IT, class NT>
template <typename _UnaryOperation>
Dcsc<IT,NT>* Dcsc<IT,NT>::Prune(_UnaryOperation __unary_op, bool inPlace)
{
    // Two-pass algorithm
    IT prunednnz = 0;
    IT prunednzc = 0;
    for(IT i=0; i<nzc; ++i)
    {
        bool colexists = false;
        for(IT j=cp[i]; j < cp[i+1]; ++j)
        {
            if(!(__unary_op(numx[j])))  // keep this nonzero
            {
                ++prunednnz;
                colexists = true;
            }
        }
        if(colexists)   ++prunednzc;
    }
    IT * oldcp = cp; 
    IT * oldjc = jc;
    IT * oldir = ir;    
    NT * oldnumx = numx;    

    cp = new IT[prunednzc+1];
    jc = new IT[prunednzc];
    ir = new IT[prunednnz];
    numx = new NT[prunednnz];

    IT cnzc = 0;
    IT cnnz = 0;
    cp[cnzc] = 0;
    for(IT i=0; i<nzc; ++i)
    {
        for(IT j = oldcp[i]; j < oldcp[i+1]; ++j)
        {
            if(!(__unary_op(oldnumx[j]))) // keep this nonzero
            {
                ir[cnnz] = oldir[j];    
                numx[cnnz++] =  oldnumx[j];
            }
        }
        if(cnnz > cp[cnzc])
        {
            jc[cnzc] = oldjc[i];
            cp[cnzc+1] = cnnz;
            ++cnzc;
        }
    }
    assert(cnzc == prunednzc);
    assert(cnnz == prunednnz);
    if (inPlace)
    {
        // delete the memory pointed by previous pointers
        DeleteAll(oldnumx, oldir, oldjc, oldcp);
        nz = cnnz;
        nzc = cnzc;
        return NULL;
    }
    else
    {
        // create a new object to store the data
        Dcsc<IT,NT>* ret = new Dcsc<IT,NT>();
        ret->cp = cp; 
        ret->jc = jc;
        ret->ir = ir;   
        ret->numx = numx;
        ret->nz = cnnz;
        ret->nzc = cnzc;

        // put the previous pointers back       
        cp = oldcp;
        jc = oldjc;
        ir = oldir;
        numx = oldnumx;
        
        return ret;
    }
}


template <class IT, class NT>
template <typename _BinaryOperation>
Dcsc<IT,NT>* Dcsc<IT,NT>::PruneColumn(NT* pvals, _BinaryOperation __binary_op, bool inPlace)
{
    // Two-pass algorithm
    IT prunednnz = 0;
    IT prunednzc = 0;
    for(IT i=0; i<nzc; ++i)
    {
        bool colexists = false;
        for(IT j=cp[i]; j < cp[i+1]; ++j)
        {
            IT colid = jc[i];
            if(!(__binary_op(numx[j], pvals[colid])))   // keep this nonzero
            {
                ++prunednnz;
                colexists = true;
            }
        }
        if(colexists)   ++prunednzc;
    }
    IT * oldcp = cp;
    IT * oldjc = jc;
    IT * oldir = ir;
    NT * oldnumx = numx;
    
    cp = new IT[prunednzc+1];
    jc = new IT[prunednzc];
    ir = new IT[prunednnz];
    numx = new NT[prunednnz];
    
    IT cnzc = 0;
    IT cnnz = 0;
    cp[cnzc] = 0;
    for(IT i=0; i<nzc; ++i)
    {
        for(IT j = oldcp[i]; j < oldcp[i+1]; ++j)
        {
            IT colid = oldjc[i];
            if(!(__binary_op(oldnumx[j], pvals[colid])))
            {
                ir[cnnz] = oldir[j];
                numx[cnnz++] =  oldnumx[j];
            }
        }
        if(cnnz > cp[cnzc])
        {
            jc[cnzc] = oldjc[i];
            cp[cnzc+1] = cnnz;
            ++cnzc;
        }
    }
    assert(cnzc == prunednzc);
    if (inPlace)
    {
        // delete the memory pointed by previous pointers
        DeleteAll(oldnumx, oldir, oldjc, oldcp);
        nz = cnnz;
        nzc = cnzc;
        return NULL;
    }
    else
    {
        // create a new object to store the data
        Dcsc<IT,NT>* ret = new Dcsc<IT,NT>();
        ret->cp = cp; 
        ret->jc = jc;
        ret->ir = ir;   
        ret->numx = numx;
        ret->nz = cnnz;
        ret->nzc = cnzc;
        
        // put the previous pointers back       
        cp = oldcp;
        jc = oldjc;
        ir = oldir;
        numx = oldnumx;
        
        return ret;
    }
}


// prune selected columns indexed by pinds
template <class IT, class NT>
template <typename _BinaryOperation>
Dcsc<IT,NT>* Dcsc<IT,NT>::PruneColumn(IT* pinds, NT* pvals, _BinaryOperation __binary_op, bool inPlace)
{
    // Two-pass algorithm
    IT prunednnz = 0;
    IT prunednzc = 0;
    IT k = 0;
    for(IT i=0; i<nzc; ++i)
    {
        bool colexists = false;
        IT colid = jc[i];
        if(colid==pinds[k]) // pinds is sorted
        {
            for(IT j=cp[i]; j < cp[i+1]; ++j)
            {
                if(!(__binary_op(numx[j], pvals[k])))   // keep this nonzero
                {
                    ++prunednnz;
                    colexists = true;
                }
            }
            k++;
        }
        else // untouched columns
        {
            colexists = true;
            prunednnz += (cp[i+1] - cp[i]);
        }
        if(colexists)   ++prunednzc;
    }
    IT * oldcp = cp;
    IT * oldjc = jc;
    IT * oldir = ir;
    NT * oldnumx = numx;
    
    cp = new IT[prunednzc+1];
    jc = new IT[prunednzc];
    ir = new IT[prunednnz];
    numx = new NT[prunednnz];
    
    IT cnzc = 0;
    IT cnnz = 0;
    cp[cnzc] = 0;
    k = 0;
    for(IT i=0; i<nzc; ++i)
    {
        IT colid = oldjc[i];
        if(colid==pinds[k]) // prunned columns
        {
            for(IT j = oldcp[i]; j < oldcp[i+1]; ++j)
            {
                if(!(__binary_op(oldnumx[j], pvals[k])))
                {
                    ir[cnnz] = oldir[j];
                    numx[cnnz++] =  oldnumx[j];
                }
            }
            k++;
        }
        else // copy other columns
        {
            for(IT j = oldcp[i]; j < oldcp[i+1]; ++j)
            {
                ir[cnnz] = oldir[j];
                numx[cnnz++] =  oldnumx[j];
            }
        }
        if(cnnz > cp[cnzc])
        {
            jc[cnzc] = oldjc[i];
            cp[cnzc+1] = cnnz;
            ++cnzc;
        }
    }
    assert(cnzc == prunednzc);
    if (inPlace)
    {
        // delete the memory pointed by previous pointers
        DeleteAll(oldnumx, oldir, oldjc, oldcp);
        nz = cnnz;
        nzc = cnzc;
        return NULL;
    }
    else
    {
        // create a new object to store the data
        Dcsc<IT,NT>* ret = new Dcsc<IT,NT>();
        ret->cp = cp;
        ret->jc = jc;
        ret->ir = ir;
        ret->numx = numx;
        ret->nz = cnnz;
        ret->nzc = cnzc;
        
        // put the previous pointers back
        cp = oldcp;
        jc = oldjc;
        ir = oldir;
        numx = oldnumx;
        
        return ret;
    }
}

template <class IT, class NT>
void Dcsc<IT,NT>::EWiseScale(NT ** scaler)  
{
    for(IT i=0; i<nzc; ++i)
    {
        IT colid = jc[i];
        for(IT j=cp[i]; j < cp[i+1]; ++j)
        {
            IT rowid = ir[j];
            numx[j] *= scaler[rowid][colid];
        }
    }
}

template <class IT, class NT>
template <typename _BinaryOperation>
void Dcsc<IT,NT>::UpdateDense(NT ** array, _BinaryOperation __binary_op) const
{
    for(IT i=0; i<nzc; ++i)
    {
        IT colid = jc[i];
        for(IT j=cp[i]; j < cp[i+1]; ++j)
        {
            IT rowid = ir[j];
            array[rowid][colid] = __binary_op(array[rowid][colid], numx[j]);
        }
    }
}

template <class IT, class NT>
IT Dcsc<IT,NT>::ConstructAux(IT ndim, IT * & aux) const
{
    float cf  = static_cast<float>(ndim+1) / static_cast<float>(nzc);
    IT colchunks = static_cast<IT> ( ceil( static_cast<float>(ndim+1) / ceil(cf)) );

    aux  = new IT[colchunks+1]; 

    IT chunksize    = static_cast<IT>(ceil(cf));
    IT reg      = 0;
    IT curchunk = 0;
    aux[curchunk++] = 0;
    for(IT i = 0; i< nzc; ++i)
    {
        if(jc[i] >= curchunk * chunksize)       // beginning of the next chunk
        {
            while(jc[i] >= curchunk * chunksize)    // consider any empty chunks
            {
                aux[curchunk++] = reg;
            }
        }
        reg = i+1;
    }
    while(curchunk <= colchunks)
    {
        aux[curchunk++] = reg;
    }
    return colchunks;
}

template <class IT, class NT>
void Dcsc<IT,NT>::Resize(IT nzcnew, IT nznew)
{
    if(nzcnew == 0)
    {
        delete[] jc;
        delete[] cp;
        nzc = 0;
    }
    if(nznew == 0)
    {
        delete[] ir;
        delete[] numx;
        nz = 0;
    }
    if ( nzcnew == 0 && nznew == 0)
    {
        return; 
    }
    if (nzcnew != nzc)  
    {
        IT * tmpcp = cp; 
        IT * tmpjc = jc;
        cp = new IT[nzcnew+1];
        jc = new IT[nzcnew];    
        if(nzcnew > nzc)    // Grow it (copy all of the old elements)
        {
            std::copy(tmpcp, tmpcp+nzc+1, cp);  // copy(first, end, result)
            std::copy(tmpjc, tmpjc+nzc, jc);
        }
        else        // Shrink it (copy only a portion of the old elements)
        {
            std::copy(tmpcp, tmpcp+nzcnew+1, cp);   
            std::copy(tmpjc, tmpjc+nzcnew, jc);
        }
        delete[] tmpcp; // delete the memory pointed by previous pointers
        delete[] tmpjc;
        nzc = nzcnew;
    }
    if (nznew != nz)
    {   
        NT * tmpnumx = numx; 
        IT * tmpir = ir;
        numx = new NT[nznew];   
        ir = new IT[nznew];
        if(nznew > nz)  // Grow it (copy all of the old elements)
        {
            std::copy(tmpnumx, tmpnumx+nz, numx);   // numx can be non-POD
            std::copy(tmpir, tmpir+nz, ir);
        }
        else    // Shrink it (copy only a portion of the old elements)
        {
            std::copy(tmpnumx, tmpnumx+nznew, numx);    
            std::copy(tmpir, tmpir+nznew, ir);
        }
        delete[] tmpnumx;   // delete the memory pointed by previous pointers
        delete[] tmpir;
        nz = nznew;
    }
}

template<class IT, class NT>
IT Dcsc<IT,NT>::AuxIndex(const IT colind, bool & found, IT * aux, IT csize) const
{
    IT base = static_cast<IT>(floor((float) (colind/csize)));
    IT start = aux[base];
    IT end = aux[base+1];

    IT * itr = std::find(jc + start, jc + end, colind);
    
    found = (itr != jc + end);
    return (itr-jc);
}

template<class IT, class NT>
void Dcsc<IT,NT>::Split(Dcsc<IT,NT> * & A, Dcsc<IT,NT> * & B, IT cut)
{
    IT * itr = std::lower_bound(jc, jc+nzc, cut);
    IT pos = itr - jc;

    if(cp[pos] == 0)
    {
        A = NULL;
    }
    else
    {
        A = new Dcsc<IT,NT>(cp[pos], pos);
        std::copy(jc, jc+pos, A->jc);
        std::copy(cp, cp+pos+1, A->cp);
        std::copy(ir, ir+cp[pos], A->ir);
        std::copy(numx, numx + cp[pos], A->numx);   // copy(first, last, result)
    }   
    if(nz-cp[pos] == 0)
    {
        B = NULL;
    }
    else
    {
        B = new Dcsc<IT,NT>(nz-cp[pos], nzc-pos);
        std::copy(jc+pos, jc+ nzc, B->jc);
        transform(B->jc, B->jc + (nzc-pos), B->jc, bind2nd(std::minus<IT>(), cut));
        std::copy(cp+pos, cp+nzc+1, B->cp);
        transform(B->cp, B->cp + (nzc-pos+1), B->cp, bind2nd(std::minus<IT>(), cp[pos]));
        std::copy(ir+cp[pos], ir+nz, B->ir);
        std::copy(numx+cp[pos], numx+nz, B->numx);  // copy(first, last, result)
    }
}

template<class IT, class NT>
void Dcsc<IT,NT>::ColSplit(std::vector< Dcsc<IT,NT>* > & parts, std::vector<IT> & cuts)
{
    IT * jcbegin = jc;
    std::vector<IT> pos; // pos has "parts-1" entries
    for(auto cutpoint = cuts.begin(); cutpoint != cuts.end(); ++cutpoint)
    {
        IT * itr = std::lower_bound(jcbegin, jc+nzc, *cutpoint);
        pos.push_back(itr - jc);
        jcbegin = itr;  // so that lower_bound searches a smaller vector
    }
    
    if(cp[pos[0]] == 0) // first piece
    {
        parts[0] = NULL;
    }
    else
    {
        parts[0] = new Dcsc<IT,NT>(cp[pos[0]], pos[0]); // Dcsc(nnz, nzc)
        std::copy(jc, jc+pos[0], parts[0]->jc);    // std::copy
        std::copy(cp, cp+pos[0]+1, parts[0]->cp);
        std::copy(ir, ir+cp[pos[0]], parts[0]->ir);
        std::copy(numx, numx + cp[pos[0]], parts[0]->numx); // copy(first, last, result)
    }
    int ncuts =  cuts.size(); // all except last piece
    for(int i=1; i< ncuts; ++i) // treat the first piece differently
    {
        if(cp[pos[i]] - cp[pos[i-1]] == 0)
        {
            parts[i] =  NULL;
        }
        else
        {
            parts[i] = new Dcsc<IT,NT>(cp[pos[i]] - cp[pos[i-1]], pos[i] - pos[i-1]); // Dcsc(nnz, nzc)
            std::copy(jc+pos[i-1], jc+pos[i], parts[i]->jc);    // std::copy
            transform(parts[i]->jc, parts[i]->jc + (pos[i]-pos[i-1]), parts[i]->jc, bind2nd(std::minus<IT>(), cuts[i-1]));  // cuts[i-1] is well defined as i>=1

            std::copy(cp+pos[i-1], cp+pos[i]+1, parts[i]->cp);
            transform(parts[i]->cp, parts[i]->cp + (pos[i]-pos[i-1]+1), parts[i]->cp, bind2nd(std::minus<IT>(), cp[pos[i-1]]));

            std::copy(ir+cp[pos[i-1]], ir+cp[pos[i]], parts[i]->ir);
            std::copy(numx+cp[pos[i-1]], numx + cp[pos[i]], parts[i]->numx);    // copy(first, last, result)
        }
    }
    if(nz - cp[pos[ncuts-1]] == 0)
    {
        parts[ncuts] = NULL;
    }
    else
    {
        parts[ncuts] = new Dcsc<IT,NT>(nz-cp[pos[ncuts-1]], nzc-pos[ncuts-1]);  // ncuts = npieces -1
        std::copy(jc+pos[ncuts-1], jc+ nzc, parts[ncuts]->jc);
        transform(parts[ncuts]->jc, parts[ncuts]->jc + (nzc-pos[ncuts-1]), parts[ncuts]->jc, bind2nd(std::minus<IT>(), cuts[ncuts-1]));
        
        std::copy(cp+pos[ncuts-1], cp+nzc+1, parts[ncuts]->cp);
        transform(parts[ncuts]->cp, parts[ncuts]->cp + (nzc-pos[ncuts-1]+1), parts[ncuts]->cp, bind2nd(std::minus<IT>(), cp[pos[ncuts-1]]));
        std::copy(ir+cp[pos[ncuts-1]], ir+nz, parts[ncuts]->ir);
        std::copy(numx+cp[pos[ncuts-1]], numx+nz, parts[ncuts]->numx);
    }
}


// Assumes A and B are not NULL
// When any is NULL, this function is not called anyway
template<class IT, class NT>
void Dcsc<IT,NT>::Merge(const Dcsc<IT,NT> * A, const Dcsc<IT,NT> * B, IT cut)
{
    assert((A != NULL) && (B != NULL));     // handled at higher level
    IT cnz = A->nz + B->nz;
    IT cnzc = A->nzc + B->nzc;
    if(cnz > 0)
    {
        *this = Dcsc<IT,NT>(cnz, cnzc);     // safe, because "this" can not be NULL inside a member function

        std::copy(A->jc, A->jc + A->nzc, jc);   // copy(first, last, result)
        std::copy(B->jc, B->jc + B->nzc, jc + A->nzc);
        transform(jc + A->nzc, jc + cnzc, jc + A->nzc, bind2nd(std::plus<IT>(), cut));

        std::copy(A->cp, A->cp + A->nzc, cp);
        std::copy(B->cp, B->cp + B->nzc +1, cp + A->nzc);
        transform(cp + A->nzc, cp+cnzc+1, cp + A->nzc, bind2nd(std::plus<IT>(), A->cp[A->nzc]));
    
        std::copy(A->ir, A->ir + A->nz, ir);
        std::copy(B->ir, B->ir + B->nz, ir + A->nz);

        // since numx is potentially non-POD, we use std::copy
        std::copy(A->numx, A->numx + A->nz, numx);
        std::copy(B->numx, B->numx + B->nz, numx + A->nz);
    }
}


template<class IT, class NT>
void Dcsc<IT,NT>::ColConcatenate(std::vector< Dcsc<IT,NT>* > & parts, std::vector<IT> & offsets)
{
    IT cnz = 0;
    IT cnzc = 0;
    size_t nmembers = parts.size();
    for(size_t i=0; i< nmembers; ++i)
    {
        cnz += parts[i]->nz;
        cnzc += parts[i]->nzc;
    }
    if(cnz > 0)
    {
        *this = Dcsc<IT,NT>(cnz, cnzc);     // safe, because "this" can not be NULL inside a member function
        
        IT run_nz = 0;
        IT run_nzc = 0;
        for(size_t i=0; i< nmembers; ++i)
        {
            std::copy(parts[i]->jc, parts[i]->jc + parts[i]->nzc, jc + run_nzc);
            transform(jc + run_nzc, jc + run_nzc + parts[i]->nzc, jc + run_nzc, bind2nd(std::plus<IT>(), offsets[i]));
            
            // remember: cp[nzc] = nnz
            std::copy(parts[i]->cp, parts[i]->cp + parts[i]->nzc, cp + run_nzc);
            transform(cp + run_nzc, cp + run_nzc + parts[i]->nzc, cp + run_nzc, bind2nd(std::plus<IT>(),run_nz));
            
            std::copy(parts[i]->ir, parts[i]->ir + parts[i]->nz, ir + run_nz);
            std::copy(parts[i]->numx, parts[i]->numx + parts[i]->nz, numx + run_nz);
            
            run_nzc += parts[i]->nzc;
            run_nz += parts[i]->nz;
        }
        // adjust the last pointer
        cp[run_nzc] = run_nz;
    }
}

template<class IT, class NT>
template<class VT>  
void Dcsc<IT,NT>::FillColInds(const VT * colnums, IT nind, std::vector< std::pair<IT,IT> > & colinds, IT * aux, IT csize) const
{
    if ( aux == NULL || (nzc / nind) < THRESHOLD)       // use scanning indexing
    {
        IT mink = std::min(nzc, nind);
        std::pair<IT,IT> * isect = new std::pair<IT,IT>[mink];
        std::pair<IT,IT> * range1 = new std::pair<IT,IT>[nzc];
        std::pair<IT,IT> * range2 = new std::pair<IT,IT>[nind];
        
        for(IT i=0; i < nzc; ++i)
        {
            range1[i] = std::make_pair(jc[i], i);   // get the actual nonzero value and the index to the ith nonzero
        }
        for(IT i=0; i < nind; ++i)
        {
            range2[i] = std::make_pair(static_cast<IT>(colnums[i]), 0); // second is dummy as all the intersecting elements are copied from the first range
        }

        std::pair<IT,IT> * itr = set_intersection(range1, range1 + nzc, range2, range2+nind, isect, SpHelper::first_compare<IT> );
        // isect now can iterate on a subset of the elements of range1
        // meaning that the intersection can be accessed directly by isect[i] instead of range1[isect[i]]
        // this is because the intersecting elements are COPIED to the output range "isect"

        IT kisect = static_cast<IT>(itr-isect);     // size of the intersection 
        for(IT j=0, i =0; j< nind; ++j)
        {
            // the elements represented by jc[isect[i]] are a subset of the elements represented by colnums[j]
            if( i == kisect || isect[i].first != static_cast<IT>(colnums[j]))
            {
                // not found, signal by setting first = second
                colinds[j].first = 0;
                colinds[j].second = 0;  
            }
            else    // i < kisect && dcsc->jc[isect[i]] == colnums[j]
            {
                IT p = isect[i++].second;
                colinds[j].first = cp[p];
                colinds[j].second = cp[p+1];
            }
        }
        DeleteAll(isect, range1, range2);
    }
    else        // use aux based indexing
    {
        bool found;
        for(IT j =0; j< nind; ++j)
        {
            IT pos = AuxIndex(static_cast<IT>(colnums[j]), found, aux, csize);
            if(found)
            {
                colinds[j].first = cp[pos];
                colinds[j].second = cp[pos+1];
            }
            else    // not found, signal by setting first = second
            {
                colinds[j].first = 0;
                colinds[j].second = 0;
            }
        }
    }
}


template <class IT, class NT>
Dcsc<IT,NT>::~Dcsc()
{
    if(nz > 0)          // dcsc may be empty
    {
        delete[] numx;
        delete[] ir;
    }
    if(nzc > 0)
    {
        delete[] jc;
        delete[] cp;
    }
}

}