DD_MATRIX.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.
 */


#include <assert.h>
#include <vector>
#include "Triplet.hpp"
#include "SparseMatrix.hpp"

template< typename T, int number_of_diagonals, int diagonal_offsets[] >
class DD_MATRIX: public SparseMatrix< T, unsigned long int > {

   private:
        
        typedef unsigned long int ULI;

   protected:

        T** nzs;

        ULI full;

        ULI d;

        bool SELF_ALLOCATED;

        size_t bytes;

   public:

        DD_MATRIX() {
                load( std::vector< Triplet< T > >(), 0, 0, 0 );
        }

        DD_MATRIX( std::string file, T zero = 0 ) {
                loadFromFile( file, zero );
        }
        
        DD_MATRIX( std::vector< Triplet< T > >& input, ULI m, ULI n, T zero ) {
                load( input, m, n, zero );
        }

        DD_MATRIX( T **nonzeroes, ULI m, ULI n, T zero ) {
                this->zero_element = zero;
                this->nor = m;
                this->noc = n;
                this->nnz = 0;

                for( ULI k=0; k<number_of_diagonals; k++ ) {
                        if( diagonal_offsets[ k ] < 0 ) {
                                this->nnz += m+diagonal_offsets[ k ] < n ? m+diagonal_offsets[ k ] : n;
                        } else {
                                this->nnz += n-diagonal_offsets[ k ] < m ? m-diagonal_offsets[ k ] : m;
                        }
                }

                d = m > n ? m : n;
                full = m > n ? m-n+1 : n-m+1;
                nzs = nonzeroes;

                bytes  = 2 * sizeof(ULI) + sizeof(T**);
                bytes += number_of_diagonals * sizeof(T*);
                bytes += number_of_diagonals * d * sizeof(T);

                SELF_ALLOCATED = false;
        }

        virtual void load( std::vector< Triplet< T > >& input, const ULI m, const ULI n, const T zero ) {

                this->zero_element = zero;
                this->nor = m;
                this->noc = n;
                this->nnz = input.size();

                nzs = new T*[ number_of_diagonals ];
                bytes = number_of_diagonals * sizeof(T*) + 2 * sizeof(ULI) + sizeof(T**);
                d = m > n ? m : n;
                full = m > n ? m-n+1 : n-m+1;

                //first look which diagonals are populated
                ULI *diags = new ULI[ m + n - 1 ];
                for( ULI r=0; r<m+n-1; r++ ) {
                        diags[ r ] = 0;
                }
                for( ULI r=0; r<this->nnz; r++ ) {
                        if( input[ r ].value != this->zero_element ) {
                                diags[input[ r ].j()-input[ r ].i()+m-1] = 1;
                        }
                }
                ULI nl1 = 0; //number of diagonals populated
                for( ULI r=0; r<m+n-1; r++ ) {
                        if( diags[ r ] == 1 ) {
                                nl1++;
                        }
                }
                
                //assume perfect human input
                ULI c = 0;
                assert( number_of_diagonals == nl1 );
                for( ULI r=0; r<m+n-1; r++ ) {
                        if( diags[ r ] == 1 ) {
                                assert( static_cast< int >( r ) - static_cast< int >( m ) + 1 == diagonal_offsets[ c++ ] );
                        }
                }

                //allocate nzs fully
                for( ULI r=0; r<number_of_diagonals; r++ ) {
                        ULI curdiag = diagonal_offsets[ r ];
                        const size_t allocLength = curdiag <= full ? d : d-(curdiag-full);
                        nzs[r] = new T[ allocLength ];
                        bytes += allocLength * sizeof(T);
                }
        
                //build diagonal code to nzs-row array
                for( ULI r=0; r<m+n-1; r++ ) {
                        diags[ r ] = static_cast< ULI >( -1 ); //signal overflow (hopefully)
                }
                for( ULI r=0; r<number_of_diagonals; r++ ) {
                        diags[ diagonal_offsets[ r ] + m - 1 ] = r;
                }
        
                for( ULI r=0; r<this->nnz; r++ ) {
                        if( input[ r ].value != this->zero_element ) {
                                ULI cur = input[ r ].j() - input[ r ].i() + m-1; //diagonal code
                                if( input[ r ].i() > input[ r ].j() )
                                        nzs[ diags[ cur ] ][ input[ r ].j() ] = input[ r ].value;
                                else
                                        nzs[ diags[ cur ] ][ input[ r ].i() ] = input[ r ].value;
                        }
                }

                //set self allocation flag
                SELF_ALLOCATED = true;

                //clean up
                delete [] diags;

                //done 
        }

        virtual void getFirstIndexPair( unsigned long int &row, unsigned long int &col ) {
                row = 0;
                col = 0;
        }

        virtual void zxa( const T* x, T* z ) {
                for( ULI j=0; j<d; j++ ) {
                        //theoretically, a compiler could unroll this inner loop perfectly
                        //and optimise out the if-statements.
                        for( ULI k=0; k<number_of_diagonals; k++ ) {
                                const ULI i = j + diagonal_offsets[ k ];
                                if( i >= this->nor ) continue;
                                if( diagonal_offsets[ k ] < 0 )
                                        z[ j ] += nzs[ k ][ j ] * x[ i ];
                                else
                                        z[ j ] += nzs[ k ][ i ] * x[ i ];
                        }
                }
        }

        virtual void zax( const T* x, T* z ) {
                //std::cout << "d= " << d << std::endl;
                for( ULI i=0; i<d; i++ ) {
                        //theoretically, a compiler could unroll this inner loop perfectly
                        //and optimise out the if-statements.
                        for( ULI k=0; k<number_of_diagonals; k++ ) {
                                const ULI j = i + diagonal_offsets[ k ];
                                //std::cout << "j>=noc: " << j << " >= " << this->noc << std::endl;
                                if( j >= this->noc ) continue;
                                //std::cout << "i,j,k = " << i << "," << j << "," << k << std::endl;
                                if( diagonal_offsets[ k ] < 0 )
                                        z[ i ] += nzs[ k ][ j ] * x[ j ];
                                else
                                        z[ i ] += nzs[ k ][ i ] * x[ j ];
                        }
                }
        }

        size_t bytesUsed() {
                return bytes;
        }

        ~DD_MATRIX() {
                if( !SELF_ALLOCATED ) return;

                for( ULI k = 0; k<number_of_diagonals; k++ )
                        delete [] nzs[ k ];
                delete [] nzs;
        }

};