f_bandmatrix.h

Go to the documentation of this file.

//----------------------------------------------------------
// version 1.1, june 1998
//
// Marco Wiedenhaus
// based on tbandmatrixcomplex.h by Jens Lenge
//
// $Id: f_bandmatrix.h,v 1.11.2.19 2019/05/28 11:13:02 garloff Exp $
//----------------------------------------------------------

#ifndef TBCI_F_BANDMATRIX_H
#define TBCI_F_BANDMATRIX_H

//#include <stdlib.h> // for use of null macro

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

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

NAMESPACE_TBCI

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

class F_BandMatErr : public NumErr
{
 public:
        F_BandMatErr() : NumErr ("Error in F_BandMatrix library") {}
        F_BandMatErr(const char* t, const long i = 0) : NumErr (t, i) {}
};
#endif

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


template <typename T> class F_BandMatrix;
template <typename T> class Matrix;

//-----------------
// F_BandMatrix<T>
//-----------------

template <typename T>
class F_BandMatrix : public Matrix_Sig<T>
{
    public:
        typedef T value_type;
        typedef T element_type;
        typedef T aligned_value_type TALIGN(MIN_ALIGN2);
        // default constructor
        F_BandMatrix()  { allocate(1,0,0); }
        // constructor
        F_BandMatrix(unsigned int dimension,
                     unsigned int superDiags = 0,
                     unsigned int subDiags = 0)  
        { allocate(dimension, superDiags, subDiags); }
        // copy constructor
        F_BandMatrix(const F_BandMatrix<T>& m)  { copy(m); }
        F_BandMatrix(const Matrix<T>& m);
        // destructor
        ~F_BandMatrix()  { destroy(); }
        void clear();

        // element access (read)
        inline T operator() (unsigned int row, unsigned int column) const;
        // element access (write)
        inline T& operator() (unsigned int row, unsigned int column);
        // element access (write)
        inline void setval(const T z, unsigned int row, unsigned int column);
        T& setval (unsigned r, unsigned c) { return this->operator () (r, c); }
        
        // columnwise access
        TVector<T> get_col(unsigned int column) const;
        void set_col( const Vector<T>& v, unsigned int column);
        void set_col( const TVector<T>& tv, unsigned int column);
        void set_col( const TSVector<T>& tsv, unsigned int column);
        
        // query matrix dimensions
        unsigned int rows() const { return dim; }
        unsigned int columns() const { return dim; }
        unsigned int size() const { return dim; }
        unsigned int numSuper() const { return super; }
        unsigned int numSub() const { return sub; }
        unsigned int ldab() const { return sub+super+1; }
        // change matrix dimensions
        inline void resize(unsigned int newDim, unsigned int newSuper, unsigned int newSub);
        
        // matrix-matrix assignment
        F_BandMatrix<T>& operator= (const F_BandMatrix<T>& m);

        // matrix-matrix comparison (==)
         friend bool operator== FGD (const F_BandMatrix<T>&, const F_BandMatrix<T>&);
        // matrix-matrix comparison (!=)
         friend bool operator!= FGD (const F_BandMatrix<T>&, const F_BandMatrix<T>&);

        // matrix negation
         friend F_BandMatrix<T> operator- FGD (const F_BandMatrix<T>& m);
        
        // matrix-vector multiplication 
         friend TVector<T> do_fbdmat_vec_mul FGD (const F_BandMatrix<T>& m, const Vector<T>& v);
        // 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;

        // matrix-scalar multiplication
         friend F_BandMatrix<T> do_fbdmat_scale FGD (const F_BandMatrix<T>&m, const T z);
         friend F_BandMatrix<T> do_fbdmat_scale FGD (const T z, const F_BandMatrix<T>&m);
        // fast in-place matrix-scalar multiplication
        F_BandMatrix<T>& operator*= (const T z);
        // matrix-scalar division
         friend F_BandMatrix<T> operator/ FGD (const F_BandMatrix<T>& m, T const z);
        // fast in-place matrix-scalar division
        F_BandMatrix<T>& operator/= (const T z);
        // matrix addition and subtraction
         friend F_BandMatrix<T> operator+ FGD (const F_BandMatrix<T>&,const F_BandMatrix<T>&);
         friend F_BandMatrix<T> operator- FGD (const F_BandMatrix<T>&,const F_BandMatrix<T>&);
        
        F_BandMatrix<T>  transposed_copy() const;
        F_BandMatrix<T>& transpose();

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

        T* const & get_fortran_matrix() const  { return comp; }
   
        F_BandMatrix<T>& swap(F_BandMatrix<T>& m);

    protected:
        unsigned int dim, super, sub;
        T* comp;
        mutable T dummy;
        // used by constructor and redim()
        void allocate(unsigned int dimension,
                      unsigned int superDiags,
                      unsigned int subDiags);
        // used by destructor and "=" operator
        void destroy();
        // used by copy constr. and "="
        void copy(const F_BandMatrix<T>& m);
        // helpers
        void find_super(const Matrix<T>& m);
        void find_sub(const Matrix<T>& m);
};

//------------------------------------------------
// F_BandMatrix<T> inline function definitions
//------------------------------------------------

// element access (read)
template <typename T>
inline T F_BandMatrix<T>::operator() (unsigned int row, unsigned int column) const
{
        BCHK(row >= dim || column >= dim, F_BandMatErr, Illegal index, (row<column ? row : column), 0);
        if (column < row)
                if (row-column > sub)
                        return 0;
        if (column > row)
                if (column-row > super)
                        return 0;
        return comp[(super+row-column) + column*(super+sub+1)];
}

// element access (write)
template <typename T>
inline T& F_BandMatrix<T>::operator() (unsigned int row, unsigned int column)
{
        BCHK(row >= dim || column >= dim, F_BandMatErr, Illegal index, (row<column ? row : column), dummy=(T)0);
        BCHK(column < row && row-column > sub, F_BandMatErr, Illegal index, (row<column ? row : column), dummy=(T)0);
        BCHK(column > row && column-row > super, F_BandMatErr, Illegal index, (row<column ? row : column), dummy=(T)0);
        return comp[super+row-column + column*(super+sub+1)];
}

// element access (write)
// compatible with sparse matrix classes
template <typename T>
inline void F_BandMatrix<T>::setval(const T z, unsigned int row, unsigned int column)
{
        BCHKNR(row >= dim || column >= dim, F_BandMatErr, Illegal index, (row<column ? row : column));
        BCHKNR(column < row && row-column > sub, F_BandMatErr, Illegal index, (row<column ? row : column));
        BCHKNR(column > row && column-row > super, F_BandMatErr, Illegal index, (row<column ? row : column));
        comp[super+row-column + column*(super+sub+1)] = z;
}

// change matrix dimensions
template <typename T>
inline void F_BandMatrix<T>::resize(unsigned int newDim, unsigned int newSuper,
                                       unsigned int newSub)
{
        destroy();
        allocate(newDim, newSuper, newSub);
}

// extract a certain column
template <typename T>
TVector<T> F_BandMatrix<T>::get_col(unsigned int column) const
{
        TVector<T> ColVec(dim);
        BCHK(column >= dim, F_BandMatErr, Illegal index, column, ColVec);
        for (unsigned int i=0; i<dim; i++)
                ColVec.set((*this)(i,column),i);
        return ColVec;
}

// replace a certain column
template <typename T>
void F_BandMatrix<T>::set_col( const Vector<T>& v, unsigned int column)
{
        BCHKNR(column >= dim, F_BandMatErr, Illegal index, column);
        for (unsigned int i=0; i<dim; ++i) {
                        if (column <= i) 
                                if (!(i-column > sub))
                        (*this)(i,column) = v(i);
                        if (column > i)
                                if (!(column-i > super))
                                        (*this)(i,column) = v(i);

        }
}

template <typename T>
inline void F_BandMatrix<T>::set_col( const TVector<T>& tv, unsigned int column)
{
        Vector<T> v(tv);
        set_col(v,column);
}

template <typename T>
inline void F_BandMatrix<T>::set_col( const TSVector<T>& tsv, unsigned int column)
{
        Vector<T> v(tsv);
        set_col(v,column);
}

// matrix-matrix assignment
template <typename T>
F_BandMatrix<T>&        F_BandMatrix<T>::operator= (const F_BandMatrix<T>& m)
{
        BCHK(dim != m.dim, F_BandMatErr, Different sizes in assignment, m.dim, *this );
        if (this == &m) return *this;
        if (dim==m.dim && super==m.super && sub==m.sub) // perform fast in-place copy
        {
                CSTD__ memcpy(comp, m.comp, sizeof(T)*(dim*(super+sub+1)));
                return *this;
        }
        destroy();
        copy(m);
        return *this;
}

// matrix-matrix comparison (==)
// also works for matrices with different super/sub
template <typename T>
bool operator== (const F_BandMatrix<T>& m1, const F_BandMatrix<T>& m2)
{
        if (m1.dim != m2.dim) return false;
        for (unsigned int j=0; j<m1.dim; j++)
                for (unsigned int i=0; i<m1.dim; i++)
                        if (m1(i,j) != m2(i,j)) return false;
        return true;
}

// matrix-matrix comparison (!=)
// also works for matrices with different super/sub
template <typename T>
bool operator!= (const F_BandMatrix<T>& m1, const F_BandMatrix<T>& m2)
{
        if (m1.dim != m2.dim) return true;
        for (unsigned int j=0; j<m1.dim; j++)
                for (unsigned int i=0; i<m1.dim; i++)
                        if (m1(i,j) != m2(i,j)) return true;
        return false;
}

// matrix negation
template <typename T>
F_BandMatrix<T> operator- (const F_BandMatrix<T>& m)
{
        F_BandMatrix<T> res(m.dim, m.super, m.sub);
        for (unsigned int j=0; j<(m.dim*(m.sub+m.super+1)); ++j)
                res.comp[j] = -m.comp[j];
        return res; // JENS: avoid copy constructor soon
}

// matrix-vector multiplication
template <typename T>
TVector<T> do_fbdmat_vec_mul(const F_BandMatrix<T>&m, const Vector<T>& v)
// matrix-vector multiplication
{
        TVector<T> res(m.dim);
        BCHK(m.dim != v.size(), F_BandMatErr, Dimension conflict, v.size(), res);
        for (unsigned int i=0; i<m.dim; i++) {
                // left band limit:
                unsigned int left = (i>m.sub ? i-m.sub : 0);
                // right band limit:
                unsigned int right = MIN(i+m.super+1,m.dim);
                T sum(0.0);
                for (unsigned int j=left; j<right; j++)
                        sum += m.comp[(m.super+i-j)+j*(m.super+m.sub+1)] * v(j);
                res.set(sum,i);
        }
        return res; // JENS: avoid copy constructor soon!
}

template <typename T>
inline TVector<T> operator* (const F_BandMatrix<T>& m, const Vector<T>& v)
{
        return do_fbdmat_vec_mul(m, v);
}


template <typename T>
inline TVector<T> operator* (const F_BandMatrix<T>& m, const TVector<T>& tv)
{
        Vector<T> v(tv);
        return do_fbdmat_vec_mul(m, v);
}

// could be optimized ...
template <typename T>
inline TVector<T> operator* (const F_BandMatrix<T>& m, const TSVector<T>& tsv)
{
        Vector<T> v(tsv);
        return do_fbdmat_vec_mul(m, v);
}

// transpose-vector multiplication
template <typename T>
TVector<T> F_BandMatrix<T>::transMult(const Vector<T>& v) const
{
        TVector<T> res(dim);
        BCHK(dim != v.size(), F_BandMatErr, Dimension conflict, v.size(), res);
        for (unsigned int j=0; j<dim; j++)
        {
                // transpose left band limit:
                unsigned int left = (j>super ? j-super : 0);
                // transpose right band limit:
                unsigned int right = MIN(j+sub+1,dim);
                T sum(0.0);
                for (unsigned int i=left; i<right; i++)
                        sum += comp[(super+i-j)+j*(super+sub+1)] * v(i);
                res.set(sum,j);
        }
        return res; // JENS: avoid copy constructor soon!
}

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

// Could be optimized ...
template <typename T>
inline TVector<T> F_BandMatrix<T>::transMult(const TSVector<T>& tsv) const
{
        Vector<T> v(tsv);
        return(this->transMult(v));
}

// matrix-scalar multiplication
template <typename T>
F_BandMatrix<T> do_fbdmat_scale (const F_BandMatrix<T>& m, const T z)
{
        F_BandMatrix<T> res(m.dim, m.super, m.sub);
        for (unsigned int j=0; j<(m.dim*(m.super+m.sub+1)); ++j)
                res.comp[j] = m.comp[j] * z;
        return res; // JENS: avoid copy constructor soon!
}


// matrix-scalar multiplication
template <typename T>
F_BandMatrix<T> do_fbdmat_scale (const T z, const F_BandMatrix<T>& m)
{
        F_BandMatrix<T> res(m.dim, m.super, m.sub);
        for (unsigned int j=0; j<(m.dim*(m.super+m.sub+1)); ++j)
                res.comp[j] = z * m.comp[j];
        return res; // JENS: avoid copy constructor soon!
}

template <typename T>
inline F_BandMatrix<T> operator* (const F_BandMatrix<T>& m, const T z)
{
        return do_fbdmat_scale(m, z);
}

template <typename T>
inline F_BandMatrix<T> operator* (const T z, const F_BandMatrix<T>& m)
{
        return do_fbdmat_scale(z, m);
}

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


// matrix-scalar division
template <typename T>
F_BandMatrix<T> operator/ (const F_BandMatrix<T>& m, const T z)
{
        F_BandMatrix<T> res(m.dim, m.super, m.sub);
        for (unsigned int j=0; j<(m.dim*(m.super+m.sub+1)); ++j)
                res.comp[j] = m.comp[j] / z;
        return res; // JENS: avoid copy constructor soon!
}

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


template <typename T>
inline void F_BandMatrix<T>::find_super(const Matrix<T>& m)
{
        for (super = dim-1; super > 0; --super)
                for (unsigned r = 0; r < dim-super; ++r)
                        if (m(r, r+super) != (T)0)
                                return;
}

template <typename T>
inline void F_BandMatrix<T>::find_sub(const Matrix<T>& m)
{
        for (sub = dim-1; sub > 0; --sub)
                for (unsigned c = 0; c < dim-sub; ++c)
                        if (m(c+sub, c) != (T)0)
                                return;
}

template <typename T>
F_BandMatrix<T>::F_BandMatrix(const Matrix<T>& m)
{
        BCHKNR(m.rows() != m.columns(), F_BandMatErr, only square matrices possible, m.rows());
        dim = m.rows();
        this->find_super(m); this->find_sub(m);
        // Allocate and copy
        this->allocate(dim, super, sub);
        for (unsigned r = 0; r < dim; ++r)
                for (unsigned c = MAX(0, (int)(r-sub));
                                c < MIN(dim, r+super+1); ++c)
                        (*this)(r,c) = m(r,c);
}

#ifndef LAPACK_INLINE
# define LAPACK_INLINE  
#endif
// allocate memory space
// used by constructor and redim()
template <typename T>
LAPACK_INLINE void F_BandMatrix<T>::allocate(unsigned int dimension,
                                unsigned int superDiags,
                                unsigned int subDiags)
{
        BCHKNR(dimension == 0, F_BandMatErr, Zero space allocating, 0);
        BCHKNR(superDiags >= dimension || subDiags >= dimension, F_BandMatErr, Dimension conflict, (subDiags<superDiags? subDiags : superDiags));
        dim = dimension;
        super = superDiags;
        sub = subDiags;
        unsigned int colSize = super + sub + 1;
                comp = new T[colSize*dim];
        BCHKNR(comp==NULL, F_BandMatErr, Out of memory, colSize*dim);
//      for(int i=0; i<(colSize*dim); i++) comp[i] = NULL;      
}

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

// copy to new memory space
// used by "=" operator and copy constructor
template <typename T>
void F_BandMatrix<T>::copy(const F_BandMatrix<T>& m)
{
        dim = m.dim;
        super = m.super;
        sub = m.sub;
        unsigned int colSize = super + sub + 1;
        comp = new T[colSize*dim];
        BCHKNR(comp==NULL, F_BandMatErr, Out of memory, colSize*dim);
        CSTD__ memcpy(comp, m.comp, sizeof(T)*(dim*colSize));
}


// matrix addition and subtraction
template <typename T>
F_BandMatrix<T> operator+ (const F_BandMatrix<T>& A,const F_BandMatrix<T>& B)
{
        unsigned int numSub(0),numSuper(0),dim(0);
        // calculate the dimations of C
        dim = MAX(A.size(), B.size());
        numSuper = MAX(A.numSuper(), B.numSuper());
        numSub = MAX(A.numSub(), B.numSub());

        // create C and full C
        F_BandMatrix<T> C(dim,numSuper,numSub);

        BCHK(A.size()!=B.size(), F_BandMatErr, Dimension conflict, A.size(), C );

        for(unsigned int i(0);i<C.size();i++)
        {
                for(int j(i-C.numSub());j<=(int)(i+C.numSuper());j++)
                {
                        if((j>=0)&&(j<(int)C.size()))
                        {
                                if( (j>=(int)(i-A.numSub())) && (j<=(int)(i+A.numSuper()))
                                &&  (j>=(int)(i-B.numSub())) && (j<=(int)(i+B.numSuper())) )
                                        C(i,j)=A(i,j)+B(i,j);
                                else if( (j>=(int)(i-B.numSub())) && (j<=(int)(i+B.numSuper())) )
                                        C(i,j)=B(i,j);
                                else if( (j>=(int)(i-A.numSub())) && (j<=(int)(i+A.numSuper())) )
                                        C(i,j)=A(i,j);                  
                        }
                }
        }
        return C;
} 

template <typename T>
F_BandMatrix<T> operator- (const F_BandMatrix<T>& A, const F_BandMatrix<T>& B)
{
        // calculate the dimensions of C
        const unsigned dim = MAX(A.size(), B.size());
        const unsigned numSuper = MAX(A.numSuper(), B.numSuper());
        const unsigned numSub = MAX(A.numSub(), B.numSub());

        // create C and fill C
        F_BandMatrix<T> C(dim, numSuper, numSub);

        BCHK(A.size()!=B.size(), F_BandMatErr, Dimension conflict, A.size(), C );

        for(unsigned int i(0);i<C.size();i++)
        {
                for(int j(i-C.numSub());j<=(int)(i+C.numSuper());j++)
                {
                        if((j>=0)&&(j<(int)C.size()))
                        {
                                if( (j>=(int)(i-A.numSub())) && (j<=(int)(i+A.numSuper()))
                                &&  (j>=(int)(i-B.numSub())) && (j<=(int)(i+B.numSuper())) )
                                        C(i,j)=A(i,j)-B(i,j);
                                else if( (j>=(int)(i-B.numSub())) && (j<=(int)(i+B.numSuper())) )
                                        C(i,j)=-B(i,j);
                                else if( (j>=(int)(i-A.numSub())) && (j<=(int)(i+A.numSuper())) )
                                        C(i,j)=A(i,j);                  
                        }
                }
        }
        return C;
}

template <typename T>
STD__ ostream& operator<< (STD__ ostream& stream, const F_BandMatrix<T>& m)
{
        for (unsigned int i = 0; i < m.rows(); i++)
        {
                for (unsigned int j = 0; j < m.columns(); j++)
                        stream << m(i,j) << " ";
                stream << STD__ endl;
        }
        return stream;
}


template <typename T>
LAPACK_INLINE void F_BandMatrix<T>::clear()
{
        for (unsigned int j=0; j<(dim*(super+sub+1)); j++)
                comp[j] = (T) 0;
}


template <typename T>
F_BandMatrix<T> F_BandMatrix<T>::transposed_copy() const
{
        F_BandMatrix<T> f(dim, sub, super);
        for (unsigned int r = 0; r < dim; ++r)
                for (unsigned int c = MAX(0, (int)(r-super)); 
                                  c < MIN(dim, r+sub+1); ++c)
                        f(r,c) = (*this)(c,r);
        return f;
}

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


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

template <typename T>
F_BandMatrix<T>& F_BandMatrix<T>::swap(F_BandMatrix<T>& m)
{
        SWAP(dim, m.dim);
        SWAP(super, m.super); SWAP(sub, m.sub);
        SWAP(comp, m.comp);
        return *this;
}

NAMESPACE_END

#endif  /* TBCI_F_BANDMATRIX_H */