csc.h

Go to the documentation of this file.

#ifndef _CSC_H_
#define _CSC_H_

#include "triple.h"
#include <iterator>
#include <array>

using namespace std;


template <class T, class ITYPE>
struct Triple;

template <class T, class ITYPE>
class Csc
{
public:
    Csc ():nz(0), m(0), n(0), logicalnz(0), issym(false) {}             // default constructor
    Csc (ITYPE size,ITYPE rows, ITYPE cols, bool isSym=false);
    Csc (const Csc<T, ITYPE> & rhs);        // copy constructor
    ~Csc();
    Csc<T, ITYPE> & operator=(const Csc<T, ITYPE> & rhs);   // assignment operator
    Csc (Triple<T, ITYPE> * triples, ITYPE size, ITYPE rows, ITYPE cols, bool isSym=false);
    Csc (ITYPE * ri, ITYPE * ci, T * val, ITYPE size, ITYPE rows, ITYPE cols, bool isSym=false);

    // we have to use another function because the compiler will reject another constructor with the same signature
    void SetPointers (ITYPE * colpointers, ITYPE * rowindices, T * vals, ITYPE size, ITYPE rows, ITYPE cols, bool fortran)  
    {
        jc = colpointers;
        ir = rowindices;
        num = vals;
        nz = size;
        m = rows;
        n = cols;
        issym = false;
        logicalnz = size;

        if(fortran)
        {
            transform(jc, jc+n+1, jc, bind2nd(minus<ITYPE>(),1));
            transform(ir, ir+nz, ir, bind2nd(minus<ITYPE>(),1));
        }
    }

    ITYPE colsize() const { return n;} 
    ITYPE rowsize() const { return m;} 
    ITYPE * getjc() const { return jc;} 
    ITYPE * getir() const { return ir;} 
    T * getnum() const { return num;} 
    ITYPE getlogicalnnz() const
    {
        return logicalnz;
    }

private:
    void Resize(ITYPE nsize);
    bool issym;

    ITYPE * jc ;    // col pointers, size n+1
    ITYPE * ir;     // row indices, size nnz 
    T * num;    // numerical values, size nnz 
    
    ITYPE logicalnz;
    ITYPE nz;
    ITYPE m;        //  number of rows
    ITYPE n;        //  number of columns

    template <class U, class UTYPE>
    friend class CsbSym;
    template <class U, class UTYPE>
    friend class BiCsb;
    template <class U, class UTYPE, unsigned UDIM>
    friend class BmCsb;
    template <class U, class UTYPE, unsigned UDIM>
    friend class BmSym;

    template <typename U, typename UTYPE>
    friend void csc_gaxpy (const Csc<U, UTYPE> & A, U * x, U * y);

    template <typename U, typename UTYPE>
    friend void csc_gaxpy_trans (const Csc<U, UTYPE> & A, U * x, U * y);
    
    template <int D, typename NT, typename IT>
    friend void csc_gaxpy_mm(const Csc<NT,IT> & A, array<NT,D> * x, array<NT,D> * y);
    
    template <int D, typename NT, typename IT>
    friend void csc_gaxpy_mm_trans(const Csc<NT,IT> & A, array<NT,D> * x, array<NT,D> * y);
};

/* y = A*x+y */
template <typename T, typename ITYPE>
void csc_gaxpy (const Csc<T, ITYPE> & A, T * x, T * y)
{
    if(A.issym)
    {
        for (ITYPE j = 0 ; j < A.n ; ++j)   // for all columns of A
        {
            for (ITYPE k = A.jc [j] ; k < A.jc [j+1] ; ++k) 
            {
                y [ A.ir[k] ] += A.num[k]  * x [j] ;
                if( j != A.ir[k] )
                    y [ j ] += A.num[k] * x[ A.ir[k] ] ;    // perform the symmetric update
            }
        }
    }
    else
    {
        for (ITYPE j = 0 ; j < A.n ; ++j)   // for all columns of A
        {
            for (ITYPE k = A.jc [j] ; k < A.jc [j+1] ; ++k) // scale jth column with x[j]
            {
                y [ A.ir[k] ] += A.num[k]  * x [j] ;
            }
        }
    }
}


/* y = A' x + y */
template <typename T, typename ITYPE>
void csc_gaxpy_trans(const Csc<T,ITYPE> & A, T * x, T * y)
{
    if(A.issym)
    {
        cout << "Trying to run A'x on a symmetric matrix doesn't make sense" << endl;
        cout << "Are you sure you're using the right data structure?" << endl;
        return; 
    }
    
    for (ITYPE j = 0; j< A.n; ++j)
    {
        for(ITYPE k= A.jc[j]; k < A.jc[j+1]; ++k)
        {
            y[j] += A.num[k] * x [ A.ir[k] ]; 
        }
    }
}


/* Y = A X + Y */
template <int D, typename NT, typename IT>
void csc_gaxpy_mm(const Csc<NT,IT> & A, array<NT,D> * x, array<NT,D> * y)
{
    if(A.issym)
    {
        cout << "Symmetric csc_gaxpy_mm not implemented yet" << endl;
        return; 
    }
    
    for (IT j = 0 ; j < A.n ; ++j)  // for all columns of A
    {
        for (IT k = A.jc[j] ; k < A.jc[j+1] ; ++k)  // scale jth column with x[j]
        {
            for(int i=0; i<D; ++i)
            {
                y[A.ir[k]][i] += A.num[k] * x[j][i];
            }
        }
    }
    
}


/* Y = A' X + Y */
template <int D, typename NT, typename IT>
void csc_gaxpy_mm_trans(const Csc<NT,IT> & A, array<NT,D> * x, array<NT,D> * y)
{
    if(A.issym)
    {
        cout << "Trying to run A'x on a symmetric matrix doesn't make sense" << endl;
        cout << "Are you sure you're using the right data structure?" << endl;
        return; 
    }
    
    for (IT j = 0; j< A.n; ++j)
    {
        for(IT k= A.jc[j]; k < A.jc[j+1]; ++k)
        {
            for(int i=0; i<D; ++i)
            {
                y[j][i] +=  A.num[k] * x[A.ir[k]][i];
            }   
        }
    }
}


#include "csc.cpp"  // Template member function definitions need to be known to the compiler
#endif