---

matrixBD.h

Go to the documentation of this file.

#ifndef MATRIXBD_H
#define MATRIXBD_H
/*------------------------------------------------------------------------------
 * matrixBD.h   A simple, easy to use and efficient Matrix library.
 *              Definition and partial implementation of class MatrixBD.
 *              
 *              Copyright 1997,1998,1999,2000 Roland Krause
 * -----------------------------------------------------------------------------
 * RCS/CVS history of changes
 *  
 * $Log: matrixBD.h,v $
 * Revision 1.3  2000/11/09 18:41:31  rokrau
 * The first commit to the sourceforge repostory. Update of all
 * intermediate changes makes this the official repository from now on.
 *
 * Revision 1.3  2000/11/09 17:51:53  rokrau
 * Removed one of the redundant multiplication methods.
 *
 * Revision 1.2  2000/10/27 17:08:09  rokrau
 * some fixes for integer matrices
 *
 * Revision 1.1  2000/10/24 06:03:04  rokrau
 * Split header files, made Ref constructor private w/ some drastic consequences for
 * Ref operators. A lot of minor changes.
 *
 * -----------------------------------------------------------------------------
 * $Date: 2000/11/09 18:41:31 $
 * $Revision: 1.3 $
 * $Author: rokrau $
 * $State: Exp $
 * -----------------------------------------------------------------------------
 */
// standard includes
// project includes
#include    "matrix.h"
#define MAT_CB(M,NR,I,J,KU) (M[(J)*((NR)-1)+(I)+(KU)])
#define REF_CB(A,I,J) MAT_CB(((A).D)->m,((A).D)->nrow,I,J,((A).D)->ku)
#define THIS_CB(I,J) REF_CB(*this,I,J)

//using namespace std;
#ifdef USE_NAMESPACE_LINAL
namespace LinAl {
#endif

class MatrixBD:public MatrixC<double>
{
public:


    explicit MatrixBD() : MatrixC<double>() {}
    explicit MatrixBD(const int nr, const int nc, const int nu, const double value=0.0)
    : MatrixC<double>(nr,nc,nu,value) {}
    explicit MatrixBD(const int nr, const int nc=1, const double value=0.0)
    : MatrixC<double>(nr,nc,0,value) {}
    MatrixBD(const MatrixBD& B) : MatrixC<double>(B) {} ;
    ~MatrixBD() {} ;
    MatrixBD& operator= (const double value)
        {this->MatrixC<double>::operator=(value); return *this;}



    int nr() const ;
    int nr(const int j) const ;
    int nc() const ;
    int nc(const int i) const ;
    void resize(const int nr, const int nc, const double value=0) ;
    virtual Ref<double> operator() (const int i, const int j=0) {return MatrixC<double>::operator()(i,j);}
    virtual double operator() (const int i, const int j=0) const {return (THIS_CB(i,j));}



    int row_begin(const int j) const {int ja=j-D->ku;return (ja>0?ja:0);}
    int row_end(const int j) const {int je=j+D->kl+1;return (je<D->ncol?je:D->ncol);}
    int col_begin(const int i) const {int ia=i-D->kl;return (ia>0?ia:0);}
    int col_end(const int i) const {int ie=i+D->ku+1;return (ie<D->ncol?ie:D->ncol);}
    int ku() const {return D->ku;}
    int kl() const {return D->kl;}
protected:
    double* ref(const int i, const int j) const {assert(i>=j-D->ku);assert((i-j)<D->nrow);return &THIS_CB(i,j);}
    double read(const int i, const int j) const {assert(i>=j-D->ku);return (((i-j)<D->nrow)?THIS_CB(i,j):0.0);}
    // need to throw an exception to warn user that he cant write
    // outside the band of the matrix
    void write(const int i, const int j, const double v) {assert(i>=j-D->ku);if ((i-j)<D->nrow) {D=D->get_own_copy();THIS_CB(i,j)=v;}}
    void add(const int i, const int j, const double v) {assert(i>=j-D->ku);if ((i-j)<D->nrow) {D=D->get_own_copy();THIS_CB(i,j)+=v;}}
    void sub(const int i, const int j, const double v) {assert(i>=j-D->ku);if ((i-j)<D->nrow) {D=D->get_own_copy();THIS_CB(i,j)-=v;}}
    void mul(const int i, const int j, const double v) {assert(i>=j-D->ku);if ((i-j)<D->nrow) {D=D->get_own_copy();THIS_CB(i,j)*=v;}}
    void div(const int i, const int j, const double v) {assert(i>=j-D->ku);if ((i-j)<D->nrow) {D=D->get_own_copy();THIS_CB(i,j)/=v;}}
public:
    virtual void stream(ostream& s);

    MatrixBD transpose() const ;
    /*
     * Multiplication of two symmetric matrices destroys the symmetry.
     * Therefore it doesnt quite make sense. Instead multiplication of a
     * symmetric matrix with a full matrix can be implemented as a member
     * function of the full matrix class.
     */
    void multiply(const MatrixBD& A, const MatrixBD& B,
                  const char TRANSA='N', const char TRANSB='N')
    {
    cout <<
        "Multiplication of two symmetric matrices destroys the symmetry. \n" <<
        "Therefore it doesnt quite make sense. Instead multiplication of \n" <<
        "a symmetric matrix with a full matrix is implemented as a member\n" <<
        "function of the full matrix class.\n" ;
    }
    void multiply(const VectorD& x, VectorD& y, const char TRANS='N',
                  const double ALPHA=1.0, const double BETA=0.0) const;
    void multiply(const double value);
    void LUdecompose(int* IPIV) ;
    void LUsubstitute(MatrixC<double>& B, int* IPIV) ;
    void LUsolve(MatrixC<double>& B) ;
    void LUinvert()
    {
    cout <<
        "Inversion of a banded matrix does in general not preserve the\n" <<
        "banded structure of the matrix. It is therefore not supported\n" <<
        "in this form.\n" ;
    }
    void CHdecompose() ;
    void CHsubstitute(MatrixC<double>& B) ;
    void CHsolve(MatrixC<double>& B) ;
    void CHinvert()
    {
    cout <<
        "Inversion of a banded matrix does in general not preserve the\n" <<
        "banded structure of the matrix. It is therefore not supported\n" <<
        "in this form.\n" ;
    }
};
#ifdef USE_NAMESPACE_LINAL
} // namespace LinAl
#endif
#endif // #ifndef MATRIXBD_H

---