MKLCRS.hpp

/*
 * Copyright (c) 2007-2014, A. N. Yzelman,   Utrecht University 2007-2011;
 *                                                    KU Leuven 2011-2014.
 *                          R. H. Bisseling, Utrecht University 2007-2014.
 * 
 * This file is part of the Sparse Library.
 * 
 * This library was developed under supervision of Prof. dr. Rob H. Bisseling at
 * Utrecht University, from 2007 until 2011. From 2011-2014, development continued 
 * at KU Leuven, where Prof. dr. Dirk Roose contributed significantly to the ideas 
 * behind the newer parts of the library code.
 * 
 *     The Sparse Library is free software: you can redistribute it and/or modify
 *     it under the terms of the GNU General Public License as published by the
 *     Free Software Foundation, either version 3 of the License, or (at your
 *     option) any later version.
 * 
 *     The Sparse Library is distributed in the hope that it will be useful, but
 *     WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
 *     or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
 *     for more details.
 * 
 *     You should have received a copy of the GNU General Public License along
 *     with the Sparse Library. If not, see <http://www.gnu.org/licenses/>.
 */


/*
 * File created by:
 *     A. N. Yzelman, Dept. of Computer Science, KU Leuven, 2012.
 */


#include "CRS.hpp"

extern "C" {
 #include <mkl_service.h>
 #include <mkl_spblas.h>
}

//#define _DEBUG

#ifndef _H_MKLCRS
#define _H_MKLCRS

template< typename T >
class MKLCRS: public CRS< T > {

  private:

  protected:

        double _one;

        char trans;

        char *descr;

        int _m;

        int _n;

        int *_col_ind;

        int *_row_start;

        void prepare() {
                _m = static_cast< int >( this->nor );
                _n = static_cast< int >( this->noc );
                _one = 1.0;

                _col_ind = new int[ this->nnz ];
                for( ULI i=0; i<this->nnz; ++i )
                        _col_ind[ i ] = static_cast< int >( this->col_ind[ i ] );

                _row_start = new int[ _m + 1 ];
                for( int i=0; i<=_m; ++i )
                        _row_start[ i ] = static_cast< int >( this->row_start[ i ] );

                trans = 'n';

                descr = new char[4];
                descr[ 0 ] = 'G';
                descr[ 1 ] = descr[ 2 ] = '*';
                descr[ 3 ] = 'C';
        }

  public:

        MKLCRS() {}

        MKLCRS( std::string file, T zero = 0 ) {
                this->loadFromFile( file, zero );

                prepare();
        }
        
        MKLCRS( const ULI number_of_nonzeros, const ULI number_of_rows, const ULI number_of_cols, T zero ) {
                this->nnz = number_of_nonzeros;
                this->nor = number_of_rows;
                this->noc = number_of_cols;
                this->zero_element = zero;
                this->row_start = new ULI[ this->nor + 1 ];
                this->ds = new T[ this->nnz ];
                this->col_ind = new ULI[ this->nnz ];

                prepare();
        }

        MKLCRS( CRS< T >& toCopy ) {
                this->zero_element = toCopy.zero_element;
                this->nnz = toCopy.nnz;
                this->nor = toCopy.nor;
                this->row_start = new ULI[ this->nor + 1 ];
                this->ds = new T[ this->nnz ];
                this->col_ind = new ULI[ this->nnz ];
                for( ULI i=0; i<this->nnz; i++ ) {
                        this->ds[ i ] = toCopy.ds[ i ];
                        this->col_ind[ i ] = toCopy.col_ind[ i ];
                }
                for( ULI i=0; i<this->nor; i++ )
                        this->row_start[ i ] = toCopy.row_start[ i ];

                prepare();
        }
        
        MKLCRS( std::vector< Triplet< T > > input, ULI m, ULI n, T zero ) {
                load( input, m, n, zero );

                prepare();
        }

#ifndef _NO_LIBNUMA

        virtual T* mv( const T* x ) {
                T* ret = (T*) numa_alloc_interleaved( this->nor * sizeof( T ) );
                for( ULI i=0; i<this->nor; i++ )
                        ret[ i ] = this->zero_element;
                this->zax( x, ret );
                return ret;
        }
#endif

        virtual void zxa( const T*__restrict__ x, T*__restrict__ z ) {
                std::cerr << "MKLCRS does not implement the zxa, sorry!" << std::endl;
                exit( EXIT_FAILURE );
        }

        virtual void zax( const T*__restrict__ x, T*__restrict__ z ) {
                //Warning: assumes T is a double
                mkl_dcsrmv( &trans, &_m, &_n, &_one, descr, this->ds, _col_ind, _row_start, &( _row_start[ 1 ] ), (T*__restrict__)x, &_one, z );
        }

        virtual ~MKLCRS() {
                delete [] _col_ind;
                delete [] _row_start;
                delete [] descr;
        }

};

#endif