matrix.h Source File

00001 #ifndef MATRIX_H
00002 #define MATRIX_H
00003 /*------------------------------------------------------------------------------
00004  * matrix.h     A simple, easy to use and efficient Matrix library.
00005  *              Definition of the pure virtual class Matrix. Definition and
00006  *              partial implementation of virtual base class MatrixC.
00007  *              
00008  *              Copyright 1997,1998,1999,2000 Roland Krause
00009  * -----------------------------------------------------------------------------
00010  * RCS/CVS history of changes 
00011  *  
00012  * $Log: matrix.h,v $
00013  * Revision 1.53  2000/11/21 17:03:07  rokrau
00014  * calls to BLAS and LAPACK use lowercase function names
00015  * if the preprocessor directive FTN_CASE_LOWER is given
00016  *
00017  * Revision 1.52  2000/11/09 18:41:31  rokrau
00018  * The first commit to the sourceforge repostory. Update of all
00019  * intermediate changes makes this the official repository from now on.
00020  *
00021  * Revision 1.52  2000/11/09 17:46:54  rokrau
00022  * More constructors and some minor updates for integer matrices
00023  *
00024  * Revision 1.51  2000/10/27 17:08:09  rokrau
00025  * some fixes for integer matrices
00026  *
00027  * Revision 1.50  2000/10/24 06:03:04  rokrau
00028  * Split header files, made Ref constructor private w/ some drastic consequences for
00029  * Ref operators. A lot of minor changes.
00030  *
00031  * -----------------------------------------------------------------------------
00032  * $Date: 2000/11/21 17:03:07 $
00033  * $Revision: 1.53 $
00034  * $Author: rokrau $
00035  * $State: Exp $
00036  * -----------------------------------------------------------------------------
00037  */
00038 // windows...
00039 #ifdef WIN32
00040 #define __STDCALL __cdecl
00041 #else
00042 #define __STDCALL
00043 #endif
00044 // standard includes
00045 #include    <iostream>
00046 #include    <cassert>
00047 #include    <cmath>
00048 #include    <cfloat>
00049 #include    <valarray>
00050 // constants
00052 #define OUTPUT_PREC 6
00054 #define EPS_EQ_NULL DBL_EPSILON
00062 #define MAT_CF(M,NR,I,J) (M[(J)*(NR)+(I)])
00063 #define REF_CF(A,I,J) MAT_CF(((A).D)->m,((A).D)->nrow,I,J)
00064 #define THIS_CF(I,J) REF_CF(*this,I,J)
00065 
00068 #ifdef USE_NAMESPACE_LINAL
00069 namespace LinAl {
00070 #endif
00071 /*
00072  * Classes
00073  */
00074 template <class T>class Ref;
00075 template <class T>class Matrix ;
00076 template <class T>class MdataC ;
00077 template <class T>class MatrixC ;
00078 class MatrixFD;
00079 class MatrixSD;
00080 class MatrixBD;
00081 class MatrixFI;
00082 class VectorD ;
00083 class VectorI ;
00101 template <class Tp> inline istream& operator>> (istream& s, Ref<Tp>& r) ;
00102 template <class Tp> inline ostream& operator>> (ostream& s, Ref<Tp>& r) ;
00103 template <class Tp>
00104 class Ref {
00105     friend istream& operator>> <> (istream& s, Ref<Tp>& r);
00106     friend ostream& operator>> <> (ostream& s, Ref<Tp>& r);
00107     friend class MatrixC<Tp>;
00108     friend class MatrixSP;
00109 public:
00111 
00112     operator Tp()  const {return M.read(i,j) ;}
00113     operator Tp*() const {return M.ref(i,j) ;}
00115 
00116 
00117     Tp operator=(const Tp v) const {M.write(i,j,v); return M.read(i,j);}
00118     Tp operator=(const Ref<Tp> R) const {M.write(i,j,R); return M.read(i,j);}
00120 
00121 
00122     Tp operator++() const {M.add(i,j,1) ;return M.read(i,j) ;}
00123     Tp operator--() const {M.sub(i,j,1) ;return M.read(i,j) ;}
00124     void operator+=(Tp v) const {M.add(i,j,v) ;}
00125     void operator-=(Tp v) const {M.sub(i,j,v) ;}
00126     void operator*=(Tp v) const {M.mul(i,j,v) ;}
00127     void operator/=(Tp v) const {M.div(i,j,v) ;}
00128     Tp operator+(const Tp v) const {return M.read(i,j)+v;}
00129     Tp operator-(const Tp v) const {return M.read(i,j)-v;}
00130     Tp operator*(const Tp v) const {return M.read(i,j)*v;}
00131     Tp operator/(const Tp v) const {return M.read(i,j)/v;}
00133 protected:
00134 private:    
00136     explicit Ref(Matrix<Tp>& _M, const int _i, const int _j=0) : M(_M),i(_i),j(_j) {}
00138     Matrix<Tp>& M ;
00140     int i ;
00142     int j ;
00143 };
00147 template <class Tp> inline Tp operator+ (const Tp& a, const Ref<Tp>& b) {Tp r;r=b+a;return r;}
00148 template <class Tp> inline Tp operator- (const Tp& a, const Ref<Tp>& b) {Tp r;r=b-a;return -r;}
00149 template <class Tp> inline Tp operator* (const Tp& a, const Ref<Tp>& b) {Tp r;r=b*a;return -r;}
00150 template <class Tp> inline Tp operator/ (const Tp& a, const Ref<Tp>& b) {Tp r;r=b/a;return 1.0f/r;}
00151 template <class Tp> inline istream& operator>> (istream& s, Ref<Tp>& r) {s>>r.M.write(r.i,r.j);return s;}
00152 template <class Tp> inline ostream& operator>> (ostream& s, Ref<Tp>& r) {s<<r.M.read(r.i,r.j);return s;}
00173 template <class Tp>
00174 class Matrix {
00175     friend Ref<Tp>;
00176 public:
00178     virtual ~Matrix() {} ;
00180     virtual int nr(const int i=0) const = 0 ;
00182     virtual int nc(const int i=0) const = 0 ;
00184     virtual int size() const = 0 ;
00186     virtual bool empty() const = 0 ;
00188     virtual Ref<Tp> operator() (const int i, const int j=0) = 0 ;
00190     virtual Tp operator() (const int i, const int j=0) const = 0 ;
00191 protected:
00193     virtual Tp* ref(const int i, const int j) const = 0 ;
00195     virtual Tp read(const int i, const int j) const = 0 ;
00197     virtual void write(const int i, const int j, const Tp v) = 0 ;
00199     virtual void add(const int i, const int j, const Tp v) = 0 ;
00201     virtual void sub(const int i, const int j, const Tp v) = 0 ;
00203     virtual void mul(const int i, const int j, const Tp v) = 0 ;
00205     virtual void div(const int i, const int j, const Tp v) = 0 ;
00206 };
00217 template <class Tp>
00218 class MdataC {
00219 public:
00221 
00222 
00223     explicit MdataC() : m(),nrow(0), ncol(0), ku(0), kl(0) {nref=1;}
00225     explicit MdataC(const int nr, const int nc, const Tp value)
00226     : m(value,nr*nc),nrow(nr), ncol(nc), ku(0), kl(0) {nref=1;}
00228     explicit MdataC(const int nr, const int nc, const int nu, const Tp value)
00229     : m(value,nr*nc) ,nrow(nr), ncol(nc), ku(nu) {kl=nr-nu-1;nref=1;}
00231     explicit MdataC(const Tp* a, const int nr, const int nc)
00232     : m(a,nr*nc), nrow(nr), ncol(nc), ku(0), kl(0) {nref=1;}
00234     explicit MdataC(const valarray<Tp>& a, const int nr, const int nc)
00235     : m(a), nrow(nr), ncol(nc), ku(0), kl(0) {nref=1;}
00237     ~MdataC() {}
00239 
00240 
00241 
00242     MdataC* get_own_copy();
00244     void assign(const Tp value);
00246     void assign(const int nr, const Tp* value);
00248     void assign(const int nr, const int nc, const Tp** value);
00250     void resize(const int nr, const int nc);
00252 
00253 
00254 
00255     valarray<Tp> m ;
00257     int nrow ;
00259     int ncol ;
00261     int nref ;
00263     int ku;
00265     int kl;
00267 private:
00269     explicit MdataC(const MdataC& D)
00270     : m(D.m), nrow(D.nrow), ncol(D.ncol), ku(D.ku), kl(D.kl) {nref=1;}
00272     MdataC& operator=(const MdataC&) ;
00273 };
00304 template <class Tp>
00305 class MatrixC:public Matrix<Tp> {
00306 public:
00308 
00309 
00310     explicit MatrixC() {D=new MdataC<Tp>();}
00312     explicit MatrixC(const int nr, const int nc, const Tp value=0)
00313         {assert((nr*nc)>0); D=new MdataC<Tp>(nr,nc,value);}
00315     explicit MatrixC(const int nr, const int nc, const int nu, const Tp value)
00316         {assert((nr*nc)>0); D=new MdataC<Tp>(nr,nc,nu,value);}
00318     explicit MatrixC(const Tp* A, const int nr, const int nc=1)
00319         {assert((nr*nc)>0); D=new MdataC<Tp>(A,nr,nc);}
00321     explicit MatrixC(const valarray<Tp>& A, const int nc=1)
00322         {assert(nc);int nr=A.size()/nc;D=new MdataC<Tp>(A,nr,nc);}
00324     MatrixC(const MatrixC& B)   {B.D->nref++ ; D=B.D ;}
00326     virtual ~MatrixC() {if (--D->nref==0) delete D;} ;
00328 
00329 
00330 
00331     MatrixC& operator=(const MatrixC& B);
00333     virtual MatrixC& operator=(const Tp value);
00335 
00339 
00340     virtual int nr(const int i=0) const {return D->nrow;} ;
00342     virtual int nc(const int i=0) const {return D->ncol;} ;
00344     virtual int size() const {return D->m.size();}
00346     virtual bool empty() const {return ((D->m.size()==0)?true:false);}
00348     virtual Ref<Tp> operator() (const int i, const int j=0) {return Ref<Tp>(*this,i,j);}
00350     virtual Tp operator() (const int i, const int j=0) const {return THIS_CF(i,j);}
00352 
00353 
00354 
00355     void get_own_copy() {D=D->get_own_copy();}
00357     void resize(const int nr, const int nc, const Tp value=0) ;
00359     virtual void stream(istream& s);
00361     virtual void stream(ostream& s);
00363     virtual MatrixC& assign(const int nc, const Tp* value);
00365     virtual MatrixC& assign(const int nr, const int nc, const Tp* value);
00367     virtual MatrixC& assign(const int nr, const int nc, const Tp** value);
00369 
00370     virtual int capacity() const {return D->m.size();}
00372     bool operator==(const MatrixC<Tp>& B) const ;
00374     bool operator!=(const MatrixC<Tp>& B) const ;
00383 
00384     operator Tp*() const {return &D->m[0];}
00386     operator Tp*() {return &D->m[0];}
00388 
00397 
00398     MatrixC<Tp> transpose() const ;
00400     void add(const MatrixC<Tp>& B);
00402     void add(const Tp value);
00404     void subtract(const MatrixC<Tp>& B);
00406     void subtract(const Tp value);
00408     void multiply(const MatrixC<Tp>& A, const MatrixC<Tp>& B,
00409                   const char TRANSA=0, const char TRANSB=0,
00410                   const Tp ALPHA=1, const Tp BETA=0);
00411     
00413     void multiply(const Tp value);
00415 protected:
00417 
00418 
00419     virtual Tp& operator[] (const int i) const {return D->m[i];}
00421     virtual Tp& operator[] (const int i) {return D->m[i];}
00423     virtual Tp* ref(const int i, const int j) const {return &THIS_CF(i,j);}
00425     virtual Tp read(const int i, const int j) const {return THIS_CF(i,j);}
00427     virtual void write(const int i, const int j, const Tp v) {D=D->get_own_copy();THIS_CF(i,j)=v;}
00429     virtual void add(const int i, const int j, const Tp v) {D=D->get_own_copy();THIS_CF(i,j)+=v;}
00431     virtual void sub(const int i, const int j, const Tp v) {D=D->get_own_copy();THIS_CF(i,j)-=v;}
00433     virtual void mul(const int i, const int j, const Tp v) {D=D->get_own_copy();THIS_CF(i,j)*=v;}
00435     virtual void div(const int i, const int j, const Tp v) {D=D->get_own_copy();THIS_CF(i,j)/=v;}
00437 
00438     MdataC<Tp>* D ;
00439 };
00444 template <class Tp>
00445 inline istream& operator>> (istream& s, MatrixC<Tp>& M) {M.stream(s);return s;}
00446 template <class Tp>
00447 inline ostream& operator<< (ostream& s, MatrixC<Tp>& M) {M.stream(s);return s;}
00449 
00453 #ifdef FTN_CASE_LOWER
00454 #define DSCAL_   dscal_
00455 #define DDOT_    ddot_
00456 #define DAXPY_   daxpy_
00457 #define DCOPY_   dcopy_
00458 #define DNRM2_   dnrm2_
00459 #define DGEMV_   dgemv_
00460 #define DSPMV_   dspmv_
00461 #define DSBMV_   dsbmv_
00462 #define DGBMV_   dgbmv_
00463 #define DGEMM_   dgemm_
00464 #define DSYMM_   dsymm_
00465 #define DGESV_   dgesv_
00466 #define DGETRF_  dgetrf_
00467 #define DGETRS_  dgetrs_
00468 #define DGETRI_  dgetri_
00469 #define DPOSV_   dposv_
00470 #define DPOTRF_  dpotrf_
00471 #define DPOTRS_  dpotrs_
00472 #define DPOTRI_  dpotri_
00473 #define DPPSV_   dppsv_
00474 #define DPPTRF_  dpptrf_
00475 #define DPPTRS_  dpptrs_
00476 #define DPPTRI_  dpptri_
00477 #define DPBSV_   dpbsv_
00478 #define DPBTRF_  dpbtrf_
00479 #define DPBTRS_  dpbtrs_
00480 #define DGBSV_   dgbsv_
00481 #define DGBTRF_  dgbtrf_
00482 #define DGBTRS_  dgbtrs_
00483 #endif
00484 extern "C" {
00485     void __STDCALL DSCAL_ (const int*,const double*,double*,const int*) ;
00486     double __STDCALL DDOT_ (const int*,const double*,const int*,double*,const int*) ;
00487     void __STDCALL DAXPY_ (const int*,const double*,const double*,const int*,double*,const int*) ;
00488     void __STDCALL DCOPY_ (const int*,const double*,const int*,double*,const int*) ;
00489     double __STDCALL DNRM2_(const int*,const double*,const int*) ;
00490     void __STDCALL DGEMV_ (const char*,const int*,const int*,const double*,const double*,const int*,const double*,const int*,const double*,const double*,const int*);
00491     void __STDCALL DSPMV_ (const char*,const int*,const double*,const double*,const double*,const int*,const double*,const double*,const int*);
00492     void __STDCALL DSBMV_ (const char*,const int*,const int*,const double*,const double*,const int*,const double*,const int*,const double*,const double*,const int*);
00493     void __STDCALL DGBMV_ (const char*,const int*,const int*,const int*,const int*,const double*,const double*,const int*,const double*,const int*,const double*,const double*,const int*);
00494     void __STDCALL DGEMM_ (const char*,const char*,const int*,const int*,const int*,const double*,const double*,const int*,const double*,const int*,const double*,double*,const int*) ;
00495     void __STDCALL DSYMM_ (const char*,const char*,const int*,const int*,const double*,const double*,const int*,const double*,const int*,const double*,const double*,const int*);
00496     void __STDCALL DGESV_ (const int*,const int*,double*,const int*,int*,double*,const int*,int*) ;
00497     void __STDCALL DGETRF_(const int*,const int*,double*,const int*,int*,int*);
00498     void __STDCALL DGETRS_(const char*,const int*,const int*,double*,const int*,int*,double*,const int*,int*) ;
00499     void __STDCALL DGETRI_(const int*,double*, const int*,int*,double*,int*,int*) ;
00500     void __STDCALL DPOSV_ (const char*,const int*,const int*,double*,const int*,double*,const int*,int*);
00501     void __STDCALL DPOTRF_(const char*,const int*,double*,const int*,int*);
00502     void __STDCALL DPOTRS_(const char*,const int*,const int*,double*,const int*,double*,const int*,int*);
00503     void __STDCALL DPOTRI_(const char*,const int*,double*,const int*,int*);
00504     void __STDCALL DPPSV_ (char*,int*,int*,double*,double*,int*,int*);
00505     void __STDCALL DPPTRF_(const char*,const int*,double*,int*);
00506     void __STDCALL DPPTRS_(const char*,const int*,const int*,double*,double*,const int*,int*);
00507     void __STDCALL DPPTRI_(const char*,const int*,double*,int*);
00508     void __STDCALL DPBSV_ (const char*,const int*,const int*,const int*,double*,const int*,double*,const int*,int*);
00509     void __STDCALL DPBTRF_(const char*,const int*,const int*,double*,const int*,int*);
00510     void __STDCALL DPBTRS_(const char*,const int*,const int*,const int*,double*,const int*,double*,const int*,int*);
00511     void __STDCALL DGBSV_   (const int*,const int*,const int*,const int*,double*,const int*,int*,double*,const int*,int*);
00512     void __STDCALL DGBTRF_(const int*,const int*,const int*,const int*,double*,const int*,int*,int*);
00513     void __STDCALL DGBTRS_(const char*,const int*,const int*,const int*,const int*,double*,const int*,int*,double*,const int*,int*);
00514 }
00515 #ifdef USE_NAMESPACE_LINAL
00516 }   // namespace LinAL
00517 #endif
00518 #endif  // #ifndef MATRIX_H