QMatrixCArray.cc

Go to the documentation of this file.

#include <TString.h>
#include "QMatrixCArray.hh"
 
ClassImp(Diana::QMatrixCArray)
 
Q_BEGIN_NAMESPACE
 
QMatrixCArray::QMatrixCArray()
: fDataOwner(true),
  fSize(0),
  fNMatrix(0),
  fNRow(0),
  fNCol(0),
  fData(0)
{}
QMatrixCArray::QMatrixCArray(const QMatrixCArray &in)
  : fDataOwner(true),
    fSize(0),
    fNMatrix(0),
    fNRow(0),
    fNCol(0),
    fData(0)
{
  Resize(in.fNMatrix,in.fNRow,in.fNCol);
  memcpy(fData,in.fData,fSize*sizeof(double));
}
QMatrixCArray::QMatrixCArray(UInt_t nMatrix,UInt_t nRow,UInt_t nCol)
  : fDataOwner(true),
    fSize(0),
    fNMatrix(0),
    fNRow(0),
    fNCol(0),
    fData(0)
{
  Resize(nMatrix,nRow,nCol);
}
    
QMatrixCArray::~QMatrixCArray() { 
  Clear(); 
}
void QMatrixCArray::Clear() {
  if(fDataOwner && fData) 
    QVector::ArrayFree(fData);
  fNMatrix=0;
  fSize=0;
  fNRow=0;
  fNCol=0;
  fData=0;
}
const QMatrixCArray &QMatrixCArray::operator=(const QMatrixCArray &in){
  Clear();
  Resize(in.fNMatrix,in.fNRow,in.fNCol);
  memcpy(fData,in.fData,fSize*sizeof(double));
  return *this;
}
void QMatrixCArray::Resize(UInt_t nMatrix,UInt_t nRow,UInt_t nCol){
  if(nMatrix == GetNMatrix() && nRow == GetNRow() && nCol == GetNCol())
    return;
  if(!fDataOwner)
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Attemp to resize QMatrixCArray when not the "
                      "data owner"));
  if(fData) QVector::ArrayFree(fData);
  fSize=nMatrix*nRow*nCol*2;
  fData=QVector::ArrayAlloc(fSize);
  fNMatrix=nMatrix;
  fNRow=nRow;
  fNCol=nCol;
}
void QMatrixCArray::Initialize(double Re,double Im) {
  for(UInt_t i=0;i<fSize;i+=2){
    fData[i]=Re;fData[i+1]=Im;
  }
}
void QMatrixCArray::Initialize(const QComplex &Val) {
  Initialize(Val.Re(),Val.Im());
}
QMatrixC QMatrixCArray::GetMatrix(UInt_t idx){
  if(idx>fNMatrix)
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      Form("Matrix index out of bounds (%ui>%ui)",
                           idx,fNMatrix)));
  QMatrixC A((gsl_complex *)fData+idx*fNCol,fNRow,fNCol,fNMatrix*fNCol);
  return A;
}
const QMatrixC QMatrixCArray::GetMatrix(UInt_t idx) const {
  if(idx>fNMatrix)
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      Form("Matrix index out of bounds (%ui>%ui)",
                           idx,fNMatrix)));
  QMatrixC A((gsl_complex *)fData+idx*fNCol,fNRow,fNCol,fNMatrix*fNCol);
  return A;
}
Diana::QVectorC QMatrixCArray::GetVector(UInt_t idx1,UInt_t idx2){
  if(idx1>fNRow || idx2>fNCol)
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      Form("Vector indices out of bounds (%ui,%ui)>"
                           "(%ui,%ui)",idx1,idx2,fNRow,fNCol)));
  QVectorC A((gsl_complex *)fData+idx1*fNCol*fNMatrix+idx2,fNMatrix,fNCol);
  return A;
}
const Diana::QVectorC QMatrixCArray::GetVector(UInt_t idx1,UInt_t idx2) const{
  if(idx1>fNRow || idx2>fNCol)
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      Form("Vector indices out of bounds (%ui,%ui)>"
                           "(%ui,%ui)",idx1,idx2,fNRow,fNCol)));
  QVectorC A((gsl_complex *)fData+idx1*fNCol*fNMatrix+idx2,fNMatrix,fNCol);
  return A;
}
void QMatrixCArray::Invert() {
  for(UInt_t i=0;i<GetNMatrix();i++){
    QMatrixC C=GetMatrix(i);
    C.Invert();
  }
}
void QMatrixCArray::CholeskyDecompose() {
  for(UInt_t i=1;i<GetNMatrix();i++){
    QMatrixC C=GetMatrix(i);
    std::cout<< "CholeskyDecomp Matrix " << i << std::endl;
    C.CholeskyDecompose();
  }
}
void QMatrixCArray
::CholeskySolveX(const Diana::QVectorC &Y,Diana::QVectorC &X) const {
  if(Y.Size()!=GetNMatrix()*GetNRow())
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "CholeskySolveX: Dimensions do not match!"));
  X.Resize(Y.Size());
  for(UInt_t i=0;i<GetNMatrix();i++){
    Diana::QVectorC subX=X.GetSubVector(GetNRow(),i*GetNRow());
    Diana::QVectorC subY=Y.GetSubVector(GetNRow(),i*GetNRow());
    Diana::QMatrixC subM=GetMatrix(i);
    subM.CholeskySolveX(Y,X);
  }
}
Q_END_NAMESPACE