csbsym.h

Go to the documentation of this file.

#ifndef _CSBSYM_H
#define _CSBSYM_H

#include <cmath>
#include <vector>
#include <algorithm>
#include <numeric>      // for std:accumulate()
#include <limits>       // C++ style numeric_limits<T>
#include <ostream>
#include <iterator>
#include "csc.h"
#include "mortoncompare.h"

using namespace std;


inline void atomicallyIncrementDouble(volatile double *target, const double by){
  asm volatile(                                            
    "movq  %0, %%rax \n\t"          // rax = *(%0)
    "xorpd %%xmm0, %%xmm0 \n\t"         // xmm0 = [0.0,0.0]
    "movsd %1, %%xmm0\n\t"          // xmm0[lo] = *(%1)
    "1:\n\t"                 
    // rax (containing *target) was last set at startup or by a failed cmpxchg
    "movq  %%rax,  %%xmm1\n\t"          // xmm1[lo] = rax
    "addsd %%xmm0, %%xmm1\n\t"          // xmm1 = xmm0 + xmm1 = by + xmm1
    "movq  %%xmm1, %%r8  \n\t"          // r8 = xmm1[lo]
    "lock cmpxchgq %%r8, %0\n\t"        // if(*(%0)==rax){ZF=1;*(%0)=r8}else{ZF=0;rax=*(%0);}
    "jnz 1b\n\t"                                // jump back if failed (ZF=0)
    : "=m"(*target)                             // outputs
    : "m"(by)                           // inputs
    : "cc", "memory", "%rax", "%r8", "%xmm0", "%xmm1" // clobbered
  );
  return;
}

/* Symmetric CSB implementation
** Only upper triangle is stored
** top[i][0] gives the ith diagonal block for every i
** Since this class works only for symmetric (hence square) matrices,
** each compressed sparse block is (lowbits+1)x(lowbits+1) and ncsb = nbr = nbc
*/
template <class NT, class IT>
class CsbSym
{
public:
    CsbSym ():nz(0), n(0), ncsb(0) {}   // default constructor (dummy)

    CsbSym (const CsbSym<NT, IT> & rhs);            // copy constructor
    ~CsbSym();
    CsbSym<NT,IT> & operator=(const CsbSym<NT,IT> & rhs);   // assignment operator
    CsbSym (Csc<NT, IT> & csc, int workers);
    
    ofstream & PrintStats(ofstream & outfile) const;
    ofstream & Dump(ofstream & outfile) const;
    IT colsize() const { return n;} 
    IT rowsize() const { return n;} 
    bool isPar() const { return ispar; }

private:
    void Init(int workers, IT forcelogbeta = 0);
    void SeqSpMV(const NT * __restrict x, NT * __restrict y) const;
    void BMult(IT** chunks, IT start, IT end, const NT * __restrict x, NT * __restrict y, IT ysize) const;

    void BlockPar(IT start, IT end, const NT * __restrict subx, const NT * __restrict subx_mirror, 
            NT * __restrict suby, NT * __restrict suby_mirror, IT rangebeg, IT rangeend, IT cutoff) const;
    void BlockTriPar(IT start, IT end, const NT * __restrict subx, NT * __restrict suby, IT rangebeg, IT rangeend, IT cutoff) const;

    void SortBlocks(pair<IT, pair<IT,IT> > * pairarray, NT * val);
    void DivideIterationSpace(IT * & lspace, IT * & rspace, IT & lsize, IT & rsize, IT size, IT d) const;

    void MultAddAtomics(NT * __restrict y, const NT * __restrict x, const IT d) const;
    void MultDiag(NT * __restrict y, const NT * __restrict x, const IT d) const;
    void MultMainDiag(NT * __restrict y, const NT * __restrict x) const;

    float Imbalance(IT d) const;
    IT nzsum(IT d) const;

    IT ** top ; // pointers array (indexed by higher-order bits of the coordinate index), size = nbr*(nbc+1)
    IT * bot;   // contains lower-order bits of the coordinate index, size nnz
    NT * num;   // contains numerical values, size nnz

    vector< pair<IT,NT> > diagonal;

    bool ispar;
    IT nz;      // # nonzeros 
    IT n;       // #{rows} = #{columns}
    IT blcrange;    // range indexed by one block

    IT ncsb;    // #{block rows) = #{block cols}
    
    IT nlowbits;    // # lower order bits (for both rows and columns)
    IT nhighbits;
    IT highmask;    // mask with the first log(n)/2 bits = 1 and the other bits = 0 
    IT lowmask;

    MortCompSym<IT> mortoncmp;  // comparison operator w.r.t. the (inverted N)-morton layout

    template <typename NU, typename IU>
    friend void csbsym_gespmv (const CsbSym<NU, IU> & A, const NU * x, NU * y);
};


#include "friends.h"
#include "csbsym.cpp"
#endif