SpImpl.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 "SpImpl.h"
#include "SpParHelper.h"
#include "PBBS/radixSort.h"
#include "Tommy/tommyhashdyn.h"
 
namespace combblas {

template <class SR, class IT, class NUM, class IVT, class OVT>
void SpImpl<SR,IT,NUM,IVT,OVT>::SpMXSpV(const Dcsc<IT,NUM> & Adcsc, int32_t mA, const int32_t * indx, const IVT * numx, int32_t veclen,  
            std::vector<int32_t> & indy, std::vector< OVT > & numy)
{
    int32_t hsize = 0;      
    // colinds dereferences A.ir (valid from colinds[].first to colinds[].second)
    std::vector< std::pair<IT,IT> > colinds( (IT) veclen);      
    Adcsc.FillColInds(indx, (IT) veclen, colinds, NULL, 0); // csize is irrelevant if aux is NULL   
 
    if(sizeof(NUM) > sizeof(OVT))   // ABAB: include a filtering based runtime choice as well?
    {
        HeapEntry<IT, OVT> * wset = new HeapEntry<IT, OVT>[veclen]; 
        for(IT j =0; j< veclen; ++j)        // create the initial heap 
        {
            while(colinds[j].first != colinds[j].second )   // iterate until finding the first entry within this column that passes the filter
            {
                OVT mrhs = SR::multiply(Adcsc.numx[colinds[j].first], numx[j]);
                if(SR::returnedSAID())
                {
                    ++(colinds[j].first);   // increment the active row index within the jth column
                }
                else
                {
                    wset[hsize++] = HeapEntry< IT,OVT > ( Adcsc.ir[colinds[j].first], j, mrhs);
                    break;  // this column successfully inserted an entry to the heap
                }
            } 
        }
        std::make_heap(wset, wset+hsize);
        while(hsize > 0)
        {
            std::pop_heap(wset, wset + hsize);          // result is stored in wset[hsize-1]
            IT locv = wset[hsize-1].runr;       // relative location of the nonzero in sparse column vector 
            if((!indy.empty()) && indy.back() == wset[hsize-1].key) 
            {
                numy.back() = SR::add(numy.back(), wset[hsize-1].num);
            }
            else
            {
                indy.push_back( (int32_t) wset[hsize-1].key);
                numy.push_back(wset[hsize-1].num);  
            }
            bool pushed = false;
            // invariant: if ++(colinds[locv].first) == colinds[locv].second, then locv will not appear again in the heap
            while ( (++(colinds[locv].first)) != colinds[locv].second ) // iterate until finding another passing entry
            {
                OVT mrhs =  SR::multiply(Adcsc.numx[colinds[locv].first], numx[locv]);
                if(!SR::returnedSAID())
                                {
                    wset[hsize-1].key = Adcsc.ir[colinds[locv].first];
                    wset[hsize-1].num = mrhs;
                    std::push_heap(wset, wset+hsize);   // runr stays the same
                    pushed = true;
                    break;
                }
            }
            if(!pushed) --hsize;
        }
        delete [] wset;
    }
    
    else
    {
        HeapEntry<IT, NUM> * wset = new HeapEntry<IT, NUM>[veclen]; 
        for(IT j =0; j< veclen; ++j)        // create the initial heap 
        {
            if(colinds[j].first != colinds[j].second)   // current != end
            {
                wset[hsize++] = HeapEntry< IT,NUM > ( Adcsc.ir[colinds[j].first], j, Adcsc.numx[colinds[j].first]);  // HeapEntry(key, run, num)
            } 
        }   
        std::make_heap(wset, wset+hsize);
        while(hsize > 0)
        {
            std::pop_heap(wset, wset + hsize);          // result is stored in wset[hsize-1]
            IT locv = wset[hsize-1].runr;       // relative location of the nonzero in sparse column vector 
            OVT mrhs = SR::multiply(wset[hsize-1].num, numx[locv]); 
        
            if (!SR::returnedSAID())
            {
                if((!indy.empty()) && indy.back() == wset[hsize-1].key) 
                {
                    numy.back() = SR::add(numy.back(), mrhs);
                }
                else
                {
                    indy.push_back( (int32_t) wset[hsize-1].key);
                    numy.push_back(mrhs);   
                }
            }
 
            if( (++(colinds[locv].first)) != colinds[locv].second)  // current != end
            {
                // runr stays the same !
                wset[hsize-1].key = Adcsc.ir[colinds[locv].first];
                wset[hsize-1].num = Adcsc.numx[colinds[locv].first];  
                std::push_heap(wset, wset+hsize);
            }
            else        --hsize;
        }
        delete [] wset;
    }
}
 
 

template <class SR, class IT, class IVT, class OVT>
void SpImpl<SR,IT,bool,IVT,OVT>::SpMXSpV(const Dcsc<IT,bool> & Adcsc, int32_t mA, const int32_t * indx, const IVT * numx, int32_t veclen,  
            std::vector<int32_t> & indy, std::vector<OVT> & numy)
{   
    IT inf = std::numeric_limits<IT>::min();
    IT sup = std::numeric_limits<IT>::max(); 
    KNHeap< IT, IVT > sHeap(sup, inf);  // max size: flops
 
    IT k = 0;   // index to indx vector
    IT i = 0;   // index to columns of matrix
    while(i< Adcsc.nzc && k < veclen)
    {
        if(Adcsc.jc[i] < indx[k]) ++i;
        else if(indx[k] < Adcsc.jc[i]) ++k;
        else
        {
            for(IT j=Adcsc.cp[i]; j < Adcsc.cp[i+1]; ++j)   // for all nonzeros in this column
            {
                sHeap.insert(Adcsc.ir[j], numx[k]); // row_id, num
            }
            ++i;
            ++k;
        }
    }
 
    IT row;
    IVT num;
    if(sHeap.getSize() > 0)
    {
        sHeap.deleteMin(&row, &num);
        indy.push_back( (int32_t) row);
        numy.push_back( num );
    }
    while(sHeap.getSize() > 0)
    {
        sHeap.deleteMin(&row, &num);
        if(indy.back() == row)
        {
            numy.back() = SR::add(numy.back(), num);
        }
        else
        {
            indy.push_back( (int32_t) row);
            numy.push_back(num);
        }
    }       
}
 

template <typename SR, typename IT, typename IVT, class OVT>
void SpImpl<SR,IT,bool,IVT,OVT>::SpMXSpV(const Dcsc<IT,bool> & Adcsc, int32_t mA, const int32_t * indx, const IVT * numx, int32_t veclen,  
            int32_t * indy, OVT * numy, int * cnts, int * dspls, int p_c)
{   
    OVT * localy = new OVT[mA];
    BitMap isthere(mA);
    std::vector< std::vector<int32_t> > nzinds(p_c);    // nonzero indices      
 
    int32_t perproc = mA / p_c; 
    int32_t k = 0;  // index to indx vector
    IT i = 0;   // index to columns of matrix
    while(i< Adcsc.nzc && k < veclen)
    {
        if(Adcsc.jc[i] < indx[k]) ++i;
        else if(indx[k] < Adcsc.jc[i]) ++k;
        else
        {
            for(IT j=Adcsc.cp[i]; j < Adcsc.cp[i+1]; ++j)   // for all nonzeros in this column
            {
                int32_t rowid = (int32_t) Adcsc.ir[j];
                if(!isthere.get_bit(rowid))
                {
                    int32_t owner = std::min(rowid / perproc, static_cast<int32_t>(p_c-1));             
                    localy[rowid] = numx[k];    // initial assignment, requires implicit conversion if IVT != OVT
                    nzinds[owner].push_back(rowid);
                    isthere.set_bit(rowid);
                }
                else    
                {
                    localy[rowid] = SR::add(localy[rowid], numx[k]);
                }   
            }
            ++i;
            ++k;
        }
    }
 
    for(int p = 0; p< p_c; ++p)
    {
        sort(nzinds[p].begin(), nzinds[p].end());
        cnts[p] = nzinds[p].size();
        int32_t * locnzinds = &nzinds[p][0];
        int32_t offset = perproc * p;
        for(int i=0; i< cnts[p]; ++i)
        {
            indy[dspls[p]+i] = locnzinds[i] - offset;   // convert to local offset
            numy[dspls[p]+i] = localy[locnzinds[i]];    
        }
    }
    delete [] localy;
}
 
 
// this version is still very good with splitters
template <typename SR, typename IT, typename IVT, typename OVT>
void SpImpl<SR,IT,bool,IVT,OVT>::SpMXSpV_ForThreading(const Dcsc<IT,bool> & Adcsc, int32_t mA, const int32_t * indx, const IVT * numx, int32_t veclen,
                                                      std::vector<int32_t> & indy, std::vector<OVT> & numy, int32_t offset)
{
    std::vector<OVT> localy(mA);
    BitMap isthere(mA);
    std::vector<uint32_t> nzinds;   // nonzero indices
    
    SpMXSpV_ForThreading(Adcsc, mA, indx, numx, veclen, indy, numy, offset, localy, isthere, nzinds);
}
 
 

template <typename SR, typename IT, typename IVT, typename OVT>
void SpImpl<SR,IT,bool,IVT,OVT>::SpMXSpV_ForThreading(const Dcsc<IT,bool> & Adcsc, int32_t mA, const int32_t * indx, const IVT * numx, int32_t veclen, std::vector<int32_t> & indy, std::vector<OVT> & numy, int32_t offset, std::vector<OVT> & localy, BitMap & isthere, std::vector<uint32_t> & nzinds)
{
    // The following piece of code is not general, but it's more memory efficient than FillColInds
    int32_t k = 0;  // index to indx vector
    IT i = 0;   // index to columns of matrix
    while(i< Adcsc.nzc && k < veclen)
    {
        if(Adcsc.jc[i] < indx[k]) ++i;
        else if(indx[k] < Adcsc.jc[i]) ++k;
        else
        {
            for(IT j=Adcsc.cp[i]; j < Adcsc.cp[i+1]; ++j)   // for all nonzeros in this column
            {
                uint32_t rowid = (uint32_t) Adcsc.ir[j];
                if(!isthere.get_bit(rowid))
                {
                    localy[rowid] = numx[k];    // initial assignment
                    nzinds.push_back(rowid);
                    isthere.set_bit(rowid);
                }
                else
                {
                    localy[rowid] = SR::add(localy[rowid], numx[k]);
                }   
            }
            ++i; ++k;
        }
    }
    int nnzy = nzinds.size();
    integerSort(nzinds.data(), nnzy);
    indy.resize(nnzy);
    numy.resize(nnzy);
    for(int i=0; i< nnzy; ++i)
    {
        indy[i] = nzinds[i] + offset;   // return column-global index and let gespmv determine the receiver's local index
        numy[i] = localy[nzinds[i]];    
    }
}
 
 
 
 

 
template <typename SR, typename IT, typename NT, typename IVT, typename OVT>
void SpMXSpV_HeapSort(const Csc<IT,NT> & Acsc, int32_t mA, const int32_t * indx, const IVT * numx, int32_t veclen, std::vector<int32_t> & indy, std::vector<OVT> & numy, int32_t offset)
{
    IT inf = std::numeric_limits<IT>::min();
    IT sup = std::numeric_limits<IT>::max();
    KNHeap< IT, OVT > sHeap(sup, inf);
    
    
    for (int32_t k = 0; k < veclen; ++k)
    {
        IT colid = indx[k];
        for(IT j=Acsc.jc[colid]; j < Acsc.jc[colid+1]; ++j)
        {
            OVT val = SR::multiply( Acsc.num[j], numx[k]);
            sHeap.insert(Acsc.ir[j], val);
        }
    }
    
    IT row;
    OVT num;
    if(sHeap.getSize() > 0)
    {
        sHeap.deleteMin(&row, &num);
        row += offset;
        indy.push_back( (int32_t) row);
        numy.push_back( num );
    }
    while(sHeap.getSize() > 0)
    {
        sHeap.deleteMin(&row, &num);
        row += offset;
        if(indy.back() == row)
        {
            numy.back() = SR::add(numy.back(), num);
        }
        else
        {
            indy.push_back( (int32_t) row);
            numy.push_back(num);
        }
    }
}
 
 
 
template <typename SR, typename IT, typename NT, typename IVT, typename OVT>
void SpMXSpV_Bucket(const Csc<IT,NT> & Acsc, int32_t mA, const int32_t * indx, const IVT * numx, int32_t veclen,
                         std::vector<int32_t> & indy, std::vector< OVT > & numy, PreAllocatedSPA<OVT> & SPA)
{
    if(veclen==0)
        return;
    
    double tstart = MPI_Wtime();
    int nthreads=1;
    int rowSplits = SPA.buckets; // SPA must be initialized as checked in SpImpl.h
#ifdef _OPENMP
#pragma omp parallel
    {
        nthreads = omp_get_num_threads();
    }
#endif
    if(rowSplits < nthreads)
    {
        std::ostringstream outs;
        outs << "Warning in SpMXSpV_Bucket: " << rowSplits << " buckets are supplied for " << nthreads << " threads\n";
        outs << "4 times the number of threads are recommended when creating PreAllocatedSPA\n";
        SpParHelper::Print(outs.str());
    }
    
    int32_t rowPerSplit = mA / rowSplits;
    
    
    //------------------------------------------------------
    // Step1: count the nnz in each rowsplit of the matrix,
    // because we don't want to waste memory
    // False sharing is not a big problem because it is written outside of the main loop
    //------------------------------------------------------
    
    std::vector<std::vector<int32_t>> bSize(rowSplits, std::vector<int32_t> ( rowSplits, 0));
    std::vector<std::vector<int32_t>> bOffset(rowSplits, std::vector<int32_t> ( rowSplits, 0));
    std::vector<int32_t> sendSize(rowSplits);
    double t0, t1, t2, t3, t4;
#ifdef BENCHMARK_SPMSPV
    t0 = MPI_Wtime();
#endif
    
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 1)
#endif
    for(int b=0; b<rowSplits; b++)
    {
        // evenly balance nnz of x among threads
        int perBucket = veclen/rowSplits;
        int spill = veclen%rowSplits;
        int32_t xstart = b*perBucket + std::min(spill, b);
        int32_t xend = (b+1)*perBucket + std::min(spill, b+1);
        std::vector<int32_t> temp(rowSplits,0);
        for (int32_t i = xstart; i < xend; ++i)
        {
            IT colid = indx[i];
            for(IT j=Acsc.jc[colid]; j < Acsc.jc[colid+1]; ++j)
            {
                uint32_t rowid = (uint32_t) Acsc.ir[j];
                int32_t splitId = rowSplits-1;
                if(rowPerSplit!=0) splitId = (rowid/rowPerSplit > rowSplits-1) ? rowSplits-1 : rowid/rowPerSplit;
                //bSize[b][splitId]++;
                temp[splitId]++;
            }
        }
        int32_t totSend = 0;
        for(int k=0; k<rowSplits; k++)
        {
            bSize[b][k] = temp[k];
            totSend += temp[k];
        }
        sendSize[b] = totSend;
    }
    
    
#ifdef BENCHMARK_SPMSPV
    t1 = MPI_Wtime() - t0;
    t0 = MPI_Wtime();
#endif
    
    
    
    // keep it sequential to avoid fault sharing
    for(int i=1; i<rowSplits; i++)
    {
        for(int j=0; j<rowSplits; j++)
        {
            bOffset[i][j] = bOffset[i-1][j] + bSize[i-1][j];
            bSize[i-1][j] = 0;
        }
    }
    
    std::vector<uint32_t> disp(rowSplits+1);
    int maxBucketSize = -1; // maximum size of a bucket
    disp[0] = 0;
    for(int j=0; j<rowSplits; j++)
    {
        int thisBucketSize = bOffset[rowSplits-1][j] + bSize[rowSplits-1][j];
        disp[j+1] = disp[j] + thisBucketSize;
        bSize[rowSplits-1][j] = 0;
        maxBucketSize = std::max(thisBucketSize, maxBucketSize);
    }
    
    
    
#ifdef BENCHMARK_SPMSPV
    double  tseq = MPI_Wtime() - t0;
#endif
    //------------------------------------------------------
    // Step2: The matrix is traversed column by column and
    // nonzeros each rowsplit of the matrix are compiled together
    //------------------------------------------------------
    // Thread private buckets should fit in L2 cache
#ifndef L2_CACHE_SIZE
#define L2_CACHE_SIZE 256000
#endif
    int THREAD_BUF_LEN = 256;
    int itemsize = sizeof(int32_t) + sizeof(OVT);
    while(true)
    {
        int bufferMem = THREAD_BUF_LEN * rowSplits * itemsize + 8 * rowSplits;
        if(bufferMem>L2_CACHE_SIZE ) THREAD_BUF_LEN/=2;
        else break;
    }
    THREAD_BUF_LEN = std::min(maxBucketSize+1,THREAD_BUF_LEN);
    
#ifdef BENCHMARK_SPMSPV
    t0 = MPI_Wtime();
#endif
    
#ifdef _OPENMP
#pragma omp parallel
#endif
    {
        int32_t* tIndSplitA = new int32_t[rowSplits*THREAD_BUF_LEN];
        OVT* tNumSplitA = new OVT[rowSplits*THREAD_BUF_LEN];
        std::vector<int32_t> tBucketSize(rowSplits);
        std::vector<int32_t> tOffset(rowSplits);
#ifdef _OPENMP
#pragma omp for schedule(dynamic,1)
#endif
        for(int b=0; b<rowSplits; b++)
        {
            
            std::fill(tBucketSize.begin(), tBucketSize.end(), 0);
            std::fill(tOffset.begin(), tOffset.end(), 0);
            int perBucket = veclen/rowSplits;
            int spill = veclen%rowSplits;
            int32_t xstart = b*perBucket + std::min(spill, b);
            int32_t xend = (b+1)*perBucket + std::min(spill, b+1);
            
            for (int32_t i = xstart; i < xend; ++i)
            {
                IT colid = indx[i];
                for(IT j=Acsc.jc[colid]; j < Acsc.jc[colid+1]; ++j)
                {
                    OVT val = SR::multiply( Acsc.num[j], numx[i]);
                    uint32_t rowid = (uint32_t) Acsc.ir[j];
                    int32_t splitId = rowSplits-1;
                    if(rowPerSplit!=0) splitId = (rowid/rowPerSplit > rowSplits-1) ? rowSplits-1 : rowid/rowPerSplit;
                    if (tBucketSize[splitId] < THREAD_BUF_LEN)
                    {
                        tIndSplitA[splitId*THREAD_BUF_LEN + tBucketSize[splitId]] = rowid;
                        tNumSplitA[splitId*THREAD_BUF_LEN  + tBucketSize[splitId]++] = val;
                    }
                    else
                    {
                        std::copy(tIndSplitA + splitId*THREAD_BUF_LEN, tIndSplitA + (splitId+1)*THREAD_BUF_LEN, &SPA.indSplitA[disp[splitId] + bOffset[b][splitId]] + tOffset[splitId]);
                        std::copy(tNumSplitA + splitId*THREAD_BUF_LEN, tNumSplitA + (splitId+1)*THREAD_BUF_LEN, &SPA.numSplitA[disp[splitId] + bOffset[b][splitId]] + tOffset[splitId]);
                        tIndSplitA[splitId*THREAD_BUF_LEN] = rowid;
                        tNumSplitA[splitId*THREAD_BUF_LEN] = val;
                        tOffset[splitId] += THREAD_BUF_LEN ;
                        tBucketSize[splitId] = 1;
                    }
                }
            }
            
            for(int splitId=0; splitId<rowSplits; ++splitId)
            {
                if(tBucketSize[splitId]>0)
                {
                    std::copy(tIndSplitA + splitId*THREAD_BUF_LEN, tIndSplitA + splitId*THREAD_BUF_LEN + tBucketSize[splitId], &SPA.indSplitA[disp[splitId] + bOffset[b][splitId]] + tOffset[splitId]);
                    std::copy(tNumSplitA + splitId*THREAD_BUF_LEN, tNumSplitA + splitId*THREAD_BUF_LEN + tBucketSize[splitId], &SPA.numSplitA[disp[splitId] + bOffset[b][splitId]] + tOffset[splitId]);
                }
            }
        }
        delete [] tIndSplitA;
        delete [] tNumSplitA;
    }
    
#ifdef BENCHMARK_SPMSPV
    t2 = MPI_Wtime() - t0;
    t0 = MPI_Wtime();
#endif
    std::vector<uint32_t> nzInRowSplits(rowSplits);
    uint32_t* nzinds = new uint32_t[disp[rowSplits]];
    
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic,1)
#endif
    for(int rs=0; rs<rowSplits; ++rs)
    {
        
        for(int i=disp[rs]; i<disp[rs+1] ; i++)
        {
            int32_t lrowid = SPA.indSplitA[i] - rs * rowPerSplit;
            SPA.V_isthereBool[rs][lrowid] = false;
        }
        uint32_t tMergeDisp = disp[rs];
        for(int i=disp[rs]; i<disp[rs+1] ; i++)
        {
            int32_t rowid = SPA.indSplitA[i];
            int32_t lrowid = rowid - rs * rowPerSplit;
            if(!SPA.V_isthereBool[rs][lrowid])// there is no conflict across threads
            {
                SPA.V_localy[0][rowid] = SPA.numSplitA[i];
                nzinds[tMergeDisp++] = rowid;
                SPA.V_isthereBool[rs][lrowid]=true;
            }
            else
            {
                SPA.V_localy[0][rowid] = SR::add(SPA.V_localy[0][rowid], SPA.numSplitA[i]);
            }
        }
        
        integerSort(nzinds + disp[rs], tMergeDisp - disp[rs]);
        nzInRowSplits[rs] = tMergeDisp - disp[rs];
        
    }
    
#ifdef BENCHMARK_SPMSPV
    t3 = MPI_Wtime() - t0;
#endif
    // prefix sum
    std::vector<uint32_t> dispRowSplits(rowSplits+1);
    dispRowSplits[0] = 0;
    for(int i=0; i<rowSplits; i++)
    {
        dispRowSplits[i+1] = dispRowSplits[i] + nzInRowSplits[i];
    }
    
#ifdef BENCHMARK_SPMSPV
    t0 = MPI_Wtime();
#endif
    int nnzy = dispRowSplits[rowSplits];
    indy.resize(nnzy);
    numy.resize(nnzy);
#ifdef BENCHMARK_SPMSPV
    tseq = MPI_Wtime() - t0;
    t0 = MPI_Wtime();
#endif
    
    int  maxNnzInSplit = *std::max_element(nzInRowSplits.begin(),nzInRowSplits.end());
    THREAD_BUF_LEN = std::min(maxNnzInSplit+1,256);
#ifdef _OPENMP
#pragma omp parallel
#endif
    {
        OVT* tnumy = new OVT [THREAD_BUF_LEN];
        int32_t* tindy = new int32_t [THREAD_BUF_LEN];
        int curSize, tdisp;
#ifdef _OPENMP
#pragma omp for schedule(dynamic,1)
#endif
        for(int rs=0; rs<rowSplits; rs++)
        {
            curSize = 0;
            tdisp = 0;
            uint32_t * thisind = nzinds + disp[rs];
            std::copy(nzinds+disp[rs], nzinds+disp[rs]+nzInRowSplits[rs], indy.begin()+dispRowSplits[rs]);
            for(int j=0; j<nzInRowSplits[rs]; j++)
            {
                
                if ( curSize < THREAD_BUF_LEN)
                {
                    tnumy[curSize++] = SPA.V_localy[0][thisind[j]];
                }
                else
                {
                    std::copy(tnumy, tnumy+curSize, numy.begin()+dispRowSplits[rs]+tdisp);
                    tdisp += curSize;
                    tnumy[0] = SPA.V_localy[0][thisind[j]];
                    curSize = 1;
                }
            }
            if ( curSize > 0)
            {
                std::copy(tnumy, tnumy+curSize, numy.begin()+dispRowSplits[rs]+tdisp);
            }
        }
        delete [] tnumy;
        delete [] tindy;
    }
    
    
#ifdef BENCHMARK_SPMSPV
    t4 = MPI_Wtime() - t0;
#endif
    
    delete[] nzinds;
    
    
    
#ifdef BENCHMARK_SPMSPV
    double tall = MPI_Wtime() - tstart;
    std::ostringstream outs1;
    outs1 << "Time breakdown of SpMSpV-bucket." << std::endl;
    outs1 << "Estimate buckets: "<< t1 << " Bucketing: " << t2 << " SPA-merge: " << t3 << " Output: " << t4  << " Total: "<< tall << std::endl;
    SpParHelper::Print(outs1.str());
#endif
    
}
 
}