00001 #ifndef MATRIXBD_H
00002 #define MATRIXBD_H
00003 00004 00005 00006 00007 00008 00009 00010 00011 00012 00013 00014 00015 00016 00017 00018 00019 00020 00021 00022 00023 00024 00025 00026 00027 00028 00029 00030 00031 00032
00033
00034
00035 #include "matrix.h"
00043 #define MAT_CB(M,NR,I,J,KU) (M[(J)*((NR)-1)+(I)+(KU)])
00044 #define REF_CB(A,I,J) MAT_CB(((A).D)->m,((A).D)->nrow,I,J,((A).D)->ku)
00045 #define THIS_CB(I,J) REF_CB(*this,I,J)
00046
00049
00050 #ifdef USE_NAMESPACE_LINAL
00051 namespace LinAl {
00052 #endif
00053
00068 class MatrixBD:public MatrixC<double>
00069 {
00070 public:
00072
00073
00074 explicit MatrixBD() : MatrixC<double>() {}
00076 explicit MatrixBD(const int nr, const int nc, const int nu, const double value=0.0)
00077 : MatrixC<double>(nr,nc,nu,value) {}
00079 explicit MatrixBD(const int nr, const int nc=1, const double value=0.0)
00080 : MatrixC<double>(nr,nc,0,value) {}
00082 MatrixBD(const MatrixBD& B) : MatrixC<double>(B) {} ;
00084 ~MatrixBD() {} ;
00086 MatrixBD& operator= (const double value)
00087 {this->MatrixC<double>::operator=(value); return *this;}
00089
00090
00091
00092 int nr() const ;
00094 int nr(const int j) const ;
00096 int nc() const ;
00098 int nc(const int i) const ;
00100 void resize(const int nr, const int nc, const double value=0) ;
00102 virtual Ref<double> operator() (const int i, const int j=0) {return MatrixC<double>::operator()(i,j);}
00104 virtual double operator() (const int i, const int j=0) const {return (THIS_CB(i,j));}
00106
00107
00108
00109 int row_begin(const int j) const {int ja=j-D->ku;return (ja>0?ja:0);}
00111 int row_end(const int j) const {int je=j+D->kl+1;return (je<D->ncol?je:D->ncol);}
00113 int col_begin(const int i) const {int ia=i-D->kl;return (ia>0?ia:0);}
00115 int col_end(const int i) const {int ie=i+D->ku+1;return (ie<D->ncol?ie:D->ncol);}
00117 int ku() const {return D->ku;}
00119 int kl() const {return D->kl;}
00121 protected:
00123 double* ref(const int i, const int j) const {assert(i>=j-D->ku);assert((i-j)<D->nrow);return &THIS_CB(i,j);}
00125 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);}
00126
00127 00129 00131 00133 00135 00137
00138 public:
00140 virtual void stream(ostream& s);
00145
00146 MatrixBD transpose() const ;
00147 00148 00149 00150 00151 00152
00154 void multiply(const MatrixBD& A, const MatrixBD& B,
00155 const char TRANSA='N', const char TRANSB='N')
00156 {
00157 cout <<
00158 "Multiplication of two symmetric matrices destroys the symmetry. \n" <<
00159 "Therefore it doesnt quite make sense. Instead multiplication of \n" <<
00160 "a symmetric matrix with a full matrix is implemented as a member\n" <<
00161 "function of the full matrix class.\n" ;
00162 }
00164 void multiply(const VectorD& x, VectorD& y, const char TRANS='N',
00165 const double ALPHA=1.0, const double BETA=0.0) const;
00167 void multiply(const double value);
00169 void LUdecompose(int* IPIV) ;
00175 void LUsubstitute(MatrixC<double>& B, int* IPIV) ;
00177 void LUsolve(MatrixC<double>& B) ;
00179 void LUinvert()
00180 {
00181 cout <<
00182 "Inversion of a banded matrix does in general not preserve the\n" <<
00183 "banded structure of the matrix. It is therefore not supported\n" <<
00184 "in this form.\n" ;
00185 }
00187 void CHdecompose() ;
00189 void CHsubstitute(MatrixC<double>& B) ;
00191 void CHsolve(MatrixC<double>& B) ;
00193 void CHinvert()
00194 {
00195 cout <<
00196 "Inversion of a banded matrix does in general not preserve the\n" <<
00197 "banded structure of the matrix. It is therefore not supported\n" <<
00198 "in this form.\n" ;
00199 }
00200 };
00201 #ifdef USE_NAMESPACE_LINAL
00202 }
00203 #endif
00204 #endif // #ifndef MATRIXBD_H