StRoot: StarRoot/TRSymMatrix.h Source File

00001 #ifndef ROOT_TRSymMatrix
00002 #define ROOT_TRSymMatrix
00003 #include "TRArray.h"
00004 #include <assert.h>
00005 class TRMatrix;
00006 class TRVector;
00007 class TDiagMatrix;
00008 class TRSymMatrix : public TRArray {
00009  public:
00010   TRSymMatrix(Int_t nrows=0) : TRArray(nrows*(nrows+1)/2), fNrows(nrows) {}
00011   TRSymMatrix(Int_t nrows,const Double_t *Array);
00012   TRSymMatrix(Int_t nrows,const Float_t *Array);
00013   TRSymMatrix(Int_t nrows,const Char_t *s);
00014   TRSymMatrix(const TRSymMatrix& W,ETRMatrixCreatorsOp kop);
00015   TRSymMatrix(ETRMatrixCreatorsOp kop,Int_t nrows);
00016   TRSymMatrix(const TRMatrix& A);
00017   TRSymMatrix(const TRMatrix& A,ETRMatrixCreatorsOp kop,const TRSymMatrix& S);
00018   TRSymMatrix(const TRSymMatrix& Q,ETRMatrixCreatorsOp kop,const TRSymMatrix& T);
00019   TRSymMatrix(const TRMatrix& A,ETRMatrixCreatorsOp kop);
00020 #ifndef __CINT__
00021   TRSymMatrix (Int_t nrows, Double_t a0, ...);
00022 #endif
00023   virtual ~TRSymMatrix() {}
00024   void  Inverse() { TCL::trsinv(fArray,fArray, fNrows);}
00025   Int_t GetNrows()  const       {return fNrows;} 
00026   Int_t GetNcols()  const       {return GetNrows();}
00027   virtual ETRMatrixType GetMatrixType() const {return kSemiPosDefinedSymMatrix;}
00028   virtual Double_t Product(const TRVector& A,ETRMatrixCreatorsOp kop=kAxSxAT);
00029   virtual Int_t    SpmInv(const TRSymMatrix &S, TRVector *B = 0);
00030   static  Int_t    spminv(Double_t *v, Double_t *b, Int_t n, 
00031                           Int_t &nrank, Double_t *diag, Bool_t *flag);
00032   virtual void    Print(Option_t *opt="") const;
00033   Double_t       &operator()(Int_t i)         {return TRArray::operator[](i);}
00034   Double_t        operator()(Int_t i) const   {return TRArray::operator[](i);}
00035   Double_t       &operator()(Int_t i,Int_t j);
00036   Double_t        operator()(Int_t i,Int_t j) const;
00037   void AddRow(const Double_t *row) {
00038     fNrows++; Set(fNrows*(fNrows+1)/2); memcpy(fArray+(fNrows-1)*fNrows/2, row, fNrows*sizeof(Double_t));
00039   }
00040   void AddRow(const Double_t row) {//Diagonal element only
00041     fNrows++; Set(fNrows*(fNrows+1)/2); fArray[(fNrows+1)*fNrows/2-1] = row;
00042   }
00043   static Int_t TrsInv(const Double_t *g, Double_t *gi, Int_t n);
00044   static Int_t TrInv(const Double_t *g, Double_t *gi, Int_t n);
00045   static Int_t TrchLU(const Double_t *g, Double_t *gi, Int_t n);
00046   static Int_t TrsmUL(const Double_t *g, Double_t *gi, Int_t n);
00047  protected:
00048   Int_t     fNrows;            // number of rows 
00049  public:
00050   ClassDef(TRSymMatrix,1)  // TRSymMatrix class (double precision)
00051 };
00052 ostream& operator<<(ostream& s,const TRSymMatrix &target);
00053 inline Double_t &TRSymMatrix::operator()(Int_t i,Int_t j){
00054   //  assert(! (j < 0 || j >= fNrows));
00055   if (j < 0 || j >= fNrows) {
00056     ::Error("TRSymMatrix::operator()", "index j %d out of bounds (size: %d, this: %p)", 
00057             j, fNrows, this); 
00058     j = 0;
00059     assert(0);
00060   }
00061   //  assert(! (i < 0 || i >= fNrows));
00062   if (i < 0 || i >= fNrows) {
00063     ::Error("TRSymMatrix::operator()", "index i %d out of bounds (size: %d, this: %p)", 
00064             i, fNrows, this); 
00065     i = 0;
00066     assert(0);
00067   }
00068   Int_t m = i;
00069   Int_t l = j;
00070   if (i > j) {m = j; l = i;}
00071   return TArrayD::operator[](m + (l+1)*l/2);
00072 }
00073 inline Double_t TRSymMatrix::operator()(Int_t i,Int_t j) const {
00074   //  assert(! (j < 0 || j >= fNrows));
00075   if (j < 0 || j >= fNrows) {
00076     ::Error("TRSymMatrix::operator()", "index j %d out of bounds (size: %d, this: %p)", 
00077             j, fNrows, this); 
00078     j = 0;
00079     assert(0);
00080   }
00081   //  assert(! (i < 0 || i >= fNrows));
00082   if (i < 0 || i >= fNrows) {
00083     ::Error("TRSymMatrix::operator()", "index i %d out of bounds (size: %d, this: %p)", 
00084             i, fNrows, this); 
00085     i = 0;
00086     assert(0);
00087   }
00088   Int_t m = i;
00089   Int_t l = j;
00090   if (i > j) {m = j; l = i;}
00091   return TArrayD::operator[](m + (l+1)*l/2);
00092 }
00093 #endif
00094