cscmatrix.h

Go to the documentation of this file.

//----------------------------------------------------------
// version 1.0, 08.1998
//
// Andreas Ahland
//
// $Id: cscmatrix.h,v 1.18.2.40 2022/11/03 17:28:10 garloff Exp $
//----------------------------------------------------------

#ifndef TBCI_CSCMATRIX_H
#define TBCI_CSCMATRIX_H

#include "tbci/vector.h"
#include "tbci/matrix_sig.h"
#include "tbci/f_matrix.h"
#include "tbci/band_matrix.h"

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

NAMESPACE_TBCI

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

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

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

template <typename T> class CSCMatrix;
template <typename T> class F_TSMatrix;
template <typename T> class F_TMatrix;
template <typename T> class F_Matrix;

//----------
// CSCMatrix
//----------
template <typename T>
class CSCMatrix : public Matrix_Sig<T>
{
        friend class F_TSMatrix<T>;

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

        CSCMatrix() { allocate(1,1); }
        CSCMatrix(const T& val, const unsigned int rows, 
                          const unsigned int columns, const unsigned int nnzeros=1) 
        {   allocate(rows, columns, nnzeros); 
                if (val != (T)0) fill(val); }
        CSCMatrix(const unsigned int rows, const unsigned int columns, 
                          const unsigned int nnzeros=1) 
        { allocate(rows, columns, nnzeros); }
        CSCMatrix(const CSCMatrix<T>& m) 
        { copy(m); }
        CSCMatrix(const F_Matrix<T>& m, const double tol=0)
        { allocate(m.rows(), m.columns());      fill(m,tol); }
        CSCMatrix(const BdMatrix<T>& m, const double tol=0)
        { allocate(m.rows(), m.columns());      fill(m,tol); }
        CSCMatrix(const   Matrix<T>& m, const double tol=0)
        { allocate(m.rows(), m.columns());      fill(m,tol); }
        ~CSCMatrix() 
        { destroy (); }
        
        operator F_TMatrix<T> () const;

        /*virtual*/ static const char* mat_info ()  { return ("CSCMatrix"); }


        unsigned int rows()    const  { return n_rows; }
        unsigned int columns() const  { return n_cols; }
        unsigned int size()    const  { return n_size; } // number of non-zero elements


        CSCMatrix<T>& resize(const unsigned int newRows, 
                             const unsigned int newColumns, 
                             const unsigned int nnzeros = 1);

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

        typename tbci_traits<T>::const_refval_type
          operator() (unsigned int row, unsigned int column) const HOT;
        typename tbci_traits<T>::const_refval_type
          get (unsigned int row, unsigned int column) const HOT;

        CSCMatrix<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); }
        T* dataPointer()                { return comp; }
        unsigned int* columnPointer()   { return pcol; }
        unsigned int* rowIndexPointer() { return irow; }

        bool operator == (const CSCMatrix<T>& m) const;
        bool operator != (const CSCMatrix<T>& m) const    
        { return !(this->operator == (m)); }
        CSCMatrix<T>& operator= (const CSCMatrix<T>& m);

        //CSCMatrix<T>& do_import (const Matrix_Sig<T>& M)
        template <typename MatType>
        CSCMatrix<T>& do_import  /*FGD*/ (const MatType& M)
        {
                destroy (); allocate(M.rows(), M.columns(), 20);
                fill (M);
                return *this;
        }
        
        //void do_export (Matrix_Sig<T>& M)
        template <typename MatType>
        void do_export /*FGD*/ (MatType& M)
        {
                M.clear (); //Matrix should be cleared!
                for (unsigned int i=0; i<n_cols; ++i)
                        //for (unsigned int j=pcol[i]; j<pcol[i+1]; ++j)  M.setval(comp[j], i,irow[j]);
                        for (unsigned int j=pcol[i]; j<pcol[i+1]; ++j)  
                                M(i,irow[j]) = comp[j];
        }

        CSCMatrix<T> operator - () const;

        CSCMatrix<T> operator + (const CSCMatrix<T>&) const;
        CSCMatrix<T> operator - (const CSCMatrix<T>&) const;

        F_TMatrix<T> operator * (const CSCMatrix<T>&) const;

        //friend  F_TMatrix<T> operator* FGD (const F_Matrix<T>& m1, const CSCMatrix<T>& m2);
        F_TMatrix<T> multf (const F_Matrix<T>&) const;
        F_TMatrix<T> mult  (const F_Matrix<T>&) const;
        F_TMatrix<T> mult1 (const F_Matrix<T>&) const;

        TVector<T> operator * (const Vector<T>& v) const HOT;
        TVector<T> operator * (const TVector<T>& tv) const HOT;
        TVector<T> operator * (TSVector<T>& tsv) const HOT;
        TVector<T> cscm_vec_mul_exact  (const Vector<T>& V) const;
        void MatVecMult (Vector<T>& res, const Vector<T>& v) const HOT;
        void MatVecMult (T* v, T* res) HOT; // res <- A*v , for ARPACK++

        TVector<T> transMult(const Vector<T>& v) const HOT;
        TVector<T> transMult(const TVector<T>& tv) const HOT;
        TVector<T> transMult(const TSVector<T>& tsv) const HOT;


        CSCMatrix<T> operator* (const T& z) const;       // matrix-scalar multiplication
        CSCMatrix<T> mult (const T&) const;

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

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

        friend STD__ ostream& operator<< FGDT (STD__ ostream& stream, const CSCMatrix<T>& m);

protected:
        unsigned int n_rows, n_cols, n_size, n_max_size;
        unsigned int* pcol;  // pointer to the beginning of each column
        unsigned int* irow;  // row index of all nonzero elements
        T* comp;             // storage vector
        mutable T dummy;

        void allocate (unsigned int rows, unsigned int columns, unsigned int nnzeros=1);
        void destroy  ();       // used by destructor and "=" operator
        void copy     (const CSCMatrix<T>& m);  // used by copy constr. and "="
        void insert   (const unsigned int column, const unsigned int pos); //dynamic inserting

        //void fill (const Matrix_Sig<T>& M, double tol=0)
        template <typename MatType>
        void fill /*FGD*/ (const MatType& M, const double tol = 0.0)
        {
                unsigned int count = 0, i, j;
                //count nonzeros
                for (j = 0; j < n_cols; ++j)
                        for (i = 0; i < n_rows; ++i)    
                                if (MATH__ fabs(M(i,j)) > tol) 
                                        ++count; 
                STD__ cout << "Matrix conversion, NNZ=" << count << STD__ endl;
                resize(M.rows(), M.columns(), count);
                for (j = 0; j < n_cols; ++j)
                        for (i = 0; i < n_rows; ++i)  
                                if (MATH__ fabs(M(i,j)) > tol) 
                                        setval(M(i,j), i,j); 
        }
};


template <typename T>
STD__ ostream& operator<< (STD__ ostream& stream, const CSCMatrix<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*/;
}

//-------------------------------
// CSCMatrix function definitions 
//-------------------------------

template <typename T>
F_TMatrix<T> CSCMatrix<T>::multf (const F_Matrix<T>& m1) const
{       
        F_TMatrix<T> res((T)0, m1.rows(), columns());
        BCHK(rows() != m1.columns(), CSCMatErr, Sizes dont match in multf, m1.columns(), res);
        if (do_exactsum2()) {
                for (unsigned int i = 0; i < m1.rows(); ++i)
                        for (unsigned int j = 0; j < columns(); ++j) {
                                REGISTER T val = 0;
                                REGISTER T corr = 0;
                                for (unsigned int x = pcol[j]; x < pcol[j+1]; ++x) {
                                        const REGISTER T y = m1(i, irow[x] )*comp[x];
                                        const REGISTER T t = val + y;
                                        corr += (t-val) - y;
                                        val = t;
                                }
                                res.setval(val-corr, i, j);
                        }
        } else {
                for (unsigned int i = 0; i < m1.rows(); ++i)
                        for (unsigned int j = 0; j < columns(); ++j) {
                                REGISTER T val = 0;
                                for (unsigned int x = pcol[j]; x < pcol[j+1]; ++x)
                                        val += m1(i, irow[x] )*comp[x];
                                res.setval(val, i, j);
                        }
        }
        return res;
}

INST(template<typename T> class CSCMatrix<T> friend F_TMatrix<T> operator* (const F_Matrix<T>&, const CSCMatrix<T>&);)
template <typename T>
inline F_TMatrix<T> operator * (const F_Matrix<T>& m1, const CSCMatrix<T>& m2)
{ return m2.multf (m1); }

template <typename T>
F_TMatrix<T> CSCMatrix<T>::mult (const F_Matrix<T>& m1) const
{       
        F_TMatrix<T> res((T)0, rows(), m1.columns());
        BCHK(columns() != m1.rows(), CSCMatErr, Sizes dont match in mult, columns(), res);
        if (do_exactsum2()) {
                for (unsigned int i = 0; i < rows(); ++i)
                        for (unsigned int j = 0; j < m1.columns(); ++j) {
                                REGISTER T val = 0;
                                REGISTER T corr = 0;
                                for (unsigned int x = 0; x < columns(); ++x) {
                                        const REGISTER T y = (*this)(i, x) * m1(x, j);
                                        const REGISTER T t = val + y;
                                        corr += (t-val) - y;
                                        val = t;
                                }
                                res.setval(val-corr, i, j);
                        }
        } else {
                for (unsigned int i = 0; i < rows(); ++i)
                        for (unsigned int j = 0; j < m1.columns(); ++j) {
                                REGISTER T val = 0;
                                for (unsigned int x = 0; x < columns(); ++x) 
                                        val += (*this)(i, x) * m1(x, j);
                                res.setval(val, i, j);
                        }
        }
        return res;
}

template <typename T>
F_TMatrix<T> CSCMatrix<T>::mult1 (const F_Matrix<T>& m1) const
{       
        F_TMatrix<T> res((T)0, rows(), m1.columns());
        BCHK(columns() != m1.rows(), CSCMatErr, Sizes dont match in mult, columns(), res);
        for (unsigned int c = 0; c < columns(); ++c) 
                for (unsigned int off = pcol[c]; off < pcol[c+1]; ++off) {
                        /* No exact summing possible */
                        const unsigned int r = irow[off];
                        const T el = comp[off]; // (*this)(r,c)
                        for (unsigned int j = 0; j < m1.columns(); ++j)
                                res.setval(r, j) += el * m1(c, j);
                }
        return res;
}

INST(template<typename T> class CSCMatrix<T> friend F_TMatrix<T> operator* (const CSCMatrix<T>&, const F_Matrix<T>&);)
template <typename T>
inline F_TMatrix<T> operator * (const CSCMatrix<T>& m1, const F_Matrix<T>& m2)
{ return m1.mult1 (m2); }

template <typename T>
F_TMatrix<T> CSCMatrix<T>::operator * (const CSCMatrix<T>& m1) const
{
        F_TMatrix<T> res((T)0, rows(), m1.columns());
        BCHK(columns() != m1.rows(), CSCMatErr, Sizes dont match in mult, columns(), res);
        if (do_exactsum2()) {
                for (unsigned int r0 = 0; r0 < rows(); ++r0)
                        for (unsigned int c1 = 0; c1 < m1.columns(); ++c1) {
                                REGISTER T val = 0;
                                REGISTER T corr = 0;
                                for (unsigned int off1 = m1.pcol[c1]; off1 < m1.pcol[c1+1]; ++off1) {
                                        const REGISTER T y = (*this)(r0, m1.irow[off1]) * m1.comp[off1];
                                        const REGISTER T t = val + y;
                                        corr += (t-val) - y;
                                        val = t;
                                }
                                res.setval(val-corr, r0, c1);
                        }
        } else {
                for (unsigned int r0 = 0; r0 < rows(); ++r0)
                        for (unsigned int c1 = 0; c1 < m1.columns(); ++c1) {
                                REGISTER T val = 0;
                                for (unsigned int off1 = m1.pcol[c1]; off1 < m1.pcol[c1+1]; ++off1)
                                        val += (*this)(r0, m1.irow[off1]) * m1.comp[off1];
                                res.setval(val, r0, c1);
                        }
        }
        return res;
}

template <typename T>
CSCMatrix<T> CSCMatrix<T>::operator + (const CSCMatrix<T>& m1) const
{
        CSCMatrix<T> res(rows(), columns(), size());
        BCHK(columns() != m1.columns(), CSCMatErr, Columns dont match in add, columns(), res);
        BCHK(rows() != m1.rows(), CSCMatErr, Rows dont match in add, rows(), res);
        for (unsigned int c = 0; c < columns(); ++c) {
                unsigned int off0 = pcol[c];
                unsigned int off1 = m1.pcol[c];
                while (off0 < pcol[c+1] || off1 < m1.pcol[c+1]) {
                        // only one is valid ...
                        if (off0 >= pcol[c+1]) {
                                res(m1.irow[off1], c) = m1.comp[off1]; 
                                ++off1;
                        } else if (off1 >= m1.pcol[c+1]) {
                                res(irow[off0], c) = comp[off0]; 
                                ++off0;
                        } else if (irow[off0] == m1.irow[off1]) {
                                // both are valid
                                res(irow[off0], c) = comp[off0] + m1.comp[off1];
                                ++off0; ++off1;
                        } else if (irow[off0] < m1.irow[off1]) {
                                res(irow[off0], c) = comp[off0];
                                ++off0;
                        } else {
                                res(m1.irow[off1], c) = m1.comp[off1];
                                ++off1;
                        }
                }
        }
        return res;
}

template <typename T>
CSCMatrix<T> CSCMatrix<T>::operator - (const CSCMatrix<T>& m1) const
{
        CSCMatrix<T> res(rows(), columns(), size());
        BCHK(columns() != m1.columns(), CSCMatErr, Columns dont match in sub, columns(), res);
        BCHK(rows() != m1.rows(), CSCMatErr, Rows dont match in sub, rows(), res);
        for (unsigned int c = 0; c < columns(); ++c) {
                unsigned int off0 = pcol[c];
                unsigned int off1 = m1.pcol[c];
                while (off0 < pcol[c+1] || off1 < m1.pcol[c+1]) {
                        // only one is valid ...
                        if (off0 >= pcol[c+1]) {
                                res(m1.irow[off1], c) = -m1.comp[off1]; 
                                ++off1; continue;
                        }
                        if (off1 >= m1.pcol[c+1]) {
                                res(irow[off0], c) = comp[off0]; 
                                ++off0; continue;
                        }
                        // both are valid
                        if (irow[off0] == m1.irow[off1]) {
                                res(irow[off0], c) = comp[off0] - m1.comp[off1];
                                ++off0; ++off1; //continue;
                        } else if (irow[off0] < m1.irow[off1]) {
                                res(irow[off0], c) = comp[off0];
                                ++off0; //continue;
                        } else {
                                res(m1.irow[off1], c) = -m1.comp[off1];
                                ++off1; //continue;
                        }
                }
        }
        return res;
}

template <typename T>
inline typename tbci_traits<T>::const_refval_type
  CSCMatrix<T>::get(unsigned int row, unsigned int col) const
{
        static T dummy = (T)0;
/*
        // Version 1 
        for (unsigned int j=pcol[col]; j<pcol[col+1]; j++)
        {
                if (irow[j] == row) return comp[j];
        }

        // Version 2
        unsigned int pos = pcol[col]; 
        // linear search ...
        while (pos < pcol[col+1] && irow[pos] < row) 
                ++pos;
        if (pos < n_size && irow[pos] == row) 
                return comp[pos];
        // Version 3

        unsigned int pos = pcol[col];
        // linear search ...
        while (pos < pcol[col+1]) {
                if (irow[pos] < row) 
                        ++pos;
                else if (irow[pos] > row)
                        break;
                else
                        return comp[pos];
        }
        // Element not found
        return (dummy = 0);
*/

/*
        // Version 4
        const unsigned int epos = pcol[col+1];
        for (unsigned int pos = pcol[col]; pos < epos; ++pos) {
                if (irow[pos] == row)
                        return comp[pos];
                else if (irow[pos] > row)
                        break;
        }
        return (dummy = 0);
*/
        // Version 5
        const unsigned int spos = pcol[col], epos = pcol[col+1];
        if (epos == spos)
                return (dummy = 0);
        if (irow[spos] > row || irow[epos-1] < row)
                return (dummy = 0);
        if (irow[spos] == row)
                return comp[spos];
        if (irow[epos-1] == row)
                return comp[epos-1];
        /* Binary search */
        if ((epos - spos) >= 32) {
                const unsigned long int off = _bin_search(irow, row, spos, epos);
                if (off != (unsigned long)(-1))
                        return comp[off];
                else
                        return (dummy = 0);
        }
        /* Linear search, be clever at which side to start */
        if (row <= (irow[spos] + irow[epos-1] + 1) / 2) {
                for (REGISTER unsigned int pos = spos; pos < epos; ++pos) {
                        if (irow[pos] == row)
                                return comp[pos];
                        else if (irow[pos] > row)
                                break;
                }
        } else {
                for (REGISTER unsigned int pos = epos; pos > spos; --pos) {
                        const unsigned po1 = pos-1;
                        if (irow[po1] == row)
                                return comp[po1];
                        else if (irow[po1] < row)
                                break;
                }
        }
        return (dummy = 0);
        // TODO: Should we do binary search?
}

template <typename T>
inline typename tbci_traits<T>::const_refval_type
  CSCMatrix<T>::operator() (const unsigned row, const unsigned col) const
{
        EXPCHK(row >= n_rows, CSCMatErr, Illegal row index, row, dummy);
        EXPCHK(col >= n_cols, CSCMatErr, Illegal col index, col, dummy);
        return get(row, col);
}

template <typename T>
CSCMatrix<T>::operator F_TMatrix<T> () const
{
        F_TMatrix<T> res((T)0, rows(), columns());
        for (unsigned int c = 0; c < columns(); ++c) 
                for (unsigned int off = pcol[c]; off < pcol[c+1]; ++off)
                        res.setval(comp[off], irow[off], c);
        return res;
}


template <typename T>
CSCMatrix<T>& CSCMatrix<T>::resize(const unsigned int newRows, 
                                   const unsigned int newColumns,
                                   const unsigned int nnzeros)
{
        if (UNLIKELY(newRows != n_rows && newColumns != n_cols)) {
                destroy ();
                allocate(newRows, newColumns, nnzeros);
        } else if (UNLIKELY(nnzeros != n_max_size)) {
                // TODO: Handle size 0
                unsigned int *newirow = new unsigned int[nnzeros];
                T            *newcomp = new T           [nnzeros];
                unsigned int i;
                const unsigned int com_sz = MIN(nnzeros, n_max_size);
                for (i = 0; i < com_sz; ++i) {
                        newirow[i] = irow[i];
                        newcomp[i] = comp[i];
                }
                for (; i < nnzeros; ++i)
                        newirow[i] = 0;
                n_max_size = nnzeros;
                delete[] irow;
                delete[] comp;
                irow = newirow;
                comp = newcomp;
        } /* else do nothing */
        return *this;
}


template <typename T>
void CSCMatrix<T>::insert(const unsigned int column, const unsigned int pos)
{
        if (UNLIKELY(n_size >= n_max_size)) {
                //not enough space left
                n_max_size += 16;
                unsigned int* new_irow = new unsigned int [n_max_size];
                BCHKNR(new_irow==NULL, CSCMatErr, Out of memory, n_max_size);
                T*            new_comp = new T            [n_max_size];
                BCHKNR(new_comp==NULL, CSCMatErr, Out of memory, n_max_size);

                // if data exist, copy data
                if (n_size > 0) {
                        unsigned int i;
                        for (i = 0; i < pos; ++i)       { 
                                // copy data up before
                                new_irow[i] = irow[i];
                                new_comp[i] = comp[i];
                        }
                        for (i = pos; i < n_size; ++i) {
                                // copy data behind insertion
                                new_irow[i+1] = irow[i];
                                new_comp[i+1] = comp[i];
                        }
                        // Following elements are not yet used
                        for (i = n_size+1; i < n_max_size; ++i) 
                                new_irow[i] = 0;
                }

                if (irow != NULL) 
                        delete[] irow;
                if (comp != NULL) 
                        delete[] comp;
                irow = new_irow;
                comp = new_comp;
        } else if (n_size > 0) { 
                // if data exist, copy data
                for (int j = n_size-1; j >= (int)pos; --j) {
                         // copy data behind insertion
                        irow[j+1] = irow[j];
                        comp[j+1] = comp[j];
                }
        }
        for (unsigned int i = column+1; i < n_cols+1; ++i)
                ++pcol[i];
        ++n_size;
}

// element access (write)
template <typename T>
CSCMatrix<T>& CSCMatrix<T>::setval(const T& z, unsigned int row, unsigned int column)
{
        setval(row, column) = z;
        return *this;
}

template <typename T>
T& CSCMatrix<T>::setval(unsigned int row, unsigned int column)
{
        EXPCHKNR(row >= n_rows, CSCMatErr, Illegal row index, row);
        EXPCHKNR(column >= n_cols, CSCMatErr, Illegal col index, column);
        unsigned int pos  = pcol[column]; 
        const unsigned int posn = pcol[column+1]; 
        while (pos < posn && irow[pos] < row) 
                ++pos;
        if (pos < posn && irow[pos] == row) 
                return comp[pos];

        // element not found -- dynamically growing required...
        insert(column, pos);

        // Write new element
        irow[pos] = row;
        comp[pos] = T(0);   // Init with zero!!

        return comp[pos];
}


// matrix-matrix assignment
template <typename T>
CSCMatrix<T>&   CSCMatrix<T>::operator= (const CSCMatrix<T>& m)
{
        if (this == &m) 
                return *this;
        destroy(); copy(m);
        return *this;
}


// KG 19980721
template <typename T>
bool CSCMatrix<T>::operator == (const CSCMatrix<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 (operator()(i,j) != cm(i,j)) 
                                 return false;
                 }
         }
         return true;
}


// matrix negation
template <typename T>
CSCMatrix<T> CSCMatrix<T>::operator - () const
{
        CSCMatrix<T> res(n_rows, n_cols);
        for (unsigned int i = 0; i < n_size; ++i) 
                res.comp[i] = -comp[i];
        return res; // avoid copy constructor soon?
}


#if 0
template <typename T>
TVector<T> CSCMatrix<T>::cscm_vec_mul_exact (const Vector<T>& v) const
{
        TVector<T> res(n_rows);
        BCHK(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
        for (unsigned row = 0; row <= n_rows; ++row) {
                REGISTER T val(0.0), corr(0.0);
                /* Search for all elements in row row */
                for (unsigned int c = 0; c < n_cols; ++c) {
                        const REGISTER T el = this->get(row, c);
                        if (el) {
                                const REGISTER T y = el*v(c);
                                const REGISTER T t = val+y;
                                corr += (t-val) - y;
                                val = t;
                        }
                }
                res.setval(row) = val - corr;
        }
        return res;
}
#else
template <typename T>
TVector<T> CSCMatrix<T>::cscm_vec_mul_exact(const Vector<T>& v) const
{
        Vector<T> res((T)0, n_rows), corr((T)0, n_rows);
        BCHK(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
        for (unsigned int i = 0; i < n_cols; ++i) {
                for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) {
                        const unsigned row = irow[j];
                        const REGISTER T y = comp[j]*v(i);
                        const REGISTER T t = res.get(row) + y;
                        corr.setval(row) += (t-res.get(row)) - y;
                        res.setval(row) = t;
                }
        }
        return res-corr;
}
#endif

// matrix-vector multiplication
template <typename T>
TVector<T> CSCMatrix<T>::operator * (const Vector<T>& v) const
{
        if (do_exactsum2())
                return this->cscm_vec_mul_exact(v);
        TVector<T> res(0, n_rows);
        BCHK(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
        /* Exact summing not possible */
        for (unsigned int i = 0; i < n_cols; ++i) {
                for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
                        res.setval(irow[j]) += comp[j]*v(i);  
        }
        return res;
}

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

/* TODO: Implement exact summation for CSCM*TSV and MatVecMult */
template <typename T>
inline TVector<T> CSCMatrix<T>::operator * (TSVector<T>& tsv) const
{
        TVector<T> res(0, n_rows);
        BCHK(n_cols != tsv.size(), CSCMatErr, Dimension conflict, tsv.size(), res);
        /* Exact summing not possible */
        for (unsigned int i = 0; i < n_cols; ++i) {
                for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
                        res.setval(irow[j]) += comp[j] * tsv.get(i);
        }
        tsv.destroy ();
        return res;
}


// matrix-vector multiplication
template <typename T>
void CSCMatrix<T>::MatVecMult(Vector<T>& res, const Vector<T>& v) const
{
        BCHKNR(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size());
        res = T(0);
        /* Exact summing not possible */
        for (unsigned int i = 0; i < n_cols; ++i) {
                for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
                        res.setval(irow[j]) += comp[j]*v(i);  
        }
}

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


template <typename T>
void CSCMatrix<T>::MatVecMult (T* v, T* res) // res <- A*v , for ARPACK++
{
        STD__ cout << "CSCMatrix<T>::MatVecMult (T*, T*): Not tested!!" << STD__ endl;
        // No Boundcheck possible!!
        unsigned int i;
        for (i = 0; i < n_cols; ++i) 
                res[i] = 0;
        /* Exact summing not possible */
        for (i = 0; i < n_cols; ++i) {
                for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
                        res[irow[j]] += comp[j]*v[i];  
        }
}


// transpose-vector multiplication
template <typename T>
TVector<T> CSCMatrix<T>::transMult(const Vector<T>& v) const
{
        TVector<T> res(n_cols);
        BCHK(n_rows != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
        if (do_exactsum2()) {
                for (unsigned int i = 0; i < n_cols; ++i) {
                        REGISTER T val = 0, corr = 0;
                        for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) {
                                const REGISTER T y = comp[j]*v(irow[j]);
                                const REGISTER T t = val + y;
                                corr += (t-val) - y;
                                val = t;
                        }
                        res.set(val-corr, i);
                }
        } else {
                for (unsigned int i = 0; i < n_cols; ++i) {
                        T val = 0;
                        for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
                                val += comp[j]*v(irow[j]);
                        res.set(val, i);
                }
        }
        return res;
}


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


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


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


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


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

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

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


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


// allocate memory space
// used by constructor and resize()
template <typename T>
void CSCMatrix<T>::allocate(unsigned int rows, unsigned int columns, 
                            unsigned int nnzeros)
{
        BCHKNR(rows <= 0 || columns <= 0, CSCMatErr, Zero space allocating, 0);
        n_rows = rows;
        n_cols = columns;
        pcol = new unsigned int[n_cols+1];
        BCHKNR(pcol==NULL, CSCMatErr, Out of memory, n_cols);
        irow = new unsigned int[nnzeros];
        BCHKNR(irow==NULL, CSCMatErr, Out of memory, nnzeros);
        comp = new T[nnzeros];
        BCHKNR(comp==NULL, CSCMatErr, Out of memory, nnzeros);

        unsigned int i;
        for (i = 0; i < n_cols+1; ++i)
                pcol[i] = 0;
        for (i = 0; i < nnzeros;  ++i)
                irow[i] = 0;

        n_size = 0;
        n_max_size = nnzeros;
}


// delete memory space
// used by "=" operator and destructor
template <typename T>
void CSCMatrix<T>::destroy ()
{
        if (irow != NULL) 
                delete[] irow;
        if (comp != NULL) 
                delete[] comp;
        if (pcol != NULL) 
                delete[] pcol;

        n_size = 0;
        n_max_size = 0;
}


// copy to new memory space
// used by "=" operator and copy constructor
template <typename T>
void CSCMatrix<T>::copy (const CSCMatrix<T>& m)
{
        n_rows = m.n_rows;
        n_cols = m.n_cols;
        n_size = m.n_size;
        n_max_size = m.n_size;
        pcol = new unsigned int[n_cols+1];
        BCHKNR(pcol==NULL, CSCMatErr, Out of memory, n_cols);
        irow = new unsigned int[n_size];
        BCHKNR(pcol==NULL, CSCMatErr, Out of memory, n_size);
        comp = new T[n_size];
        BCHKNR(comp==NULL, CSCMatErr, Out of memory, n_size);

        for (unsigned int j = 0; j < n_cols+1; ++j)
                pcol[j] = m.pcol[j];
        for (unsigned int i = 0; i < n_size;   ++i) {
                irow[i] = m.irow[i];
                comp[i] = m.comp[i];
        }
}



// set all elements to a certain value
template <typename T>
inline CSCMatrix<T>& CSCMatrix<T>::fill (const T& val)
{
        if (val == (T)0)
                return this->clear();
        //for (unsigned int i = 0; i < n_size; ++i)
        //      comp[i] = val;
        for (unsigned c = 0; c < columns(); ++c)
                for (unsigned r = 0; r < rows(); ++r)
                        (*this)(r, c) = val;
        return *this;
}

// clear all elements
template <typename T>
inline CSCMatrix<T>& CSCMatrix<T>::clear()
{ CSTD__ memset (comp, 0, sizeof(T) * n_size); return *this; }

template <typename T>
CSCMatrix<T>& CSCMatrix<T>::setunit (const T& val)
{
        BCHK (n_rows != n_cols, CSCMatErr, 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>
CSCMatrix<T>& CSCMatrix<T>::swap(CSCMatrix<T>& m)
{       
        SWAP(n_rows, m.n_rows); SWAP(n_cols, m.n_cols);
        SWAP(n_size, m.n_size); SWAP(n_max_size, m.n_max_size);
        SWAP(pcol, m.pcol); SWAP(irow, m.irow);
        SWAP(comp, m.comp);
        return *this;
}

template <typename T>
CSCMatrix<T> CSCMatrix<T>::transposed_copy() const
{
        CSCMatrix<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 CSCMatrix<T>& CSCMatrix<T>::transpose()
{
        CSCMatrix<T> tr(this->transposed_copy());
        return swap(tr);
}

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

NAMESPACE_END

#endif  /* TBCI_CSCMATRIX_H */