QMatrix.cc

Go to the documentation of this file.

/*
* 
* Class QMatrix
*
*/
 
#include "QMatrix.hh"
#include "QVector.hh"
#include "QError.hh"
#include <cmath>
#include <TString.h>
#include <sstream>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
 
using std::cout;
using std::endl;
using std::map;
using std::ostream;
 
#ifdef _HAVE_FFTW3_
#include <fftw3.h>
#endif
 
QObjectImp(Diana::QMatrix);
 
Q_BEGIN_NAMESPACE
void QMatrix::MatrixCopy(const double *orig,const UInt_t Nrow,
                         const UInt_t Ncol,double *dest,
                         UInt_t origTda,UInt_t destTda){
  if(orig && dest){
 
    if(origTda<Ncol) origTda=Ncol;
    if(destTda<Ncol) destTda=origTda;
 
    if(origTda==Ncol && destTda==Ncol)
      memcpy(dest,orig,Nrow*Ncol*sizeof(double));
    else
      for(UInt_t i=0;i<Nrow;i++)
        memcpy(dest+i*destTda,orig+i*origTda,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 QMatrix from origin to destination "
                      "before memory allocation!"));
}
 
QMatrix::QMatrix()
  : fDataOwner(true),
    fSize(0),
    fNCol(0),
    fTda(0),
    fAllocSize(0),
    fData(0),
    fMat(0)
{}
 
 
QMatrix::QMatrix(UInt_t nrow,UInt_t ncol)  
  : fDataOwner(true),
    fSize(0),
    fNCol(0),
    fTda(0),
    fAllocSize(0),
    fData(0),
    fMat(0)
{
  Resize(nrow,ncol);
}
 
QMatrix::QMatrix(const QVector& vec)
  : fDataOwner(true),
    fSize(0),
    fNCol(0),
    fTda(0),
    fAllocSize(0),
    fData(0),
    fMat(0)
{
  Resize(vec.Size(),1);
  QVector::ArrayCopy(vec.fData, vec.fStride, fData, 1, fSize);
}
 
QMatrix::QMatrix(const QMatrix& 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.MemRowElements(),
             MemRowElements());
}
 
QMatrix::QMatrix(const gsl_matrix* 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,mat->tda,MemRowElements());
}
QMatrix::QMatrix(double *data,const UInt_t nrow,const UInt_t ncol,
                 const UInt_t tda)
  : fDataOwner(false),
    fSize(nrow*ncol),
    fNCol(ncol),
    fTda(tda),
    fAllocSize(nrow*ncol),
    fData(data),
    fMat(0)
{
  if(fTda<fNCol)
    fTda = fNCol;
  else
    fTda = tda;
  fMat=NULL;
}
QMatrix::~QMatrix(){
  if(fMat) free(fMat);
  if(fData && fDataOwner) QVector::ArrayFree(fData);
  fMat=NULL;
  fData=NULL;
  fSize=0;
  fNCol=0;
  fTda=0;
  fAllocSize=0;
}
void QMatrix::Clear()
{
  Resize(0,0);
}
void QMatrix::Resize(UInt_t nrow,UInt_t ncol)
{
  if(GetNRow() == nrow && GetNCol() == ncol) return;
  if(fDataOwner){
    if(nrow*ncol > fAllocSize) {
      if(fData) QVector::ArrayFree(fData);
      fAllocSize=nrow*ncol;
      fData=QVector::ArrayAlloc(fAllocSize);
    }
    fTda = ncol;
    fSize = nrow*ncol;
    fNCol = ncol;
    if(fMat) free(fMat);
    fMat = NULL;
  }
  else
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Attempt to resize a QMatrix which does not own "
                      "the memory."));
}
void QMatrix::Initialize(double val)
{
  // We must take care here not to overwrite the data in the trailing
  // dimensions.
  for(UInt_t i=0;i<(UInt_t) GetNRow(); i++)
    for(UInt_t j=0;j<(UInt_t) GetNCol(); j++)
      fData[MemRowElements()*i+j]=val;
 
}
 
void QMatrix::SetToIdentity()
{
  SetMat();
  gsl_matrix_set_identity(fMat);
}
 
void QMatrix::SetCol(UInt_t n,const QVector& vec)
{
  if(vec.Size()==GetNRow()) {
    SetMat();
    vec.SetMathVector();
    gsl_matrix_set_col(fMat,n,vec.fMathVec);
  }
  else
    {
      DianaThrow (QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to set a column of a QMatrix with a different size QVector"));
    }
}
 
void QMatrix::SetRow(UInt_t n,const QVector& vec)
{
  if(vec.Size()==GetNCol()) {
    SetMat();
    //  assert(n<=nrow()&&vec.size()==ncol());
    vec.SetMathVector();
    gsl_matrix_set_row(fMat,n,vec.fMathVec);
  }
  else {
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to set a row of a QMatrix with a different size QVector"));
  }
}
 
void QMatrix::Set(const map<int,int> & p){
 
  map<int,int>::const_iterator iter;
  int n_ele = p.size();
 
  this->Resize(n_ele,2);
 
  int i=0;
 
  SetMat();
  for(iter = p.begin(); iter!= p.end(); iter++) {
    gsl_matrix_set(fMat,i,0,iter->first);
    gsl_matrix_set(fMat,i,1,iter->second);    
    i++;
  }
}
QVector QMatrix::GetCol(UInt_t n)
{
  return QVector(fData+n,GetNRow(),MemRowElements());
}
QVector QMatrix::GetCol(UInt_t n) const
{
  QVector res(GetNRow());
  res.SetMathVector();
  SetMat();
  gsl_matrix_get_col(res.fMathVec,fMat,n);
  return res;
}
 
QVector QMatrix::GetRow(UInt_t n)
{
  return QVector(fData+n*MemRowElements(),fNCol,1);
}
QVector QMatrix::GetRow(UInt_t n) const
{
  QVector res(GetNCol());
  res.SetMathVector();
  SetMat();
  gsl_matrix_get_row(res.fMathVec,fMat,n);
  return res;
}
QVector QMatrix::GetDiaganol()
{
  return QVector(fData,GetNRow(),MemRowElements()+1);
}
QVector QMatrix::GetDiaganol() const 
{
  QVector res(GetNCol());
  for(UInt_t i=0;i<GetNCol();i++)
    res[i]=(*this)(i,i);
  return res;
}
//e che e' sto coso?!??
QVector QMatrix::GetRowByColumnIntValue(UInt_t col,int val) const
{
  for(UInt_t i = 0; i < GetNRow(); i++) {
    if(lround(operator()(i,col)) == val) 
      return GetRow(i);
  }
  QVector dummy;
  return dummy; 
}
 
double &QMatrix::GetElement(UInt_t i,UInt_t j)
{
  SetMat();
  if( !(i<GetNRow() && j<GetNCol()) ) {
    std::stringstream stream;
    stream<<"N row: "<<GetNRow()<<" Index: "<<i;  
    stream<<"N col: "<<GetNCol()<<" Index: "<<j;
    DianaThrow( QError(QERR_OUT_OF_RANGE,__FILE__,__LINE__,stream.str()));
  }
  return fData[i*MemRowElements()+j];
}
 
const double &QMatrix::GetElement(UInt_t i, UInt_t j) const
{
  SetMat();
  if( !(i<GetNRow() && j<GetNCol()) ) {
    std::stringstream stream;
    stream<<"N row: "<<GetNRow()<<" Index: "<<i;  
    stream<<"N col: "<<GetNCol()<<" Index: "<<j;
    DianaThrow( QError(QERR_OUT_OF_RANGE,__FILE__,__LINE__,stream.str()));
  }
  return fData[i*MemRowElements()+j];
}
 
 
const QMatrix& QMatrix::operator=(const QMatrix& orig) 
{
  Resize(orig.GetNRow(),orig.GetNCol());
  if(fSize)
    QMatrix::MatrixCopy(orig.fData,GetNRow(),GetNCol(),fData,
                        orig.MemRowElements(),MemRowElements());
  if(fMat) free(fMat);
  fMat=NULL;
  return *this;
}
 
QMatrix QMatrix::operator-() 
{
  return operator*=(-1.);
}
 
const QMatrix& QMatrix::operator*=(double t) 
{
  SetMat();
  gsl_matrix_scale(fMat,t);
  return *this;
}
 
const QMatrix& QMatrix::operator*=(const QMatrix& other) 
{  
  // multiply a matrix ths=this*other
  //gsl does not have matrix multiplications?!? 
  //  gsl_matrix* tmpM=gsl_matrix_alloc(other.GetNRow(),other.GetNCol()); 
  QMatrix tmp(*this);
  Resize(GetNRow(),other.GetNCol());
  SetMat();
  other.SetMat();
  tmp.SetMat();
  gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,tmp.fMat,other.fMat,0.0,fMat);
  return *this;
}
 
const QMatrix& QMatrix::operator/=(double t) 
{
  return operator*=(1/t);
}
 
const QMatrix& QMatrix::operator+=(const QMatrix& other) 
{
  SetMat();
  other.SetMat();
  gsl_matrix_add(fMat,other.fMat);
  return *this;
}
 
const QMatrix& QMatrix::operator-=(const QMatrix& other) 
{
  SetMat();
  other.SetMat();
  gsl_matrix_sub(fMat,other.fMat);
  return *this;
}
 
const QMatrix& QMatrix::Mult(const QMatrix& other) 
{
  SetMat();
  other.SetMat();
  gsl_matrix_mul_elements(fMat,other.fMat);
  return *this;
}
 
const QMatrix& QMatrix::Div(const QMatrix& other) 
{
  SetMat();
  other.SetMat();
  gsl_matrix_div_elements(fMat,other.fMat);
  return *this;
}
 
QMatrix QMatrix::operator*(const QMatrix& mat) const
{
  QMatrix tmp(*this);
  tmp*=mat;
  return tmp;
}
 
QMatrix QMatrix::operator*(double t) const
{ 
   QMatrix tmp(*this);
   tmp*=t;
   return tmp;
}
 
QMatrix QMatrix::operator/(double t) const
{
  QMatrix tmp(*this);  
  tmp/=t;
  return tmp;
}
 
QMatrix QMatrix::operator+(const QMatrix& mat) const
{
  QMatrix tmp(*this);
  tmp+=mat;
  return tmp;
}
 
QMatrix QMatrix::operator-(const QMatrix& mat) const
{
  QMatrix tmp(mat);
  tmp-=mat;
  return tmp;
}
 
QVector QMatrix::operator*(const QVector& vec) const
{
  QVector v; // Allocate space for the response, because the matrix m will disappear
  QMatrix m(vec);
  v = ((*this)*m).GetCol(0);
  return v;
}
 
QMatrix QMatrix::T() const 
{
  QMatrix tmp(*this);
  tmp.Transpose();
  return tmp;
}
 
const QMatrix &QMatrix::Transpose()
{
  if(!fDataOwner && GetNRow()!=GetNCol())
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
              "Can not transpose matrix if not data owner or square!"));
  SetMat();
  gsl_matrix_transpose(fMat);
  if(GetNRow()!=GetNCol())
    fTda=fNCol=GetNRow();
  return *this;
}
 
QMatrix QMatrix::Inv() const
{
  QMatrix tmp(*this);
  tmp.Invert();
  return tmp;
}
 
const QMatrix &QMatrix::Invert()
{
  SetMat();
  int err=0;
  int sign=0; 
  gsl_permutation * p = gsl_permutation_calloc(GetNRow());
  gsl_matrix * LU = gsl_matrix_calloc(GetNRow(), fNCol);
  int * signum = &sign;
  gsl_matrix_memcpy(LU, fMat);
  err=gsl_linalg_LU_decomp(LU, p, signum);
  if(err)cout<<" ERROR: failed to decompose matrix "<<err<<endl;
  err=gsl_linalg_LU_invert(LU, p, fMat);
  if(err)std::cout<<" ERROR: failed to invert LU matrix "<<err<<std::endl;
  gsl_permutation_free(p);
  gsl_matrix_free(LU);
  return *this;
}
 
double QMatrix::Det() const
{
  SetMat();
  int err=0;
  int sign=0;
  gsl_permutation * p = gsl_permutation_calloc(GetNRow());
  gsl_matrix * tmp_ptr = gsl_matrix_calloc(GetNRow(),fNCol);
  int * signum = &sign;
  gsl_matrix_memcpy(tmp_ptr, fMat);
  err=gsl_linalg_LU_decomp(tmp_ptr, p, signum);
  if(err)cout<<" ERROR: failed to compute matrix determinant "<<err<<endl;
  double res = gsl_linalg_LU_det(tmp_ptr, *signum);
  gsl_permutation_free(p);
  gsl_matrix_free(tmp_ptr);
  return res;
}
const QMatrix &QMatrix::CholeskyDecompose() {
  SetMat();
  int err=gsl_linalg_cholesky_decomp(fMat);
  if(err) std::cout << " ERROR: Failed to Cholesky decompose the Matrix " 
                    << err << std::endl;
  return *this;
}
QMatrix QMatrix::GetCholeskyDecomposition() const {
  QMatrix out(*this);
  out.CholeskyDecompose();
  return out;
}
void QMatrix::CholeskySolveX(const QVector &Y,QVector &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_cholesky_solve(fMat,Y.fMathVec,X.fMathVec);
  if(err) std::cout << " ERROR: Failed to cholesky solve for X!" << err
                    << std::endl;
}
QMatrix QMatrix::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())));
  double *firstElem = fData+startI*MemRowElements()+startJ;
  QMatrix ret(firstElem,nI,nJ,fTda);
  return ret;
}
const QMatrix QMatrix::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())));
  double *firstElem = fData+startI*MemRowElements()+startJ;
  QMatrix ret(firstElem,nI,nJ,fTda);
  return ret;
}
 
QMatrix QMatrix::GetSummedAreaMatrix() const {
  QMatrix 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 QMatrix::SetMat() const
{
  if(!fMat){        
    fMat = (gsl_matrix*) malloc(sizeof(gsl_matrix));
    fMat->size1=GetNRow();
    fMat->size2=fNCol;
    fMat->owner=1;
    fMat->data=fData;
    fMat->block=0;
    fMat->tda=fTda;//CMT CHECK
  }
}
 
Q_END_NAMESPACE
 
/*************************************
 *      NON MEMBER FUNCTIONS
 * **********************************/
 
Diana::QMatrix operator*(double t, const Diana::QMatrix&mat){
  Diana::QMatrix tmp(mat);
  tmp*=t;
  return tmp;
}
 
ostream& operator<<(ostream&s,const Diana::QMatrix&mat){
  for (UInt_t i=0;i<mat.GetNRow();i++){
    for (UInt_t j=0;j<mat.GetNCol();j++) s<<mat(i,j)<<" ";
    s<<endl;
  }
  return s;
}