QMatrixC.hh

Go to the documentation of this file.

 
#ifndef __QMATRIXC_HH_
#define __QMATRIXC_HH_
 
 
#include <gsl/gsl_complex.h>
#include <gsl/gsl_matrix_complex_double.h>
#include <map>
#include <iostream>
#include <sstream>
#include <QObject.hh>
#include <QMatrix.hh>
#include <QComplex.hh>
 
Q_BEGIN_NAMESPACE
 
class QVectorC;
 
class QMatrixC : public QObject {
public:
  QMatrixC() ;
  explicit QMatrixC(const UInt_t nrow,const UInt_t ncol);
  explicit QMatrixC(const QVectorC& vec);
  QMatrixC(const QMatrixC& orig);
  QMatrixC(const QMatrix& orig);
  
  virtual ~QMatrixC();
  void Resize(const UInt_t nrow,const UInt_t ncol);
  UInt_t Size() const { return fSize/2; };
  
 
  void BuildTest(UInt_t nI,UInt_t nJ) {
    Resize(nI,nJ);
    for(UInt_t i=0;i<nI;i++){
      for(UInt_t j=i;j<nJ;j++){
        (*this)(i,j)=QComplex(10*(i+j)*(i+j)+1,3.*(1.*i-1.*j));
        (*this)(j,i)=(*this)(i,j).Conj();
        if(i==j)
          (*this)(i,j)+=1000;
      }
    }
  }
 
  void Clear();
  inline UInt_t GetNRow() const {return fNCol ? fSize/(fNCol*2) : 0;};
  inline UInt_t GetNCol() const {return fNCol;};
  void Initialize(const QComplex& val);
  void Initialize(const double re=0,const double im=0){Initialize(QComplex(re,im));};
  void SetToIdentity();
  void SetCol(const UInt_t ncol,const QVectorC& vec);
  void SetRow(const UInt_t nrow,const QVectorC& vec);
  //ma che e' ao
  //        void Set(const std::map<int,int>& p);
  QVectorC GetCol(const UInt_t ncol) const;
  QVectorC GetRow(const UInt_t nrow) const;
  QVectorC GetDiagonal() const;
  
  //        QVectorC GetRowByColumnIntValue(int col,int val) const;
  QComplex operator()(const UInt_t i,const UInt_t j);
  const QComplex operator() (const UInt_t i,const UInt_t j) const;
  const QMatrixC& operator=(const QMatrixC& orig);
  QMatrixC   operator-() const;
  const QMatrixC& operator*=(const QComplex& t);
  const QMatrixC& operator*=(const double t);
  const QMatrixC & operator*=(const QMatrixC& other);
  const QMatrixC & operator/=(const QComplex& t);
  const QMatrixC& operator/=(const double t);
  const QMatrixC & operator+=(const QMatrixC& mat);
  const QMatrixC & operator-=(const QMatrixC& mat);
  const QMatrixC& Mult(const QMatrixC& mat);
  const QMatrixC& Div(const QMatrixC& mat);
  
  QMatrixC operator*(const QMatrixC& mat) const;
  QMatrixC operator*(const QComplex& t) const;
  QMatrixC operator*(const double t) const;
  QMatrixC operator/(const QComplex& t) const;
  QMatrixC operator/(const double t) const;
  QMatrixC operator+(const QMatrixC& mat) const;
  QMatrixC operator-(const QMatrixC& mat) const;
  QVectorC operator*(const QVectorC& vec) const;
  QMatrixC T() const; // returns a different matrix transposed wrt the original
  const QMatrixC &Transpose(); // transposes the matrix itself
  const QMatrixC & H();
  QMatrixC GetH() const;
  const QMatrixC &Conjugate();
  
  QMatrixC GetConjugate() const;
  
  QMatrix Magnitude() const;
 
  QMatrixC Inv() const; // returns a different matrix inverted wrt the original
  const QMatrixC &Invert(); // inverts the matrix itself
  QComplex Det() const;
 
  QVector HermitianEigenValues();
 
  void HermitianEigenSystem(QVector &Val,QMatrixC &Vec);
 
  virtual const QMatrixC &CholeskyDecompose();
 
  virtual QMatrixC GetCholeskyDecomposition() const;
  //  QMatrixC GetCholeskyInverse() const;
 
  //  QMatrixC &CholeskyInvert();
 
  virtual void CholeskySolveX(const QVectorC &Y,QVectorC &X) const;
  QMatrixC GetSubMatrix(UInt_t nI,UInt_t nJ,UInt_t startI=0,UInt_t startJ=0);
 
  const QMatrixC GetSubMatrix(UInt_t nI,UInt_t nJ,UInt_t startI=0,
                              UInt_t startJ=0) const;
 
  QMatrixC GetSummedAreaMatrix() const;
 
 
  void Print(const char *opt="") const;
 
protected: 
  QMatrixC (const gsl_matrix_complex* mat);
  QMatrixC(const double* orig);
  
  QMatrixC(gsl_complex *data,const UInt_t nrow,const UInt_t ncol,
           const UInt_t tda=0);
  
  void SetMat() const;
  
  QComplex GetElement(const UInt_t i,const UInt_t j){
    if( i<GetNRow() && j<fNCol ) return fData+2*(i*MemRowCElements()+j);
    std::stringstream stream;
    stream<<"N row: "<<GetNRow()<<" Index: "<<i;
    stream<<"N col: "<<fNCol<<" Index: "<<j;
    DianaThrow( QError(QERR_OUT_OF_RANGE,__FILE__,__LINE__,stream.str()));
    return QComplex(Q_DOUBLE_DEFAULT);
  }
  
  const QComplex GetElement(const UInt_t i,const UInt_t j) const {
    if( i<GetNRow() && j<fNCol ) return fData+2*(i*MemRowCElements()+j);
    std::stringstream stream;
    stream<<"N row: "<<GetNRow()<<" Index: "<<i;
    stream<<"N col: "<<fNCol<<" Index: "<<j;
    DianaThrow( QError(QERR_OUT_OF_RANGE,__FILE__,__LINE__,stream.str()));
    return QComplex(Q_DOUBLE_DEFAULT);
  }
 
  UInt_t MemRowCElements() const {
    return (fNCol>fTda ? fNCol : fTda);
  }
 
  static void MatrixCopy(double *orig,UInt_t nrow,UInt_t ncol,double *dest,
                         UInt_t origTda=0,UInt_t destTda=0);
 
  //private:
protected:
  Bool_t fDataOwner; 
  UInt_t fSize; 
  UInt_t fNCol;
  UInt_t fTda; 
  UInt_t fAllocSize;
  double *fData;//[fSize]
 
  mutable gsl_matrix_complex *fMat;
  
  friend class QVectorC;
  friend class  QMatrixCArray;
  QObjectDef(QMatrixC,4);
};
 
// operations other than member functions
//QMatrix operator*(double t, const QMatrix&mat);
 
// formatted printout
 
Q_END_NAMESPACE
std::ostream& operator<<(std::ostream&s,const Diana::QMatrixC& m);
 
#endif