QMatrixC.cc

Go to the documentation of this file.

/*
* 
* Class QMatrixC
*
*/
 
#include "QMatrixC.hh"
#include "QVectorC.hh"
#include "QComplex.hh"
#include "QError.hh"
#include <TMath.h>
#include <cmath>
#include <TString.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
 
QObjectImp(Diana::QMatrixC);
 
#ifdef _HAVE_FFTW3_
#include <fftw3.h>
#endif
 
Q_BEGIN_NAMESPACE
 
void QMatrixC::MatrixCopy(double *orig,UInt_t nrow,UInt_t ncol,
                          double *dest,UInt_t origCTda,UInt_t destCTda)
{
  if(nrow == ncol && ncol == 0) return;
  if(orig && dest){
    
    if(origCTda<ncol) origCTda=ncol;
    if(destCTda<ncol) destCTda=origCTda;
    
    if(origCTda==ncol && destCTda==ncol)
      memcpy(dest,orig,2*nrow*ncol*sizeof(double));
    else
      for(UInt_t i=0;i<nrow;i++)
        memcpy(dest+2*i*destCTda,orig+2*i*origCTda,2*ncol*sizeof(double));
  }
  else if (!orig && !dest){
    // If neither of them are allocated, do nothing?
  }
  else
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Can not copy QMatrixC from origin to destination "
                      "before memory allocation!"));
 
}
 
QMatrixC::QMatrixC()
  : fDataOwner(true),
    fSize(0),
    fNCol(0),
    fTda(0),
    fAllocSize(0),
    fData(0),
    fMat(0)
{}
 
 
QMatrixC::QMatrixC(const UInt_t nrow,const UInt_t ncol)  
  : fDataOwner(true),
    fSize(0),
    fNCol(0),
    fTda(0),
    fAllocSize(0),
    fData(0),
    fMat(0)
{
  Resize(nrow,ncol);
}
 
QMatrixC::QMatrixC(const QVectorC& vec)
  : fDataOwner(true),
    fSize(0),
    fNCol(0),
    fTda(0),
    fAllocSize(0),
    fData(0),
    fMat(0)
{
  Resize(vec.Size(),1);
  if(vec.fStride==1)
    QVector::ArrayCopy(vec.fData, vec.fStride, fData, 1,fAllocSize);
  else{
    for(UInt_t i=0;i<vec.Size();i++)
      this->operator()(i,0)=vec[i];
  }
}
 
QMatrixC::QMatrixC(const QMatrixC& orig)
  : fDataOwner(true),
    fSize(0),
    fNCol(0),
    fTda(0),
    fAllocSize(0),
    fData(0),
    fMat(0)
{
  Resize(orig.GetNRow(),orig.GetNCol());
  MatrixCopy(orig.fData,GetNRow(),GetNCol(),fData,orig.MemRowCElements(),
             MemRowCElements());
}
 
QMatrixC::QMatrixC(const QMatrix& orig)
  : fDataOwner(true),
    fSize(0),
    fNCol(0),
    fTda(0),
    fAllocSize(0),
    fData(0),
    fMat(0)
{
  Resize(orig.GetNRow(),orig.GetNCol());
  for(UInt_t i = 0; i < orig.GetNRow(); i++) {
    for(UInt_t j = 0; j < orig.GetNCol(); j++) {
      this->operator()(i,j).SetRe(orig(i,j)); 
      this->operator()(i,j).SetIm(0); 
    }
  }
}
QMatrixC::QMatrixC(const gsl_matrix_complex* mat)
  : fDataOwner(true),
    fSize(0),
    fNCol(0),
    fTda(0),
    fAllocSize(0),
    fData(0),
    fMat(0)
{
  Resize(mat->size1,mat->size2);
  MatrixCopy(mat->data,GetNRow(),GetNCol(),fData,MemRowCElements(),mat->tda);
}
 
QMatrixC::QMatrixC(gsl_complex *data,const UInt_t nrow,const UInt_t ncol,
                   const UInt_t tda)
  : fDataOwner(false),
    fSize(2*ncol*nrow),
    fNCol(ncol),
    fTda(tda),
    fAllocSize(2*ncol*nrow),
    fData((double *)data),
    fMat(0)
{
  if(fTda<fNCol)
    fTda = fNCol;
  else
    fTda = tda;
}
QMatrixC::~QMatrixC()
{
  Clear();
}
void QMatrixC::Clear()
{
  if(fMat) free(fMat);
  if(fData && fDataOwner) QVector::ArrayFree(fData);
  fMat=NULL;
  fData=NULL;
  fSize=0;
  fNCol=0;
  fTda=0;
  fAllocSize=0;
}
 
void QMatrixC::Resize(const UInt_t nrow,const UInt_t ncol)
{
  if(GetNRow() == nrow && GetNCol() == ncol) return;
  if(fDataOwner){
    if(nrow*ncol > fAllocSize/2) {
      if(fData) QVector::ArrayFree(fData);
      fAllocSize=nrow*ncol*2;
      fData=QVector::ArrayAlloc(fAllocSize);
    }
    fTda = ncol;
    fSize = nrow*ncol*2;
    fNCol = ncol;
    if(fMat) free(fMat);
    fMat=NULL; 
  }
  else
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Attempt to resize a QMatrixC which does not own "
                      "the memory."));
}
void QMatrixC::Initialize(const QComplex& val)
{
  for(UInt_t i=0;i<GetNRow();i++) {
    for(UInt_t j=0;j<(UInt_t)(fNCol);j++) {
      operator()(i,j)=val;
    }
  }
}
void QMatrixC::SetToIdentity()
{
  SetMat();
  gsl_matrix_complex_set_identity(fMat);
}
void QMatrixC::SetCol(const UInt_t n,const QVectorC& vec)
{
 if(vec.Size()==GetNRow()){
   SetMat();
   vec.SetMathVector();
   gsl_matrix_complex_set_col(fMat,n,vec.fMathVec);
 }
 else {
     DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to set a column of a QMatrixC with a different size QVectorC"));
   }
}
 
void QMatrixC::SetRow(const UInt_t n,const QVectorC& vec)
{
  if(vec.Size()==GetNCol()){
    SetMat();
    vec.SetMathVector();
    gsl_matrix_complex_set_row(fMat,n,vec.fMathVec);
  }
  else {
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to set a row of a QMatrixC with a different size QVectorC"));
  }
}
QVectorC QMatrixC::GetCol(const UInt_t n) const
{
  QVectorC res(GetNRow());
  res.SetMathVector();
  SetMat();
  gsl_matrix_complex_get_col(res.fMathVec,fMat,n);
  return res;
}
 
QVectorC QMatrixC::GetRow(const UInt_t n) const
{
  QVectorC res(GetNCol());
  res.SetMathVector();
  SetMat();
  gsl_matrix_complex_get_row(res.fMathVec,fMat,n);
  return res;
}
QVectorC QMatrixC::GetDiagonal() const
{
  if(GetNCol() != GetNRow()) {
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"GetDiagonal can be called only on square matrices"));
  }
  QVectorC res(GetNCol());
  res.SetMathVector();
  SetMat();
  for(UInt_t i = 0; i < GetNCol(); i++) {
    res[i] = this->operator()(i,i);
  }
  return res;
}
QComplex QMatrixC::operator()(const UInt_t i,const UInt_t j)
{
  return GetElement(i,j);
}
 
const QComplex QMatrixC::operator()(const UInt_t i,const UInt_t j) const
{
  return GetElement(i,j);
}
 
const QMatrixC& QMatrixC::operator=(const QMatrixC& orig) 
{
  Resize(orig.GetNRow(),orig.fNCol);
  MatrixCopy(orig.fData,GetNRow(),GetNCol(),fData,
             orig.MemRowCElements(),MemRowCElements());
  if(fMat) free(fMat);
  fMat = NULL;
  return *this;
}
 
QMatrixC QMatrixC::operator-() const 
{
  return (*this)*(-1.);
}
 
const QMatrixC& QMatrixC::operator*=(const QComplex& t) 
{
  SetMat();
  gsl_matrix_complex_scale(fMat,t.Get_gsl_complex());
  return *this;
}
const QMatrixC& QMatrixC::operator*=(const double t) 
{
  QComplex tmp(t,0);
  return operator*=(tmp);
}
const QMatrixC& QMatrixC::operator*=(const QMatrixC& other) {  
  // multiply a matrix ths=this*other
  //gsl does not have matrix multiplications?!? 
  //  gsl_matrix* tmpM=gsl_matrix_alloc(other.GetNRow(),other.GetNCol()); 
 
  QMatrixC tmp(*this);
  Resize(GetNRow(),other.GetNCol());
  gsl_complex uno;
  uno.dat[0]=1;
  uno.dat[1]=0;
  gsl_complex zero;
  zero.dat[0]=0;
  zero.dat[1]=0;
  SetMat();
  other.SetMat();
  tmp.SetMat();
  gsl_blas_zgemm(CblasNoTrans,CblasNoTrans,uno,tmp.fMat,other.fMat,zero,fMat);
 
  return *this;
}
const QMatrixC& QMatrixC::operator/=(const QComplex& t) 
{
  return operator*=(1/t);
}
const QMatrixC& QMatrixC::operator/=(const double t) 
{
  return operator*=(1/t);
}
 
const QMatrixC& QMatrixC::operator+=(const QMatrixC& other) 
{
  SetMat();
  other.SetMat();
  gsl_matrix_complex_add(fMat,other.fMat);
  return *this;
}
 
const QMatrixC& QMatrixC::operator-=(const QMatrixC& other) 
{
  SetMat();
  other.SetMat();
  gsl_matrix_complex_sub(fMat,other.fMat);
  return *this;
}
 
const QMatrixC& QMatrixC::Mult(const QMatrixC& other) 
{
  SetMat();
  other.SetMat();
  gsl_matrix_complex_mul_elements(fMat,other.fMat);
  return *this;
}
 
const QMatrixC& QMatrixC::Div(const QMatrixC& other) 
{
  SetMat();
  other.SetMat();
  gsl_matrix_complex_div_elements(fMat,other.fMat);
  return *this;
}
QMatrixC QMatrixC::operator*(const QMatrixC& mat) const
{
  QMatrixC tmp(*this);
  tmp*=mat;
  return tmp;
}
QMatrixC QMatrixC::operator*(const double t) const
{ 
   QMatrixC tmp(*this);
   tmp*=t;
   return tmp;
}
QMatrixC QMatrixC::operator*(const QComplex& t) const
{
  QMatrixC tmp(*this);
  return tmp*=t;
}
QMatrixC QMatrixC::operator/(const QComplex& t) const
{
  QMatrixC tmp(*this);
  return tmp/=t;
}
QMatrixC QMatrixC::operator/(const double t) const
{
  QMatrixC tmp(*this);  
  tmp/=t;
  return tmp;
}
 
QMatrixC QMatrixC::operator+(const QMatrixC& mat) const
{
  QMatrixC tmp(*this);
  tmp+=mat;
  return tmp;
}
 
QMatrixC QMatrixC::operator-(const QMatrixC& mat) const
{
  QMatrixC tmp(mat);
  tmp-=mat;
  return tmp;
}
QVectorC QMatrixC::operator*(const QVectorC& vec) const
{
  QMatrixC v(vec);
  QMatrixC tmp((*this)*v);
  return tmp.GetCol(0);
}
QMatrixC QMatrixC::T() const 
{
  QMatrixC tmp(*this);
  tmp.Transpose();
  return tmp;
}
const QMatrixC &QMatrixC::Transpose()
{
  Resize(GetNCol(),GetNRow());
  SetMat();
  gsl_matrix_complex_transpose(fMat);
  return *this;
}
const QMatrixC &QMatrixC::H()
{
  Transpose();
  Conjugate();
  return *this;
}
QMatrixC QMatrixC::GetH() const {
  QMatrixC tmp(*this);
  tmp.H();
  return tmp;
}
const QMatrixC &QMatrixC::Conjugate()
{
  for(UInt_t i=0;i<GetNRow();i++){
    for(UInt_t j=0;j<GetNCol();j++){
      fData[2*(i*MemRowCElements()+j)+1]*=-1;
    }
  }
  return *this;
}
QMatrixC QMatrixC::GetConjugate() const
{
  QMatrixC tmp(*this);
  tmp.Conjugate();
  return tmp;
}
QMatrix QMatrixC::Magnitude() const
{
  SetMat();
  QMatrix mag(GetNRow(),GetNCol());
  for(UInt_t i = 0; i < GetNRow(); i++) {
    for(UInt_t j = 0; j < GetNCol(); j++) {
      mag(i,j) = this->operator()(i,j).GetMagnitude();
    }
  }
  return mag;
}
QMatrixC QMatrixC::Inv() const
{
  QMatrixC tmp(*this);
  tmp.Invert();
  return tmp;
}
const QMatrixC &QMatrixC::Invert()
{
  SetMat();
  int err=0;
  int sign=0; 
  gsl_permutation * p = gsl_permutation_calloc(GetNRow());
  gsl_matrix_complex *LU = gsl_matrix_complex_calloc(GetNRow(),fNCol);
  gsl_matrix_complex_memcpy(LU,fMat);
  int * signum = &sign;
  err=gsl_linalg_complex_LU_decomp(LU, p, signum);
  if(err)std::cout<<" ERROR: failed to decompose matrix "<<err<<std::endl;
  err=gsl_linalg_complex_LU_invert(LU, p, fMat);
  if(err)std::cout<<" ERROR: failed to decompose matrix "<<err<<std::endl;
  gsl_permutation_free(p);
  gsl_matrix_complex_free(LU);
  return *this;
}
QComplex QMatrixC::Det() const
{
  SetMat();
  int err=0;
  int sign=0;
  gsl_permutation * p = gsl_permutation_calloc(GetNRow());
  gsl_matrix_complex * tmp_ptr = gsl_matrix_complex_calloc(GetNRow(),fNCol);
  int * signum = &sign;
  gsl_matrix_complex_memcpy(tmp_ptr, fMat);
  err=gsl_linalg_complex_LU_decomp(tmp_ptr, p, signum);
  if(err)std::cout<<" ERROR: failed to compute matrix determinant "<<err<<std::endl;
  gsl_complex res = gsl_linalg_complex_LU_det(tmp_ptr, *signum);
  gsl_permutation_free(p);
  gsl_matrix_complex_free(tmp_ptr);
  return QComplex(&res);
}
QVector QMatrixC::HermitianEigenValues()  
{
  if(GetNRow()!=GetNCol())
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Can not calculate eigen values of non square matrix"
                      " using HermitianEigenValues()."));
  QVector Ret(GetNCol());
  QVectorC diag=GetDiagonal();
  SetMat();
  Ret.SetMathVector();
  gsl_eigen_herm_workspace *wSpace = gsl_eigen_herm_alloc(GetNCol());
  int err = gsl_eigen_herm(fMat,Ret.fMathVec,wSpace);
  if(err) 
    std::cout << " ERROR: Could not calculate eigenvalues of hermitian "
      "matrix. " << err << std::endl;
  for(UInt_t i=0;i<GetNRow();i++){
    for(UInt_t j=0;j<i;j++){
      (*this)(i,j)=(*this)(j,i).Conj();
    }
    (*this)(i,i)=diag[i];
  }
  gsl_eigen_herm_free(wSpace);
  Ret.Sort(false);
  return Ret;
}
void QMatrixC::HermitianEigenSystem(QVector &Val,QMatrixC &Vec)
{
  if(GetNRow()!=GetNCol())
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Can not calculate eigensystem of non square matrix"
                      " using HermitianEigenSystem()."));
  Val.Resize(GetNRow());
  Vec.Resize(GetNRow(),GetNCol());
  QVectorC diag=GetDiagonal();
  SetMat();
  Val.SetMathVector();
  Vec.SetMat();
  gsl_eigen_hermv_workspace *wSpace = gsl_eigen_hermv_alloc(GetNCol());
  int err = gsl_eigen_hermv(fMat,Val.fMathVec,Vec.fMat,wSpace);
  if(err) 
    std::cout << " ERROR: Could not calculate eigenvalues of hermitian "
      "matrix. " << err << std::endl;
  for(UInt_t i=0;i<GetNRow();i++){
    for(UInt_t j=0;j<i;j++){
      (*this)(i,j)=(*this)(j,i).Conj();
    }
    (*this)(i,i)=diag[i];
  }
  gsl_eigen_hermv_sort(Val.fMathVec,Vec.fMat,GSL_EIGEN_SORT_ABS_DESC);
  gsl_eigen_hermv_free(wSpace);
}
const QMatrixC &QMatrixC::CholeskyDecompose() {
  if(GetNRow()!=GetNCol())
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Can not Cholsky decompose non-square matrix!"));
  SetMat();
  int err=gsl_linalg_complex_cholesky_decomp(fMat);
  if(err) std::cout << " ERROR: Failed to Cholesky decompose matrix! " << err
                    << std::endl;
  return *this;
}
QMatrixC QMatrixC::GetCholeskyDecomposition() const {
  QMatrixC out(*this);
  out.CholeskyDecompose();
  return out;
}
void QMatrixC::CholeskySolveX(const QVectorC &Y,QVectorC &X) const {
  if(Y.Size()!=GetNRow())
    throw QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                 "Vector and matrix dimensions do not match "
                 "in CholeskySolveX");
  X.Resize(Y.Size());
  SetMat();
  X.SetMathVector();
  Y.SetMathVector();
  int err = gsl_linalg_complex_cholesky_solve(fMat,Y.fMathVec,X.fMathVec);
  if(err) std::cout << " ERROR: Failed to cholesky solve for X!" << err 
                    << std::endl;
}
void QMatrixC::SetMat() const
{
  if(!fMat) {        
    fMat = (gsl_matrix_complex*) malloc(sizeof(gsl_matrix_complex));
    
    fMat->size1=GetNRow();
    fMat->size2=fNCol;
    fMat->owner=1;
    fMat->data=(double *)fData;
    fMat->block=0;
    fMat->tda=(fTda>fNCol ? fTda : fNCol);
  }
}
 
QMatrixC QMatrixC::GetSubMatrix(UInt_t nI,UInt_t nJ,UInt_t startI,
                                UInt_t startJ)
{
  if(startI + nI > GetNRow() || startJ + nJ > GetNCol())
    throw(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                 Form("Cannot extract a submatrix of size %ix%i starting from"
                      " element %i,%i from a matrix of size %ix%i",
                      nI,nJ,startI,startJ,GetNRow(),GetNCol())));
  gsl_complex *firstElem =
    (gsl_complex *)fData+startI*MemRowCElements()+startJ;
  QMatrixC ret(firstElem,nI,nJ,fTda);
  return ret;
}
const QMatrixC QMatrixC::GetSubMatrix(UInt_t nI,UInt_t nJ,UInt_t startI,
                                      UInt_t startJ) const
{
  if(startI + nI > GetNRow() || startJ + nJ > GetNCol())
    throw(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                 Form("Cannot extract a submatrix of size %ix%i starting from"
                      " element %i,%i from a matrix of size %ix%i",
                      nI,nJ,startI,startJ,GetNRow(),GetNCol())));
  gsl_complex *firstElem =
    (gsl_complex *)fData+startI*MemRowCElements()+startJ;
  QMatrixC ret(firstElem,nI,nJ,fTda);
  return ret;
}
 
QMatrixC QMatrixC::GetSummedAreaMatrix() const {
  QMatrixC Ret(*this);
  for(size_t i=0;i<GetNRow();i++) {
    for(size_t j=0;j<GetNCol();j++) {
      if(i>0) Ret(i,j) += Ret(i-1,j);
      if(j>0) Ret(i,j) += Ret(i,j-1);
      if(i>0 && j>0) Ret(i,j) -= Ret(i-1,j-1);
    }
  }
  return Ret;
}
 
void QMatrixC::Print(const char *opt) const {
  std::cout << "fSize " << fSize << std::endl;
  std::cout << "fNRow " << GetNRow() << std::endl;
  std::cout << "fNCol " << fNCol << std::endl;
  std::cout << "fTda "  <<  fTda << std::endl;
  std::cout << "fData " << fData << std::endl;
  std::cout << "fData    ";
  for(UInt_t i=0;i<TMath::Min(GetNRow(),(UInt_t)3);i++){
    for(UInt_t j=0;j<TMath::Min(GetNCol(),(UInt_t)3);j++){
      std::cout << (*this)(i,j) << "\t";
    }
    std::cout << std::endl;
  }
  std::cout << std::endl;
}
 
Q_END_NAMESPACE
 
/*************************************
 *      NON MEMBER FUNCTIONS
 * **********************************/
 
  Diana::QMatrixC operator*(double t, const Diana::QMatrixC&mat){
  Diana::QMatrixC tmp(mat);
  tmp*=t;
  return tmp;
}
 
std::ostream& operator<<(std::ostream&s,const Diana::QMatrixC&mat){
  for (UInt_t i=0;i<mat.GetNRow();i++){
    for (UInt_t j=0;j<mat.GetNCol();j++) s<<mat(i,j)<<" ";
    s<<std::endl;
  }
  return s;
}