TDSPMatrix.h

Go to the documentation of this file.

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef TDSPMATRIX_H
#define TDSPMATRIX_H

#include <TH2.h>
#include <TArrayI.h>
#include <TDSPVector.h>



#define LINALG_INDEX_OPT 1
#define LINALG_MAT_ROWWISE 1


//
enum eMatrixApplyMode { 
  kMatrixApplyReal,      
  kMatrixApplyConjA,     
  kMatrixApplyConjB,     
  kMatrixApplyConjAB     
};

class TDSPMatrix : public TDSPVector  {

  static TH2   *fHistogram2;
  char*         fZTitle; 
 protected:
  
  UInt_t fRows;   
  UInt_t fCols;   
  
 public: 
  
  friend istream    &operator>>(istream&,TDSPMatrix&);
  friend ostream    &operator<<(ostream&,TDSPMatrix&);
  friend istream    &operator>>(istream&,TDSPMatrix*);
  friend ostream    &operator<<(ostream&,TDSPMatrix*);

  TDSPMatrix();
  TDSPMatrix(Int_t rows, Int_t cols);

  ~TDSPMatrix();
  
  TComplex&   Swap(Int_t, Int_t, Int_t, Int_t);  
  void        SetSize(Int_t rows, Int_t cols);
  UInt_t      GetRows() const { return fRows; }; 
  UInt_t      GetCols() const { return fCols; }; 
  void        SetRows(UInt_t, Bool_t savedata=kFALSE); 
  void        SetCols(UInt_t, Bool_t savedata=kFALSE); 
  virtual void SetLen(Int_t num, Bool_t savedata=kFALSE); 

  
  Int_t       Clo() const { return 0;};
  Int_t       Chi() const { return fCols-1;};
  
  Int_t       Rlo() const { return 0;};
  Int_t       Rhi() const { return fRows-1;};  

  TDSPMatrix* Dup() { return Copy(NULL);};                
  TDSPMatrix* Copy(TDSPMatrix* n);                   
  TDSPVector* Apply(TDSPVector*, TDSPVector *out=NULL, eMatrixApplyMode m=kMatrixApplyReal);
  TDSPVector* RightApply(TDSPVector*, TDSPVector* out =NULL); 
  TDSPMatrix* RightApply(TDSPMatrix*, TDSPMatrix* out =NULL); 
  TDSPVector* LeftApply( TDSPVector*, TDSPVector* out =NULL); 
  TDSPMatrix* LeftApply( TDSPMatrix*, TDSPMatrix* out =NULL); 

  TDSPVector* AdjointRightApply(TDSPVector*, TDSPVector* out =NULL); 
  TDSPMatrix* AdjointRightApply(TDSPMatrix*, TDSPMatrix* out =NULL); 
  TDSPVector* AdjointLeftApply( TDSPVector*, TDSPVector* out =NULL); 
  TDSPMatrix* AdjointLeftApply( TDSPMatrix*, TDSPMatrix* out =NULL); 
  TDSPVector* TransposedRightApply(TDSPVector*, TDSPVector* out =NULL); 
  TDSPVector* TransposedLeftApply( TDSPVector*, TDSPVector* out =NULL); 

  TDSPMatrix* Transform(TDSPMatrix* H,TDSPMatrix* result=NULL); 
  TDSPMatrix* AdjointTransform(TDSPMatrix* H,TDSPMatrix* result=NULL); 

  TDSPMatrix* Transposed(TDSPMatrix *r=NULL); 
  void        TransposeMe(); 
  TDSPMatrix* Adjoint(TDSPMatrix *r=NULL); 
  TComplex    Trace(); 
  Double_t    Trace2(); 

  Bool_t      DiagZero(); 
  Bool_t      DiagZeros(); 
  Bool_t      IsDiag(); 
  Bool_t      IsSquare() const { return fRows==fCols;}; 
  Bool_t      IsHermitian(); 

  TDSPMatrix* ReversedColumns(TDSPMatrix *r=NULL); 
  TDSPMatrix* ReversedRows(TDSPMatrix *r=NULL); 
  void        ReverseColumns(); 
  void        ReverseRows(); 

  TDSPVector* GetDiag(TDSPVector *a=NULL); 
  void        SetDiag(TDSPVector *a); 
  void        SetDiag(TComplex a); 
  void        SetNonDiag(TComplex a); 
  void        AddDiag(TComplex a); 
  void        AddDiag(TDSPVector* a); 
  void        MultDiag(TComplex a); 
  void        MultDiag(TDSPVector* a); 

  Double_t    Redundancy() const { return ((Double_t)fRows)/((Double_t)fCols); };
  
  TDSPMatrix *SubMatrix(Int_t atRow, Int_t atCol, 
            Int_t m, Int_t n, TDSPMatrix*out=NULL); 
  TDSPMatrix *TransposedSubMatrix(Int_t atRow, Int_t atCol, 
                  Int_t m, Int_t n, TDSPMatrix*out=NULL); 
  TComplex   *MoveVec(Int_t m, Int_t n);       
  TDSPVector* GetColumn(Int_t c=0, TDSPVector *out=NULL); 
  TDSPVector* GetRow(Int_t r=0, TDSPVector *out=NULL);    
  void        SetColumn(UInt_t c, TDSPVector *in, UInt_t ifrom=0, UInt_t ito=0); 
  void        SetRow(UInt_t r, TDSPVector *in, UInt_t ifrom=0, UInt_t ito=0);  
  void        RowConvolute(TComplex* inp, TComplex *out,
               Int_t N=0,
               Bool_t cyclic=kTRUE ); 

  TDSPVector* RowConvolute(TDSPVector*inp, TDSPVector*out=NULL, 
               Bool_t cyclic=kTRUE,
               Bool_t zerosbefore=kTRUE); 

  TH2*        Draw(Option_t *option = "hist surf4", 
           Double_t dx=1, Double_t xoff=0,
           Double_t dy=1, Double_t yoff=0,
           TH2 *h = NULL);
  void        SetZTitle(char *name); 
  void        FFTShiftMe(Int_t dim); 

  //
  //
  

#ifdef HAVE_MATLAB
  TDSPMatrix &operator=(TmxArray &m);
#endif
  
  
  TDSPMatrix  &operator=(TDSPMatrix& x) { x.Copy(this);return *this;};
  void        Print();
  void        Input();
  void        Unit(TComplex diag=1.0) { Delta(diag);}; 
  void        Delta(TComplex diag=1.0);   
  void        Fourier(); 
  void        Hilbert(); 
  
  TComplex&   Element(Int_t i, Int_t j) { return fVec[i*fCols+j];};
  TComplex&   operator()(UInt_t row, UInt_t col);

  
  // Routines based on the LU-Decomposition
  // (see LUComp.cpp)
  
  void        LUDecompose(TArrayI& Index, Int_t &d);
  void        LUBacksubst(TArrayI& Index, TDSPVector& B);
  void        LUBacksubst(TArrayI& Index, TDSPMatrix& B);
  void        LUSolveLinear(TDSPVector& B); 
  void        LUSolveLinear(TDSPMatrix& B); 
  
  TDSPMatrix* Inverse(TDSPMatrix *result=NULL);  
  TComplex    Det(); 

  // other external routines
  //(see SVD.cpp)

  TDSPMatrix* PseudoInverse(TDSPMatrix *result=NULL);
  TDSPVector* SVD(TDSPVector *S=NULL,
          TDSPMatrix *U=NULL, 
          TDSPMatrix *VStar=NULL); 
  ClassDef(TDSPMatrix,1)

};

#include <MatrixRoutines.h>



inline TDSPVector* TDSPMatrix::Apply(TDSPVector *inp, TDSPVector *out, eMatrixApplyMode m) {
  switch(m) {
  case kMatrixApplyReal :
    return RightApply(inp,out);
  case kMatrixApplyConjA :
    return AdjointRightApply(inp,out);
  case kMatrixApplyConjB :
    Error("Apply","error");
    return NULL;
  case  kMatrixApplyConjAB :
    Error("Apply","error");
    return NULL;
  }
  return NULL;
}


inline TDSPMatrix* TDSPMatrix::Transform(TDSPMatrix* H, TDSPMatrix *result) {
  
  TDSPMatrix* b = LeftApply(H);
  result=H->AdjointLeftApply(b,result);
  delete b;
  
  return result;

}


inline TDSPMatrix* TDSPMatrix::AdjointTransform(TDSPMatrix* H, TDSPMatrix *result) {
  
  TDSPMatrix* b = H->AdjointRightApply(this);
  result = H->LeftApply(b,result);
  delete b;
  return result;

}

inline TDSPMatrix* TDSPMatrix::LeftApply( TDSPMatrix *inp, TDSPMatrix *out) {
  return inp->RightApply(this,out);
}





inline TComplex& TDSPMatrix::operator()(UInt_t row, UInt_t col) {
  if (fVec) {
    if ((row<fRows)&&(row>=0)&&(col<fCols)&&(col>=0))
      return Element(row,col);
    else 
      Error("operator()","index (%d,%d) out of bounds (%d..%d,%d..%d)",row,col,0,fRows-1,0,fCols-1);
  } else 
    Error("operator()","matrix not allocated (this->fVec=0x0)!");
  return ComplexZero;
}

inline Bool_t TDSPMatrix::DiagZero() {
  if (!IsSquare()) return kTRUE;
  for(register UInt_t i=0;i<fRows; ++i) 
    if (Element(i,i)==0.0) return kTRUE;
  return kFALSE;
}

inline Bool_t TDSPMatrix::DiagZeros() {
  Int_t min=TMath::Min(fRows,fCols);
  for(register Int_t i=0;i<min; ++i) 
    if (Element(i,i)!=0.0) return kFALSE;
  return kTRUE;
}

inline Bool_t TDSPMatrix::IsDiag() {
  for(register UInt_t i=0;i<fRows;++i)
    for(register UInt_t j=0;j<fCols;++j)
      if ((i!=j)&&(Element(i,j)!=0.0)) return kFALSE;
  return kTRUE;
}

inline TDSPVector* TDSPMatrix::GetDiag(TDSPVector *a) {
  Int_t m=TMath::Min(fRows,fCols);
  if (!a) a=new TDSPVector();
  a->SetLen(m);
  for(register Int_t i=0;i<m;++i)
    a->Element(i)=Element(i,i);
  return a;
}


inline TDSPMatrix* TDSPMatrix::Transposed(TDSPMatrix *r) {
  
  if (r==this) { 
    Error("Transposed","Can not be the same matrix - use TDSPMatrix::TransposeMe() !");
    return NULL;
  }

  if (!r) r=new TDSPMatrix(); 
  r->SetSize(fCols,fRows);
  for(register UInt_t i=0;i<fRows;++i)
    for(register UInt_t j=0;j<fCols;++j)
      r->Element(j,i)=Element(i,j);
  return r;

}

inline void TDSPMatrix::TransposeMe() {
  if (fCols!=fRows) 
    Error("TransposeMe","non square matrix - use TDSPMatrix::Transposed()!");
  else
    for(register UInt_t i=0;    i<fRows;++i)
      for(register UInt_t j=i+1;j<fCols;++j)
    Swap(j,i,i,j);
}

inline TDSPMatrix* TDSPMatrix::Adjoint(TDSPMatrix *r) {
  
  if (!r) r=new TDSPMatrix(fCols,fRows); 

  for(register UInt_t i=0;i<fRows;++i)
    for(register UInt_t j=0;j<fCols;++j) {
      r->Element(j,i).fRe = Element(i,j).fRe;
      r->Element(j,i).fIm = -Element(i,j).fIm;
    }

  return r;

}
inline Double_t Dist(TDSPMatrix *s1, TDSPMatrix *s2) {
  return s1->Dist(s2);
}
inline Double_t Dist2(TDSPMatrix *s1, TDSPMatrix *s2) {
  return s1->Dist2(s2);
}
inline Double_t Dot(TDSPMatrix *s1, TDSPMatrix *s2) {
  return s1->Dot(s2);
}
inline Double_t Dot2(TDSPMatrix *s1, TDSPMatrix *s2) {
  return s1->Dot2(s2);
}
inline Double_t Dist(TDSPMatrix& s1, TDSPMatrix& s2) {
  return s1.Dist(s2);
}
inline Double_t Dist2(TDSPMatrix& s1, TDSPMatrix& s2) {
  return s1.Dist2(s2);
}
inline Double_t Dot(TDSPMatrix& s1, TDSPMatrix& s2) {
  return s1.Dot(s2);
}
inline Double_t Dot2(TDSPMatrix& s1, TDSPMatrix& s2) {
  return s1.Dot2(s2);
}
inline Bool_t IsDiag(TDSPMatrix*s) { return s->IsDiag();}; 
inline Bool_t IsDiag(TDSPMatrix&s) { return s.IsDiag();}; 


inline TDSPMatrix*  TDSPMatrix::Copy(TDSPMatrix *output){
  
  if (!output) output = new TDSPMatrix();
  output->SetSize(fRows,fCols);

  TComplex *ovec = output->GetVec();

  for(register Int_t i=0;i<Num;++i)  
    ovec[i] = fVec[i];

  return output;  
}

//  Swaps the Elements and return a reference to the second one element which contains now the first element 
//
inline TComplex& TDSPMatrix::Swap(Int_t i, Int_t j, Int_t m, Int_t n) {
#ifdef LINALG_MAT_ROWWISE
  UInt_t    l1=i*fCols+j;
  UInt_t    l2=m*fCols+n;
#else
  UInt_t    l1=j*fRows+i;
  UInt_t    l2=n*fRows+m;
#endif
  TComplex tmp=fVec[l1];
  fVec[l1] = fVec[l2];
  fVec[l2] = tmp;
  return fVec[l2];
}

#endif

---