src.main/CjMatrix.hh

/* -*- Mode: C; indent-tabs-mode: t; c-basic-offset: 4; tab-width: 4 -*- */
/*
 * mitie
 * Copyright (C) Jason Spencer 2009 <mitie@jasonspencer.org>
 * 
 * mitie 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.
 * 
 * mitie 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 this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef _CJMATRIX_H_
#define _CJMATRIX_H_ 1

#include "CjMatrixException.hh"
#include "CjMatrixImpl.hh"
#include "CjMatrixImpl_1D.hh"

#include <iterator>

#include <ostream> // for dumpMatrix
#include <iomanip>      // for setw

#include <cassert>
#include <typeinfo> // for use in CjMatrixSubVector_const_iterator::operator<<(ostream..)
// T must be default constructable

//fwd decl
template <typename, bool> class CjMatrixSubVector_const_iterator;
template <typename, bool> class CjMatrixSubVector_iterator;
template <typename, bool> class CjMatrixSubVector;
template <typename, bool> class CjMatrixSubVectorC;


// JSPNOTE: it might be a good idea to include the matrix dimensions in the template which would prevent attempts of binary operations (+*..) on matrices which are of different dimeneions at compile time. XbyY * YbyX can be handled by free friend functions. Loading unknown sized matrices from a file/stream may be a little harder.
template <typename T> class CjMatrix {

private:
CjMatrixImpl<T> * m_impl;

protected:
CjMatrixImpl<T> * get_impl() { return m_impl; }

public:


CjMatrix (const unsigned int width, const unsigned int height, const T * const defaultValue = 0):m_impl( new CjMatrixImpl_1D<T>( width, height, defaultValue ) ) { }


CjMatrix (CjMatrixImpl<T> * impl) throw ():m_impl( impl ) { }   // the ownership and responsibility for destruction of impl is taken over by CjMatrix


explicit  CjMatrix (const CjMatrix & cjm, const unsigned int width, const unsigned int height, const T * const defaultValue = 0 ):m_impl( new CjMatrixImpl_1D<T>( cjm, width, height, defaultValue ) ) { };


explicit CjMatrix (const CjMatrix & cjm):m_impl(cjm.m_impl->clone()) { }

~CjMatrix() { delete m_impl; }


CjMatrix<T> & operator=(const CjMatrix<T> & cjm) {
/*      // TODO this is a simple solution, but maybe it would be easier to resize current impl and copy values? No chance of a throw if a re-alloc isn't required... and it'll maintain the impl type
        CjMatrixImpl<T> * new_impl = cjm.m_impl->clone();
        std::swap(new_impl,m_impl);
        delete new_impl;
        return *this;*/
        if (this != &cjm) {
                m_impl->resize_with_copyImpl(cjm.getWidth(), cjm.getHeight(), &cjm);
        }
        return *this;
        }


void resize_with_copy(unsigned int width, unsigned int height, const CjMatrix<T> * cjm = (CjMatrix<T> *) 0) { m_impl->resize_with_copyImpl(width, height, cjm); };

// these next four functions could (should?) be free functions with STL copy/fill to CjMatrix iterators, but this is potentially faster and leads to easier to read code

CjMatrix<T> & setColumnValues(unsigned int col, const T & val) { m_impl->setColumnValuesImpl(col, val); return *this; }

CjMatrix<T> & setRowValues(unsigned int row, const T & val) { m_impl->setRowValuesImpl(row, val); return *this; }

CjMatrix<T> & setAllValues(const T val) { m_impl->setAllValuesImpl(val); return *this; }

CjMatrix<T> & setDiagonal(const T & val) { m_impl->setDiagonalImpl(val); return *this; }


CjMatrix<T> & setDefaultValue(const T & val) { m_impl->setDefaultValueImpl( val ); return *this; }
T getDefaultValue() const { return m_impl->getDefaultValueImpl(); }

bool getDefaultValueSet() const { return m_impl->getDefaultValueSetImpl(); }


bool operator==(const CjMatrix<T> & cjm) const { return m_impl->equalElements(cjm); }

bool operator!=(const CjMatrix<T> & cjm) const { return !operator==(cjm); }

CjMatrix<T> & transpose(void) { m_impl->transposeImpl(); return *this; }        // transpose the matrix


bool isSymmetric() const { return m_impl->isSymmetricImpl(); }
unsigned int getWidth() const { return m_impl->getWidthImpl(); }
unsigned int getHeight() const { return m_impl->getHeightImpl(); }


inline T & at (const unsigned int row, const unsigned int col) { return m_impl->atImpl(row,col); }

inline const T & at (const unsigned int row, const unsigned int col) const { return m_impl->atImpl(row,col); }

inline T & operator() (const unsigned int row, const unsigned int col) { return m_impl->getImpl(row,col); }

inline const T & operator() (const unsigned int row, const unsigned int col) const { return m_impl->getImpl(row,col); }
 
// subvectors..
typedef CjMatrixSubVector<T, false>     row_t;
typedef CjMatrixSubVector<T, true>              column_t;
typedef CjMatrixSubVectorC<T, false>    const_row_t;
typedef CjMatrixSubVectorC<T, true>     const_column_t;

inline row_t row(const int r) { return row_t( r, this); }
inline column_t column(const int c) { return column_t( c, this); }
inline const_row_t row(const int r) const { return const_row_t( r, this); }
inline const_column_t column(const int c) const { return const_column_t( c, this); }

typename column_t::iterator rows_begin() { return column(0).begin(); }
typename column_t::iterator rows_end() { return column(0).end(); }
typename const_column_t::const_iterator rows_cbegin() { return column(0).cbegin(); }
typename const_column_t::const_iterator rows_cend() { return column(0).cend(); }
typename row_t::iterator columns_begin() { return row(0).begin(); }
typename row_t::iterator columns_end() { return row(0).end(); }
typename const_row_t::const_iterator columns_cbegin() { return row(0).cbegin(); }
typename const_row_t::const_iterator columns_cend() { return row(0).cend(); }

// friends..
//friend void std::swap<> ( CjMatrix<T> & a, CjMatrix<T> & b ); //nothrow
// romans..
// countrymen..
void swap (CjMatrix<T> & other) { std::swap(m_impl, other.m_impl); }
};

template<typename T, bool is_col> class CjMatrixSubVector {
protected:      
        CjMatrix <T> * const m_mat;
        unsigned int m_row_col;         // this vector represents this row/column [0..n-1]
public:
        typedef CjMatrixSubVector_iterator<T, is_col>                   iterator;
        typedef CjMatrixSubVector_const_iterator<T, is_col>     const_iterator;
    
        CjMatrixSubVector( const unsigned int row_col, CjMatrix <T> * const mat):m_mat(mat),m_row_col(row_col) { }
        CjMatrixSubVector( const CjMatrixSubVector<T, is_col> & orig):m_mat(orig.m_mat),m_row_col(orig.m_row_col) { }
        // the default shallow copy constructor is fine.

        iterator begin() { return iterator(m_mat, m_row_col, false); }
        iterator end() { return iterator(m_mat, m_row_col, true); }
        const_iterator cbegin() { return const_iterator(m_mat, m_row_col, false); }
        const_iterator cend() { return const_iterator(m_mat, m_row_col, true); }
        
        template<typename TP, bool is_colP> friend class CjMatrixSubVectorC;
};

template<typename T, bool is_col> class CjMatrixSubVectorC {
protected:
        const CjMatrix <T> * const m_mat;
        unsigned int m_row_col;
public:
        typedef CjMatrixSubVector_const_iterator<T, is_col>     const_iterator;

        CjMatrixSubVectorC (const unsigned int row_col, const CjMatrix <T> * const mat):m_mat(mat),m_row_col(row_col) { }
        CjMatrixSubVectorC( const CjMatrixSubVectorC<T, is_col> & orig):m_mat(orig.m_mat),m_row_col(orig.m_row_col) { }
        CjMatrixSubVectorC( const CjMatrixSubVector<T, is_col> & orig):m_mat(orig.m_mat),m_row_col(orig.m_row_col) { }

        const_iterator cbegin() { return const_iterator(m_mat, m_row_col, false); }
        const_iterator cend() { return const_iterator(m_mat, m_row_col, true); }
};

template <typename U, bool is_col> class CjMatrixSubVector_iterator : public std::iterator < std::random_access_iterator_tag, U, int > {
protected:
        CjMatrix<U> * m_mat;
        int m_n;                // the offset into the row or column - whether reverse or not, it's always the offset from start - if it goes beyond the bounds then noone cares.
        int m_row_col;  // the row or column number that we are representing

public:
        CjMatrixSubVector_iterator( CjMatrix <U> * const mat, const unsigned int row_col, const bool is_end ):
                m_mat(mat),
                m_n( is_end ? (is_col?mat->getHeight():mat->getWidth()) : 0 ),
                m_row_col(row_col) { }

        template <typename V, bool is_col_prime> friend class CjMatrixSubVector_iterator;       // to enable easy conversion with !is_col type
        template <typename V, bool is_col_prime> friend class CjMatrixSubVector_const_iterator; // to allow operator- on non-const type

        CjMatrixSubVector_iterator( const CjMatrixSubVector_iterator<U, is_col> & i ): m_mat(i.m_mat), m_n( i.m_n ), m_row_col(i.m_row_col) { }
        explicit CjMatrixSubVector_iterator( const CjMatrixSubVector_iterator<U, !is_col> & i ): m_mat(i.m_mat), m_n( i.m_row_col ), m_row_col(i.m_n) { }
        CjMatrixSubVector_iterator<U,is_col> & operator=(const CjMatrixSubVector_iterator<U,is_col> & it) {     m_mat = it.m_mat; m_n=it.m_n; m_row_col=it.m_row_col; return *this; }

        CjMatrixSubVector_iterator<U,!is_col> convertType() const { return CjMatrixSubVector_iterator<U,!is_col>(*this); }
        bool atEnd() const { return m_n == (is_col?(int)m_mat->getHeight():(int)m_mat->getWidth()); }   // TODO should be >= rather than == ?    nope - test is whether it is end iterator, not whether it is past end
        CjMatrixSubVector_iterator<U,is_col> end_iterator() { return CjMatrixSubVector_iterator<U,is_col>(m_mat, m_row_col, true); }
        CjMatrixSubVector_iterator<U,is_col> transpose( CjMatrix <U> * const mat) { if( (getRow() > (int)mat->getHeight()) || (getColumn() > (int)mat->getWidth()) ) throw E_CjMatrix_BoundsError("CjMatrixSubVector_iterator::transpose(CjMatrix): Invalid element location being transposed to new matrix"); CjMatrixSubVector_iterator tmp(*this); tmp.m_mat = mat; return tmp; }

        int getRowCol() const { return m_row_col; }
        int getOffset() const { return m_n; }
        int getRow() const { return (is_col?m_n:m_row_col); }
        int getColumn() const { return (is_col?m_row_col:m_n); }
        bool isSameMatrix(const CjMatrix <U> * const mat) const { return mat==m_mat; }
        bool atDiagonal() const { return m_n==m_row_col; }

        CjMatrixSubVector_iterator & operator++ () {  ++m_n; return *this; }
        CjMatrixSubVector_iterator operator++ (int) { CjMatrixSubVector_iterator<U, is_col> old_val(*this); ++(*this); return old_val; }
        CjMatrixSubVector_iterator & operator-- () {  --m_n; return *this; }
        CjMatrixSubVector_iterator operator-- (int) { CjMatrixSubVector_iterator<U, is_col> old_val(*this); --(*this); return old_val; }
        typename CjMatrixSubVector_iterator<U, is_col>::reference operator[] (const int n) { return m_mat->operator()( is_col?(m_n+n):m_row_col, is_col?m_row_col:(m_n+n) ); }

        bool operator==(const CjMatrixSubVector_iterator<U, true> & i) const { return ( (m_mat==i.m_mat) && ((!is_col)?(m_n==i.m_row_col):(m_n==i.m_n)) && ((!is_col)?(m_row_col==i.m_n):(m_row_col==i.m_row_col) )); }
        bool operator==(const CjMatrixSubVector_iterator<U, false> & i) const { return ( (m_mat==i.m_mat) && (is_col?(m_n==i.m_row_col):(m_n==i.m_n)) && (is_col?(m_row_col==i.m_n):(m_row_col==i.m_row_col) )); }
        bool operator==(const CjMatrixSubVector_const_iterator<U, true> & i) const { return ( (m_mat==i.m_mat) && ((!is_col)?(m_n==i.m_row_col):(m_n==i.m_n)) && ((!is_col)?(m_row_col==i.m_n):(m_row_col==i.m_row_col) )); }
        bool operator==(const CjMatrixSubVector_const_iterator<U, false> & i) const { return ( (m_mat==i.m_mat) && (is_col?(m_n==i.m_row_col):(m_n==i.m_n)) && (is_col?(m_row_col==i.m_n):(m_row_col==i.m_row_col) )); }

        bool operator!=(const CjMatrixSubVector_iterator<U, is_col> & i) const { return !(*this == i ); }
        bool operator!=(const CjMatrixSubVector_const_iterator<U, is_col> & i) const { return !(*this == i ); }

        typename CjMatrixSubVector_iterator<U, is_col>::reference operator* () { return m_mat->operator()( is_col?m_n:m_row_col, is_col?m_row_col:m_n ); }
        typename CjMatrixSubVector_iterator<U, is_col>::pointer operator->() const { return &(m_mat->operator()( is_col?m_n:m_row_col, is_col?m_row_col:m_n )); }
        CjMatrixSubVector_iterator<U, is_col> & operator+= (const int n) { m_n += n; return *this; }
        CjMatrixSubVector_iterator<U, is_col> & operator-= (const int n) { m_n -= n; return *this; }
        
        CjMatrixSubVector_iterator<U, is_col> operator+ (const int n) const { CjMatrixSubVector_iterator<U, is_col> ret(*this); ret.m_n = m_n+n; return ret; }
        CjMatrixSubVector_iterator<U, is_col> friend operator+ ( const int n, const CjMatrixSubVector_iterator<U, is_col> & rhs ) { CjMatrixSubVector_iterator<U, is_col> ret(rhs); ret.m_n += n; return ret; }
        CjMatrixSubVector_iterator<U, is_col> operator- (const int n) const { CjMatrixSubVector_iterator<U, is_col> ret(*this); ret.m_n = m_n-n; return ret; }
        typename CjMatrixSubVector_iterator<U, is_col>::difference_type operator- (const CjMatrixSubVector_iterator<U, is_col> & iter) const { return (m_n-iter.m_n); }
        typename CjMatrixSubVector_iterator<U, is_col>::difference_type operator- (const CjMatrixSubVector_const_iterator<U, is_col> & iter) const { return (m_n-iter.m_n); }

        bool operator< (const CjMatrixSubVector_iterator<U, is_col> & iter) const { return (m_n < iter.m_n); }
        bool operator> (const CjMatrixSubVector_iterator<U, is_col> & iter) const { return (m_n > iter.m_n); }
        bool operator<= (const CjMatrixSubVector_iterator<U, is_col> & iter) const { return (m_n <= iter.m_n); }
        bool operator>= (const CjMatrixSubVector_iterator<U, is_col> & iter) const { return (m_n >= iter.m_n); }
        bool operator< (const CjMatrixSubVector_const_iterator<U, is_col> & iter) const { return (m_n < iter.m_n); }
        bool operator> (const CjMatrixSubVector_const_iterator<U, is_col> & iter) const { return (m_n > iter.m_n); }
        bool operator<= (const CjMatrixSubVector_const_iterator<U, is_col> & iter) const { return (m_n <= iter.m_n); }
        bool operator>= (const CjMatrixSubVector_const_iterator<U, is_col> & iter) const { return (m_n >= iter.m_n); }
        
        friend std::ostream & operator<<(std::ostream & out, const CjMatrixSubVector_iterator<U, is_col> & a ) {
                out << "CjMatrixSubVector_iterator<" << typeid(U).name() << "," << (is_col?"true":"false") << ">{m_mat=" << a.m_mat << ",m_n=" << a.m_n << ",m_row_col=" << a.m_row_col << "}";
                return out;
                }
};

template <typename V, bool is_col> class CjMatrixSubVector_const_iterator : public std::iterator < std::random_access_iterator_tag, const V, int > {
protected:
        const CjMatrix<V> * m_mat;
        int m_n;                // the offset into the row or column - whether reverse or not, it's always the offset from start - if it goes beyond the bounds then noone cares.
        int m_row_col;  // the row or column number that we are representing - should be unsigned but it will be compared to signed ints and the code is more readable if the conversion is done in ctor rather than throughout the code
public:
        CjMatrixSubVector_const_iterator( const CjMatrix <V> * const mat, const unsigned int row_col, const bool is_end ):
                m_mat(mat),
                m_n( is_end ? (is_col?mat->getHeight():mat->getWidth()) : 0 ),
                m_row_col((int)row_col) { }

        template <typename U, bool is_col_prime> friend class CjMatrixSubVector_const_iterator; // to enable easy conversion with !is_col type
        template <typename U, bool is_col_prime> friend class CjMatrixSubVector_iterator;       // to allow operator- on non-const type
        // copy constructor
        CjMatrixSubVector_const_iterator( const CjMatrixSubVector_const_iterator<V, is_col> & i ): m_mat(i.m_mat), m_n( i.m_n ), m_row_col(i.m_row_col) { }
        // conversion from non-const constructor
        CjMatrixSubVector_const_iterator<V, is_col>( const CjMatrixSubVector_iterator<V, is_col> & i ): m_mat(i.m_mat), m_n( i.m_n ), m_row_col(i.m_row_col) { }
        explicit CjMatrixSubVector_const_iterator<V, is_col>( const CjMatrixSubVector_const_iterator<V, !is_col> & i ): m_mat(i.m_mat), m_n( i.m_row_col ), m_row_col(i.m_n) { }
        CjMatrixSubVector_const_iterator<V,is_col> & operator=(const CjMatrixSubVector_const_iterator<V,is_col> & it) { m_mat = it.m_mat; m_n=it.m_n; m_row_col=it.m_row_col; return *this; }
        
        int getRowCol() const { return m_row_col; }
        int getOffset() const { return m_n; }
        int getRow() const { return (is_col?m_n:m_row_col); }
        int getColumn() const { return (is_col?m_row_col:m_n); }
        bool isSameMatrix(const CjMatrix <V> * const mat) const { return mat==m_mat; }
        bool atDiagonal() const { return m_n==m_row_col; }

        // flips the direction from row traversal to column traversal, or vice versa
        CjMatrixSubVector_const_iterator<V,!is_col> convertType() const { return CjMatrixSubVector_const_iterator<V,!is_col>(*this); }
        bool atEnd() const { return m_n == (is_col?(int)m_mat->getHeight():(int)m_mat->getWidth()); } // TODO should be >= rather than == ?
        CjMatrixSubVector_const_iterator<V,is_col> end_iterator() { return CjMatrixSubVector_const_iterator<V,is_col>(m_mat, m_row_col, true); }
        CjMatrixSubVector_const_iterator<V,is_col> transpose(const CjMatrix <V> * const mat) { if( (getRow() > (int)mat->getHeight()) || (getColumn() > (int)mat->getWidth()) ) throw E_CjMatrix_BoundsError("CjMatrixSubVector_const_iterator::transpose(CjMatrix): Invalid element location being transposed to new matrix"); CjMatrixSubVector_const_iterator tmp(*this); tmp.m_mat = mat; return tmp; }
        
        CjMatrixSubVector_const_iterator & operator++ () {  ++m_n; return *this; }
        CjMatrixSubVector_const_iterator operator++ (int) { CjMatrixSubVector_const_iterator<V, is_col> old_val(*this); ++(*this); return old_val; }
        CjMatrixSubVector_const_iterator & operator-- () {  --m_n; return *this; }
        CjMatrixSubVector_const_iterator operator-- (int) { CjMatrixSubVector_const_iterator<V, is_col> old_val(*this); --(*this); return old_val; }
        const typename CjMatrixSubVector_const_iterator<V, is_col>::reference operator[] (const int n) { return m_mat->operator()( is_col?(m_n+n):m_row_col, is_col?m_row_col:(m_n+n) ); }

        bool operator==(const CjMatrixSubVector_const_iterator<V, true> & i) const { return ( (m_mat==i.m_mat) && ((!is_col)?(m_n==i.m_row_col):(m_n==i.m_n)) && ((!is_col)?(m_row_col==i.m_n):(m_row_col==i.m_row_col) )); }
        bool operator==(const CjMatrixSubVector_const_iterator<V, false> & i) const { return ( (m_mat==i.m_mat) && (is_col?(m_n==i.m_row_col):(m_n==i.m_n)) && (is_col?(m_row_col==i.m_n):(m_row_col==i.m_row_col) )); }
        bool operator==(const CjMatrixSubVector_iterator<V, true> & i) const { return ( (m_mat==i.m_mat) && ((!is_col)?(m_n==i.m_row_col):(m_n==i.m_n)) && ((!is_col)?(m_row_col==i.m_n):(m_row_col==i.m_row_col) )); }
        bool operator==(const CjMatrixSubVector_iterator<V, false> & i) const { return ( (m_mat==i.m_mat) && (is_col?(m_n==i.m_row_col):(m_n==i.m_n)) && (is_col?(m_row_col==i.m_n):(m_row_col==i.m_row_col) )); }

        bool operator!=(const CjMatrixSubVector_const_iterator<V, is_col> & i) const { return !(*this == i ); }
        bool operator!=(const CjMatrixSubVector_iterator<V, is_col> & i) const { return !(*this == i ); }

        const typename CjMatrixSubVector_const_iterator<V, is_col>::reference operator* () { return m_mat->operator()( is_col?m_n:m_row_col, is_col?m_row_col:m_n ); }
        const typename CjMatrixSubVector_const_iterator<V, is_col>::pointer operator->() const { return &(m_mat->operator()( is_col?m_n:m_row_col, is_col?m_row_col:m_n )); }
        CjMatrixSubVector_const_iterator<V, is_col> & operator+= (const int n) { m_n += n; return *this; }
        CjMatrixSubVector_const_iterator<V, is_col> & operator-= (const int n) { m_n -= n; return *this; }
        
        CjMatrixSubVector_const_iterator<V, is_col> operator+ (const int n) const { CjMatrixSubVector_const_iterator<V, is_col> ret(*this); ret.m_n = m_n+n; return ret; }
        CjMatrixSubVector_const_iterator<V, is_col> friend operator+ ( const int n, const CjMatrixSubVector_const_iterator<V, is_col> & rhs ) { CjMatrixSubVector_const_iterator<V, is_col> ret(rhs); ret.m_n += n; return ret; }
        CjMatrixSubVector_const_iterator<V, is_col> operator- (const int n) const { CjMatrixSubVector_const_iterator<V, is_col> ret(*this); ret.m_n = m_n-n; return ret; }
        typename CjMatrixSubVector_const_iterator<V, is_col>::difference_type operator- (const CjMatrixSubVector_const_iterator<V, is_col> & iter) const { return (m_n-iter.m_n); }
        typename CjMatrixSubVector_const_iterator<V, is_col>::difference_type operator- (const CjMatrixSubVector_iterator<V, is_col> & iter) const { return (m_n-iter.m_n); }

        bool operator< (const CjMatrixSubVector_const_iterator<V, is_col> & iter) const { return (m_n < iter.m_n); }
        bool operator> (const CjMatrixSubVector_const_iterator<V, is_col> & iter) const { return (m_n > iter.m_n); }
        bool operator<= (const CjMatrixSubVector_const_iterator<V, is_col> & iter) const { return (m_n <= iter.m_n); }
        bool operator>= (const CjMatrixSubVector_const_iterator<V, is_col> & iter) const { return (m_n >= iter.m_n); }
        bool operator< (const CjMatrixSubVector_iterator<V, is_col> & iter) const { return (m_n < iter.m_n); }
        bool operator> (const CjMatrixSubVector_iterator<V, is_col> & iter) const { return (m_n > iter.m_n); }
        bool operator<= (const CjMatrixSubVector_iterator<V, is_col> & iter) const { return (m_n <= iter.m_n); }
        bool operator>= (const CjMatrixSubVector_iterator<V, is_col> & iter) const { return (m_n >= iter.m_n); }
        friend std::ostream & operator<<(std::ostream & out, const CjMatrixSubVector_const_iterator<V, is_col> & a ) {
                out << "CjMatrixSubVector_iterator<" << typeid(V).name() << "," << (is_col?"true":"false") << ">{m_mat=" << a.m_mat << ",m_n=" << a.m_n << ",m_row_col=" << a.m_row_col << "}";
                return out;
                }
};

// free functions

namespace std {
template <typename T> void swap (CjMatrix<T> & a, CjMatrix<T> & b) { a.swap(b); }
}

template <typename T> std::string dumpMatrix( const CjMatrix<T> & a, const unsigned int cellWidth, const std::string & rowPrefix, const bool blankDiagonal = false, const unsigned int max_col = 0, const unsigned int max_row = 0 ) {
std::stringstream ss(std::ios_base::out);
dumpMatrixTmpl(ss, a, cellWidth, rowPrefix, blankDiagonal, max_col, max_row);
return ss.str();
}

template <typename T> std::ostream & operator<<(std::ostream & out, const CjMatrix<T> & a ){
dumpMatrixTmpl(out, a, 0, "", false, 0, 0);
return out;
}

template <typename T> void dumpMatrixTmpl(std::ostream &os, const CjMatrix<T> & a, const unsigned int cellWidth, const std::string & rowPrefix, const bool blankDiagonal, const unsigned int max_col , const unsigned int max_row ) {

unsigned int h = (max_row==0)?a.getHeight():std::min(max_row,a.getHeight());
unsigned int w = (max_col==0)?a.getWidth():std::min(max_col,a.getWidth());
std::cout << "dumpMatrixTmpl: height is " << h << " and width is " << w << std::endl;
std::cout << "dumpMatrixTmpl: rowPrefix is " << rowPrefix << std::endl;
for( unsigned int i = 0; i < h; ++i) {
        os << std::setw(0);
        os << rowPrefix.c_str();
        os << std::right;
        typename CjMatrix<T>::const_row_t row = a.row(i);
        unsigned int j = 0;
        for( typename CjMatrix<T>::const_row_t::const_iterator it = row.cbegin(); (it != row.cend()) && (j!=w) ; ++it, ++j) {
                os << ((cellWidth==0)?"\t":"")  << std::setw(cellWidth) << ((blankDiagonal&&(i==j))?0:(*it));
        }
        if(i!=(h-1)) os << std::endl;
        }
}

#endif // _CJMATRIX_H_