SparseMatrix.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 Mathematics, Utrecht University, 2009.
 */


#ifndef _H_SM
#define _H_SM

#include "ULIdef.hpp"
#include "Matrix.hpp"
#include "FileToVT.hpp"
#include "alignment.hpp"

template< typename T, typename IND >
class SparseMatrix: public Matrix< T > {
        private:
        
        protected:

                IND nor;

                IND noc;

                IND nnz;

        public:

                T zero_element;

                SparseMatrix() {}

                SparseMatrix( const IND nzs, const IND nr, const IND nc, const T zero ):
                        nnz( nzs ), nor( nr ), noc( nc ), zero_element( zero ) {}

                virtual ~SparseMatrix() {}

                virtual void load( std::vector< Triplet< T > >& input, const IND m, const IND n, const T zero ) = 0;

                void loadFromFile( const std::string file, const T zero=0 ) {
                        ULI m, n;
                        const size_t pos = file.find_last_of( '.' );
                        const std::string ext = file.substr( pos + 1, file.length() );
                        std::vector< Triplet< T > > VT;
                        if( ext.compare( "trp" ) == 0 ) {
                                VT = Triplet< T >::load( file, m, n );
                        } else if( ext.compare( "crs" ) == 0 || ext.compare( "csr" ) == 0 ) {
                                VT = Triplet< T >::loadCRS( file, m, n );
                        } else //assume default matrix-market format
                                VT = FileToVT::parse( file, m, n );
                        this->load( VT, m, n, zero );
                }

                virtual unsigned long int m() {
                        return static_cast< unsigned long int >( nor );
                }

                virtual unsigned long int n() {
                        return static_cast< unsigned long int >( noc );
                }

                virtual unsigned long int nzs() {
                        return static_cast< unsigned long int >( nnz );
                }

                virtual void getFirstIndexPair( IND &row, IND &col ) = 0;

                virtual T* mv( const T* x ) {
                        T* ret = NULL;
                        //allocate output vector, catch error
                        if( posix_memalign( (void**)&ret, 64, nor * sizeof( T ) ) != 0 ) {
                                std::cerr << "SparseMatrix::mv(): could not allocate output vector; returning NULL..." << std::endl;
                                return NULL;
                        }
                        //initialise output vector
                        for( IND i=0; i<nor; i++ ) {
                                ret[ i ] = zero_element;
                        }
                        //call in-place z=Ax
                        zax( x, ret );
                        //return new output vector
                        return ret;
                }

                virtual void zax( const T*__restrict__ x, T*__restrict__ z ) = 0;

                virtual void zxa( const T*__restrict__ x, T*__restrict__ z ) = 0;

};

#endif //_H_SM