StRoot: StarRoot/TRDiagMatrix.cxx Source File

00001 #include <iomanip>
00002 #include "TRVector.h"
00003 #include "TRDiagMatrix.h"
00004 #include "TError.h"
00005 ClassImp(TRDiagMatrix);
00006 //________________________________________________________________________________
00007 TRDiagMatrix::TRDiagMatrix(Int_t nrows,Double_t a0, ...) : 
00008   TRArray(nrows), fNrows(nrows) {
00009   __VA_LIST__(a0);
00010 }
00011 //________________________________________________________________________________
00012 TRDiagMatrix::TRDiagMatrix(const TRDiagMatrix& S,ETRMatrixCreatorsOp kop) {
00013   Int_t i;
00014   switch (kop) {
00015   case kInvertedPosDef:
00016   case kInverted:
00017   case kInvertedA:
00018     fNrows = S.GetNcols();
00019     Set(fNrows*(fNrows+1)/2);
00020     for (i = 0; i < fNrows;  i++) if (fArray[i] != 0) fArray[i] = 1./fArray[i];
00021     break;
00022   default:
00023     Error("TRDiagMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
00024   }
00025 }
00026 //________________________________________________________________________________
00027 Double_t TRDiagMatrix::Product(const TRVector& A,ETRMatrixCreatorsOp kop) {
00028   Int_t M, N; // N == 1
00029   Double_t Value = 0;
00030   Int_t i;
00031   switch (kop) { // 
00032   case kAxSxAT: //A[M,N]*D[N,N]*AT[M,N] => R[M,M]; 
00033   case kATxSxA: //BT[N,M]*S[N,N]*B[N,M] => R[M,M]; 
00034     M = A.GetNrows();
00035     N = GetNrows();
00036     assert(N == A.GetNcols() || M == N);
00037     for (i = 0; i < N; i++) Value += A[i]*fArray[i]*A[i];
00038     break;
00039     break;
00040   default:
00041     Error("TRDiagMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
00042   }
00043   return Value;
00044 }
00045 //________________________________________________________________________________
00046 ostream& operator<<(ostream& s,const TRDiagMatrix &target) {
00047   static const int width = 10;
00048   static const Double_t zero = 0;
00049   Int_t Nrows = target.GetNrows();
00050   const Double_t *Array = target.GetArray();
00051   s << "Diagonal Matrix Size \t[" 
00052     << Nrows << "," << Nrows << "]" << endl;
00053   if (Array) {
00054     s.setf(ios::fixed,ios::scientific);
00055     s.setf(ios::showpos);
00056     for (int i = 0; i< Nrows; i++) {
00057       for (int j = 0; j <= i; j++)
00058         if (i == j)
00059           s << std::setw(width) << std::setprecision(width-3) << Array[i] << ":\t";
00060         else
00061           s << std::setw(width) << std::setprecision(width-3) << zero << ":\t";
00062       s << endl;
00063     }
00064     s.unsetf(ios::showpos);
00065   }
00066   else s << " Empty";   
00067   return s;
00068 }
00069 //________________________________________________________________________________
00070 void TRDiagMatrix::Print(Option_t *opt) const {if (opt) {}; cout << *this << endl;}