crmatrix.h

Go to the documentation of this file.

//----------------------------------------------------------
// version 1.1, may 1998
//
// Marco Wiedenhaus
// Basis: TCRMatrixComplex by Jens Lenge
//
// $Id: crmatrix.h,v 1.11.2.19 2019/05/28 11:13:02 garloff Exp $
//----------------------------------------------------------

#ifndef TBCI_CRMATRIX_H
#define TBCI_CRMATRIX_H  

#include "tbci/matrix_sig.h"
#include "tbci/matrix.h"
#include "tbci/vector.h"

//#include <cstdlib> // for use of null macro

// Avoid -fguiding-decls
#if !defined(NO_GD) && !defined(AUTO_DECL)
# include "tbci/crmatrix_gd.h"
#endif

NAMESPACE_TBCI

// exception class
#ifndef TBCI_DISABLE_EXCEPT
//# include "except.h" is included in basics.h

class CRMatErr : public NumErr
{
    public:
        CRMatErr() 
        : NumErr("Error in crmatrix library") {}
        CRMatErr(const char* t, const long i = 0) 
        : NumErr(t, i) {}
        CRMatErr(const CRMatErr& ce)
        : NumErr(ce) {}
        virtual ~CRMatErr() throw() {}
};
#endif


#ifdef PRAGMA_I
# pragma interface "crmatrix.h"
#endif


//----------
// CRMatrix
//----------

template <typename T>
class CRMatrix : public Matrix_Sig<T>
{
protected:
        unsigned int n_rows, n_cols;
        unsigned int* rowSize;   // number of elements actually stored in each row
        unsigned int** colIndex; // column index of each stored element
        T** comp;
        mutable T dummy;
        // used by constructor and resize()
        void allocate(unsigned int rows, unsigned int columns);
        // used by destructor and "=" operator
        void destroy();
        // used by copy constr. and "="
        void copy(const CRMatrix<T>& m);

public:
        typedef T value_type;
        typedef T element_type;
        typedef T aligned_value_type TALIGN(MIN_ALIGN2);

        CRMatrix() { allocate(1,1); }
        CRMatrix(unsigned int rows, unsigned int columns) { allocate(rows,columns); }
        CRMatrix(const CRMatrix<T>& m) { copy(m); }
        CRMatrix(const Matrix/*_Sig*/<T>& mat);
        ~CRMatrix() { destroy(); }

        // allow instantiation (Matrix_Sig)
        /*virtual*/ static const char* mat_info ()  { return ("CRMatrix"); }

        CRMatrix<T>& clear ();
        CRMatrix<T>& fill (const T&);
        CRMatrix<T>& setunit (const T& = (T)1);

        inline const T& operator () (unsigned int row, unsigned int column) const;
        inline const T& get (unsigned int row, unsigned int column) const
        { return (*this) (row, column); }
        CRMatrix<T>& setval (const T& z, unsigned int row, unsigned int column);
        T& setval (unsigned int row, unsigned int column);
        T& operator () (unsigned int row, unsigned int column)
        { return setval (row, column); }
        
        // rowwise access
        TVector<T> get_row(unsigned int row) const;
        void set_row(const Vector<T>& v, unsigned int row);
        void set_row(const TVector<T>& tv, unsigned int row);
        void set_row(const TSVector<T>& tsv, unsigned int row);
        // replace a certain row using pre-compressed array storage
        void set_row( T* value, unsigned int row, unsigned int count, unsigned int* index);
        
        // query matrix dimensions
        unsigned int rows()    const {return n_rows;}
        unsigned int columns() const {return n_cols;}
        unsigned int size()    const; // number of non-zero elements

        // change matrix dimensions
        inline void resize(unsigned int newRows, unsigned int newColumns);

        // matrix-matrix assignment and comparison
        CRMatrix<T>& operator= (const CRMatrix<T>& m);
        bool operator == (const CRMatrix<T>& m) const;
        bool operator != (const CRMatrix<T>& m) const     { return !(this->operator == (m)); }

        // matrix negation
        CRMatrix<T> operator - () const;

        // matrix - vector multiplication
        TVector<T> operator * (const Vector<T>& v) const;
        TVector<T> operator * (const TVector<T>& tv) const;
        TVector<T> operator * (TSVector<T>& tsv) const;
        
        // For friend void MatVecMult FGD (Vector<T>& res, const CRMatrix<T>& m, const Vector<T>& v);
        void MatVecMult (Vector<T>& res, /*ourself ,*/ const Vector<T>& v) const;
        void MatVecMult (T* v, T* res); // res <- A*v , for ARPACK++

        // transpose-vector multiplication
        TVector<T> transMult(const Vector<T>& v) const;
        TVector<T> transMult(const TVector<T>& tv) const;
        TVector<T> transMult(const TSVector<T>& tsv) const;

        CRMatrix<T> operator* (const T& z) const;        // matrix-scalar multiplication
        CRMatrix<T>& operator*= (const T& z);          // fast in-place matrix-scalar multiplication
        CRMatrix<T> operator/ (const T& z) const;        // matrix-scalar division
        CRMatrix<T>& operator/= (const T& z);        // fast in-place matrix-scalar division

        // for friend CrMatrix,T> operator* (const T&, const CRMatrix<T>&);
        CRMatrix<T> mult (const T& z) const;

        CRMatrix<T>& swap (CRMatrix<T>&);
        CRMatrix<T>  transposed_copy() const;
        CRMatrix<T>& transpose();

        // Output operations
        friend STD__ ostream& operator<< FGD (STD__ ostream& stream, const CRMatrix<T>& m);
};


//-------------------------------
// CRMatrix function definitions 
//-------------------------------


template <typename T>
CRMatrix<T>::CRMatrix(const Matrix/*_Sig*/<T>& mat)
{
        allocate (mat.rows(), mat.columns());
        for (unsigned i = 0; i < n_rows; ++i)
                for (unsigned j = 0; j < n_cols; ++j)
                        (*this) (i, j) = mat( i, j);
}

// element access (read)
template <typename T>
inline const T& CRMatrix<T>::operator() (unsigned int row, unsigned int column) const
{
        BCHK(row >= n_rows || column >= n_cols, CRMatErr, Illegal index, (row<column ? row : column), dummy = (T)0);
        for (unsigned int j=0; j<rowSize[row]; j++)
                if (colIndex[row][j] == column)
                        return comp[row][j];
        return (dummy = (T)0);
}


// number of non-zero elements
template <typename T>
unsigned int CRMatrix<T>::size() const
{
        unsigned int erg=0;
        for (unsigned int i=0; i<n_rows; ++i)
                erg += rowSize[i];
        return erg;
}



// change matrix dimensions
template <typename T>
inline void CRMatrix<T>::resize(unsigned int newRows, unsigned int newColumns)
{
        if (newRows != n_rows || newColumns != n_cols)
                destroy();
        allocate(newRows, newColumns);
}


// element access (write)
template <typename T>
CRMatrix<T>& CRMatrix<T>::setval(const T& z, unsigned int row, unsigned int column)
{
        BCHK(row >= n_rows, CRMatErr, Illegal row index, row, *this);
        BCHK(column >= n_cols, CRMatErr, Illegal col index, column, *this);
        unsigned int j;
        for (j=0; j<rowSize[row]; ++j)
                if (colIndex[row][j] == column) {
                        comp[row][j] = z;
                        return *this;
                }
        
        // element not found -- dynamically growing required...
        if (z==0.0) 
                return *this; // no need for growing
        // create new column...
        unsigned int* newIndex = new unsigned int[rowSize[row]+1];
        BCHK( newIndex==NULL, CRMatErr, Out of Memory, rowSize[row]+1, *this);
        T* newComp = new T[rowSize[row]+1];
        BCHK( newComp==NULL, CRMatErr, Out of Memory, rowSize[row]+1, *this);
        // store "old" elements before new element
        for (j=0; j<rowSize[row]; ++j)
                if (colIndex[row][j] < column) {
                        newIndex[j] = colIndex[row][j];
                        newComp[j] = comp[row][j];
                } else
                        break;
        // store new element
        newIndex[j] = column;
        newComp[j] = z;
        // store "old" elementes behind new element
        for (unsigned int i=j; i<rowSize[row]; i++) {
                newIndex[i+1] = colIndex[row][i];
                newComp[i+1] = comp[row][i];
        }
        // kill old column
        if (rowSize[row] != 0) {
                delete[] colIndex[row];
                delete[] comp[row];
        }
        rowSize[row]++;
        colIndex[row] = newIndex;
        comp[row] = newComp;
        return *this;
}

template <typename T>
T& CRMatrix<T>::setval(unsigned int row, unsigned int column)
{
        BCHK(row >= n_rows || column >= n_cols, CRMatErr, "Illegal index", (row<column ? row : column), comp[0][0]);
        unsigned int j;
        for (j=0; j<rowSize[row]; j++)
                if (colIndex[row][j] == column)
                {
                        return comp[row][j];
                }
        
        // element not found -- dynamically growing required...
        // create new column...
        unsigned int* newIndex = new unsigned int[rowSize[row]+1];
        BCHK( newIndex==NULL, CRMatErr, Out of Memory, rowSize[row]+1, comp[0][0]);
        T* newComp = new T[rowSize[row]+1];
        BCHK( newComp==NULL, CRMatErr, Out of Memory, rowSize[row]+1, comp[0][0]);
        // store "old" elements before new element
        for (j=0; j<rowSize[row]; j++)
                if (colIndex[row][j] < column)
                {
                        newIndex[j] = colIndex[row][j];
                        newComp[j] = comp[row][j];
                }
                else break;
        // store new element
        unsigned int pos_store(j);
        newIndex[j] = column;
//      newComp[j] = z;
        // store "old" elementes behind new element
        for (unsigned int i=j; i<rowSize[row]; i++)
        {
                newIndex[i+1] = colIndex[row][i];
                newComp[i+1] = comp[row][i];
        }
        // kill old column
        if (rowSize[row] != 0)
        {
                delete[] colIndex[row];
                delete[] comp[row];
        }
        rowSize[row]++;
        colIndex[row] = newIndex;
        comp[row] = newComp;

        // Init with zero!!
        comp[row][pos_store]=0;
        return comp[row][pos_store];
}
        
// extract a certain row
template <typename T>
TVector<T> CRMatrix<T>::get_row(unsigned int row) const
{
        TVector<T> RowVec(n_cols);
        BCHK(row >= n_rows, CRMatErr, Illegal index, row, RowVec);
        for (unsigned int j=0; j<n_cols; j++)
                RowVec.set((*this)(row,j),j);
        return RowVec;
}

// replace a certain row
template <typename T>
void CRMatrix<T>::set_row( const Vector<T>& v, unsigned int row)
{
        BCHKNR(row >= n_rows, CRMatErr, Illegal index, row);
        unsigned int count = 0;
        unsigned int j;
        for (j=0; j<v.size(); j++)
                if (v(j) != 0.0) count++;
        if (rowSize[row] != count)
        {
                if (rowSize[row]!=0) // kill old row...
                {
                        delete[] colIndex[row];
                        delete[] comp[row];
                }
                if ((rowSize[row]=count) != 0) // create new row...
                {
                        colIndex[row] = new unsigned int[count];
                        BCHKNR(colIndex[row]==NULL, CRMatErr, Out of memory, count);
                        comp[row] = new T[row];
                        BCHKNR(comp[row]==NULL, CRMatErr, Out of memory, row);
                }
                else
                {
                        colIndex[row] = NULL;
                        comp[row] = NULL;
                }
        }
        count = 0;
        for (j=0; j<v.size(); j++) // copy values...
        if (v(j)!=0.0)
        {
                colIndex[row][count] = j;
                comp[row][count] = v(j);
                count++;
        }
}

template <typename T>
void CRMatrix<T>::set_row( const TVector<T>& tv, unsigned int row)
{
        Vector<T> v(tv);
        set_row(v,row);
}

template <typename T>
void CRMatrix<T>::set_row( const TSVector<T>& tsv, unsigned int row)
{
        Vector<T> v(tsv);
        set_row(v,row);
}

// replace a certain row
// using pre-compressed array storage
template <typename T>
void CRMatrix<T>::set_row( T* value,unsigned int row, unsigned int count, unsigned int* index)
{
        BCHKNR(row >= n_rows || count >= n_cols, CRMatErr, Illegal index, (row<count ? row : count));
        if (rowSize[row] != count)
        {
                if (rowSize[row]!=0) // kill old row...
                {
                        delete[] colIndex[row];
                        delete[] comp[row];
                }
                if ((rowSize[row]=count) != 0) // create new row...
                {
                        colIndex[row] = new unsigned int[count];
                        BCHKNR(colIndex[row]==NULL, CRMatErr, Out of memory, count);
                        comp[row] = new T[row];
                        BCHKNR(comp[row]==NULL, CRMatErr, Out of memory, row);
                }
                else
                {
                        colIndex[row] = NULL;
                        comp[row] = NULL;
                }
        }
        for (unsigned int j=0; j<count; j++) // copy values...
        {
                colIndex[row][j] = index[j];
                comp[row][j] = value[j];
        }
}

// matrix-matrix assignment
template <typename T>
CRMatrix<T>&    CRMatrix<T>::operator= (const CRMatrix<T>& m)
{
        if (this == &m) return *this;
        if (n_rows==m.n_rows) // perform rowwise copy
        {
                for (unsigned int i=0; i<n_rows; i++)
                {
                        if (rowSize[i]!=m.rowSize[i]) // different row sizes
                        {
                                if (rowSize[i] != 0) // kill old row...
                                {
                                        delete[] colIndex[i];
                                        delete[] comp[i];
                                }
                                if ((rowSize[i]=m.rowSize[i]) != 0) // create new row...
                                {
                                        colIndex[i] = new unsigned int[rowSize[i]];
                                BCHK(colIndex[i]==NULL, CRMatErr, Out of memory, rowSize[i], *this);
                                        comp[i] = new T[rowSize[i]];
                                BCHK(comp[i]==NULL, CRMatErr, Out of memory, rowSize[i], *this);
                                }
                                else
                                {
                                        colIndex[i] = NULL;
                                        comp[i] = NULL;
                                }
                        }
                        for (unsigned int j=0; j<rowSize[i]; j++) // copy values...
                        {
                                colIndex[i][j] = m.colIndex[i][j];
                                comp[i][j] = m.comp[i][j];
                        }
                } // for each row...
                return *this;
        } // end of rowwise copy...
        destroy();
        copy(m);
        return *this;
}


// KG 980721
template <typename T>
bool CRMatrix<T>::operator == (const CRMatrix<T>& cm) const
{
        if (n_rows != cm.n_rows || n_cols != cm.n_cols) 
                return false;
        for (unsigned int i=0; i<n_rows; ++i)
                for (unsigned int j=0; j<n_cols; ++j)
                        if ((*this)(i,j) != cm(i,j)) 
                                return false;
        return true;
}


// matrix negation
template <typename T>
CRMatrix<T> CRMatrix<T>::operator - () const
{
        CRMatrix<T> res(n_rows, n_cols);
        for (unsigned int i=0; i<n_rows; ++i)
                for (unsigned int j=0; j<n_cols; ++j)
                        if ((*this)(i,j) != (T)0) 
                                res.setval(-(*this)(i,j),i,j);
        return res; // JENS: avoid copy constructor soon
}


// matrix-vector multiplication
template <typename T>
TVector<T> CRMatrix<T>::operator * (const Vector<T>& v) const
{
        TVector<T> res(n_rows);
        BCHK(n_cols != v.size(), CRMatErr, Dimension conflict, v.size(), res);
        for (unsigned int i=0; i<n_rows; ++i) {
                T var = 0.0;
                for (unsigned int j=0; j<rowSize[i]; ++j)
                        var += comp[i][j] * v(colIndex[i][j]);
                res.set(var,i);
        }
        return res;
}


template <typename T>
inline TVector<T> CRMatrix<T>::operator * (const TVector<T>& tv) const
{
        Vector<T> v(tv);
        return (this->operator * (v));
}


template <typename T>
inline TVector<T> CRMatrix<T>::operator * (TSVector<T>& tsv) const
{
        TVector<T> res(n_rows);
        BCHK(n_cols != tsv.size(), CRMatErr, Dimension conflict, tsv.size(), res);
        for (unsigned int i=0; i<n_rows; ++i) {
                T var = 0.0;
                for (unsigned int j=0; j<rowSize[i]; ++j)
                        var += comp[i][j] * tsv.get(colIndex[i][j]);
                res.set(var,i);
        }
        tsv.destroy ();
        return res;
}


// matrix-vector multiplication
template <typename T>
void CRMatrix<T>::MatVecMult(Vector<T>& res, const Vector<T>& v) const
{
        BCHKNR(n_cols != v.size(), CRMatErr, Dimension conflict, v.size());
        for (unsigned int i=0; i<n_rows; ++i) {
                T var = 0.0;
                for (unsigned int j=0; j<rowSize[i]; ++j)
                        var += comp[i][j] * v(colIndex[i][j]);
                res(i) = var;
        }
}

INST(template <typename T> class CRMatrix friend void MatVecMult(Vector<T>&, const CRMatrix<T>&, const Vector<T>&);)
template <typename T>
inline void MatVecMult(Vector<T>& res, const CRMatrix<T>& m, const Vector<T>& v)
{ m.MatVecMult (res, v); }

template <typename T>
void CRMatrix<T>::MatVecMult (T* v, T* res) // res <- A*v , for ARPACK++
{
        // No Boundcheck possible!!
        for (unsigned int i=0; i<this->n_rows; ++i) {
                T var = 0.0;
                for (unsigned int j=0; j<this->rowSize[i]; ++j)
                        var += this->comp[i][j] * v[this->colIndex[i][j]];
                res[i] = var;
        }
}


// transpose-vector multiplication
template <typename T>
TVector<T> CRMatrix<T>::transMult(const Vector<T>& v) const
{
        // Optimzed version: 7.8.98 a. Ahland
        TVector<T> res(0, n_cols); // Zero initialisation is very important!!
        BCHK(n_rows != v.size(), CRMatErr, Dimension conflict, v.size(), res);

        unsigned int pos;
        for (unsigned int i=0; i<n_rows; ++i) {
                for (unsigned int j=0; j<rowSize[i]; ++j) {
                        pos = colIndex[i][j];
                        res.setval(pos)  += comp[i][j] * v(i);
                }
        }
        return res;
}


template <typename T>
TVector<T> CRMatrix<T>::transMult(const TVector<T>& tv) const
{
        Vector<T> v(tv);
        return(this->transMult(v));
}


template <typename T>
TVector<T> CRMatrix<T>::transMult(const TSVector<T>& tsv) const
{
        Vector<T> v(tsv);
        return(this->transMult(v));
}


// matrix-scalar multiplication
template <typename T>
CRMatrix<T> CRMatrix<T>:: operator* (const T& z) const
{
        CRMatrix<T> res(*this);
        for (unsigned int i=0; i<this->n_rows; ++i)
                for (unsigned int j=0; j<this->rowSize[i]; ++j)
                        res.comp[i][j] = this->comp[i][j] * z;
        return res; // JENS: avoid copy constructor soon!
}


// matrix-scalar multiplication
template <typename T>
CRMatrix<T> CRMatrix<T>::mult (const T& z) const
{
        CRMatrix<T> res(*this);
        for (unsigned int i=0; i<n_rows; ++i)
                for (unsigned int j=0; j<rowSize[i]; ++j)
                        res.comp[i][j] = z * comp[i][j];
        return res;
}

INST(template <typename T> class CRMatrix friend CRMatrix<T> operator * (const T&, const CRMatrix<T>&);)
template <typename T>
inline CRMatrix<T> operator* (const T& z, const CRMatrix<T>& m)
{ return m.mult (z); }


// fast in-place matrix-scalar multiplication
template <typename T>
CRMatrix<T>& CRMatrix<T>::operator*= (const T& z)
{
        for (unsigned int i=0; i<n_rows; ++i)
                for (unsigned int j=0; j<rowSize[i]; ++j)
                        comp[i][j] *= z;
        return *this;
}

// matrix-scalar division
template <typename T>
CRMatrix<T> CRMatrix<T>::operator/ (const T& z) const
{
        CRMatrix<T> res(*this);
        for (unsigned int i=0; i<n_rows; ++i)
                for (unsigned int j=0; j<this->rowSize[i]; ++j)
                        res.comp[i][j] = this->comp[i][j] / z;
        return res; // JENS: avoid copy constructor soon!
}


// fast in-place matrix-scalar division
template <typename T>
CRMatrix<T>& CRMatrix<T>::operator/= (const T& z)
{
        for (unsigned int i=0; i<n_rows; i++)
                for (unsigned int j=0; j<rowSize[i]; j++)
                        comp[i][j] /= z;
        return *this;
}

// allocate memory space
// used by constructor and resize()
template <typename T>
void CRMatrix<T>::allocate(unsigned int rows, unsigned int columns)
{
        BCHKNR(rows <= 0 || columns <= 0, CRMatErr, Zero space allocating, 0);
        n_rows = rows; dummy = 0;
        n_cols = columns;
        rowSize = new unsigned int[n_rows];
        BCHKNR(rowSize==NULL, CRMatErr, Out of memory, n_rows);
        colIndex = new unsigned int*[n_rows];
        BCHKNR(colIndex==NULL, CRMatErr, Out of memory, n_rows);
        comp = new T*[n_rows];
        BCHKNR(comp==NULL, CRMatErr, Out of memory, n_rows);
        for (unsigned int i=0; i<n_rows; i++)
        {
                rowSize[i] = 0;
                colIndex[i] = NULL;
                comp[i] = NULL;
        }
}

// delete memory space
// used by "=" operator and destructor
template <typename T>
void CRMatrix<T>::destroy()
{
        for (unsigned int i=0; i<n_rows; ++i) {
                if (colIndex[i]!=NULL) delete[] colIndex[i];
                if (comp[i]!=NULL) delete[] comp[i];
        }
        if (rowSize != NULL)
                delete[] rowSize;
        if (comp != NULL) 
                delete[] comp;
        if (colIndex != NULL) 
                delete[] colIndex;
}

// copy to new memory space
// used by "=" operator and copy constructor
template <typename T>
void CRMatrix<T>::copy(const CRMatrix<T>& m)
{
        n_rows = m.n_rows;
        n_cols = m.n_cols; dummy = 0;
        rowSize = new unsigned int[n_rows];
        BCHKNR(rowSize==NULL, CRMatErr, Out of memory, n_rows);
        colIndex = new unsigned int*[n_rows];
        BCHKNR(colIndex==NULL, CRMatErr, Out of memory, n_rows);
        comp = new T*[n_rows];
        BCHKNR(comp==NULL, CRMatErr, Out of memory, n_rows);
        for (unsigned int i=0; i<n_rows; ++i) {
                if ((rowSize[i]=m.rowSize[i]) != 0) {
                        colIndex[i] = new unsigned int[rowSize[i]];
                        BCHKNR(colIndex[i]==NULL, CRMatErr, Out of memory, rowSize[i]);
                        comp[i] = new T[rowSize[i]];
                        BCHKNR(comp[i]==NULL, CRMatErr, Out of memory, rowSize[i]);
                        for (unsigned int j=0; j<rowSize[i]; ++j) {
                                colIndex[i][j] = m.colIndex[i][j];
                                comp[i][j] = m.comp[i][j];
                        }
                } else {
                        colIndex[i] = NULL;
                        comp[i] = NULL;
                }
        }
}



template <typename T>
inline CRMatrix<T>& CRMatrix<T>::fill (const T& val)
{
        for (unsigned int i=0; i<n_rows; ++i)
                for (unsigned int j=0; j<rowSize[i]; ++j)
                        comp[i][j] = val;
        return *this;
}

template <typename T>
inline CRMatrix<T>& CRMatrix<T>::clear ()
{ return fill ((T)0); }

template <typename T>
inline CRMatrix<T>& CRMatrix<T>::setunit (const T& val)
{
        BCHK (n_rows != n_cols, CRMatErr, setunit is meaningless for non-square matrices, n_cols, *this);
        clear ();
        for (unsigned int i = 0; i < n_cols; ++i)
                setval (val, i, i);
        return *this;
}

template <typename T>
CRMatrix<T>& CRMatrix<T>::swap(CRMatrix<T>& m)
{       
        SWAP(n_rows, m.n_rows); SWAP(n_cols, m.n_cols);
        SWAP(rowSize, m.rowSize); SWAP(colIndex, m.colIndex);
        SWAP(comp, m.comp);
        return *this;
}

template <typename T>
CRMatrix<T> CRMatrix<T>::transposed_copy() const
{
        CRMatrix<T> m(columns(), rows());
        for (unsigned r = 0; r < rows(); ++r)
                for (unsigned c = 0; c < columns(); ++c)
                        m(c,r) = (*this)(r,c);
        return m;
}

template <typename T>
inline CRMatrix<T>& CRMatrix<T>::transpose()
{
        CRMatrix<T> tr(this->transposed_copy());
        return swap(tr);
}

INST(template <typename T> class CRMatrix friend CRMatrix<T> transpose(const CRMatrix<T>& crm);)
template <typename T>
inline CRMatrix<T> transpose(const CRMatrix<T>& crm)
{
        return crm.transposed_copy();
}


template <typename T>
STD__ ostream& operator<< (STD__ ostream& stream, const CRMatrix<T>& m)
{
        for (unsigned int i = 0; i < m.rows(); ++i) {
                stream << "|\t";
                for (unsigned int j = 0; j < m.columns(); ++j)
                        stream << m(i,j) << "\t";
                stream << "|\n";
        }
        return stream << STD__ flush;
}

NAMESPACE_END

#endif  /* TBCI_CRMATRIX_H */