StRoot: StarRoot/TRMatrix.cxx Source File

00001 #include <iomanip>
00002 #include "TRMatrix.h"
00003 #include "TMath.h"
00004 #include "TError.h"
00005 #include "TString.h"
00006 ClassImp(TRMatrix);
00007 //________________________________________________________________________________
00008 TRMatrix::TRMatrix(const TRMatrix &S, Int_t NI, Int_t NJ, Int_t I, Int_t J){
00009   if (NI == 0) NI = S.NI();
00010   if (NJ == 0) NJ = S.NJ();
00011   if (NI > S.NI()) NI = S.NI();
00012   if (NJ > S.NJ()) NJ = S.NJ();
00013   if (I == 0) {::Error("TRMatrix::TRMatrix(const TRMatrix &)", "index i %d out of bounds (size: %d, this: %p)", 
00014                        I, S.NI(), this); I = 1;}
00015   if (J == 0) {::Error("TRMatrix::TRMatrix(const TRMatrix &)", "index j %d out of bounds (size: %d, this: %p)", 
00016                        J, S.NJ(), this); J = 1;}
00017   if (I+NI-1 > S.NI()) {::Error("TRMatrix::TRMatrix(const TRMatrix &)", "index i %d out of bounds (size: %d, this: %p)", 
00018                        I+NI-1, S.NI(), this); I = 1;}
00019   if (J+NJ-1 > S.NJ()) {::Error("TRMatrix::TRMatrix(const TRMatrix &)", "index j %d out of bounds (size: %d, this: %p)", 
00020                        J+NJ-1, S.NJ(), this); J = 1;}
00021   fNrows = NI; 
00022   fNcols = NJ;
00023   Set(fNrows*fNcols);
00024   const Double_t *Array = S.GetArray();
00025   for (int i = 0; i < fNrows; i++) 
00026     for (int j = 0; j < fNcols; j++) 
00027       fArray[j + i*fNcols] = Array[J+j-1 + (I+i-1)*S.NJ()];
00028 }
00029 //________________________________________________________________________________
00030 TRMatrix &TRMatrix::operator=(const TRMatrix &rhs) {
00031     if (this != &rhs) SetMatrix(rhs.GetNrows(),rhs.GetNcols(),rhs.GetArray());
00032     return *this;
00033   }
00034 //________________________________________________________________________________
00035 TRMatrix::TRMatrix(Int_t nrows,Int_t ncols,Double_t a0, ...) : 
00036   TRArray(nrows*ncols), fNrows(nrows), fNcols(ncols) {
00037   __VA_LIST__(a0);
00038 }
00039 //________________________________________________________________________________
00040 void TRMatrix::SetMatrix(Int_t nrows,Int_t ncols,const  Double_t *array) {   
00041   fNrows = nrows; fNcols = ncols; Set(fNrows*fNcols,array);
00042 }
00043 //________________________________________________________________________________
00044 TRMatrix::TRMatrix(ETRMatrixCreatorsOp kop,Int_t nrows):
00045 TRArray(nrows*nrows), fNrows(nrows), fNcols(nrows) {
00046   switch (kop) {
00047   case kZero:
00048     break;
00049   case kUnit:
00050     for (int i=0; i<fNrows; i++) fArray[i+fNrows*i] = 1;
00051     break;
00052   default:
00053     Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
00054   }
00055 }
00056 //________________________________________________________________________________
00057 TRMatrix::TRMatrix(const TRMatrix& A, ETRMatrixCreatorsOp kop) {
00058   switch (kop) {
00059   case kTransposed:
00060     fNrows = A.GetNcols();
00061     fNcols = A.GetNrows();
00062     Set(fNrows*fNcols);
00063     TCL::mxtrp(A.GetArray(),fArray,fNcols,fNrows);
00064     break;
00065   default:
00066     Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
00067   }
00068 }
00069 //________________________________________________________________________________
00070 TRMatrix::TRMatrix(const TRMatrix& A, ETRMatrixCreatorsOp kop,const TRMatrix& B):  TRArray(0){
00071   Int_t NI, NJ, NK;
00072   switch (kop) {
00073   case kAxB:
00074     NI = A.GetNrows(); fNrows = NI;
00075     NJ = A.GetNcols();
00076     assert (NJ == B.GetNrows());
00077     NK = B.GetNcols(); fNcols = NK;
00078     Set(NI*NK);
00079     TCL::mxmpy(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00080     break;
00081   case kAxBT:
00082     NI = A.GetNrows(); fNrows = NI;
00083     NJ = A.GetNcols();
00084     assert (NJ == B.GetNcols());
00085     NK = B.GetNrows(); fNcols = NK;
00086     Set(NI*NK);
00087     TCL::mxmpy1(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00088     break;
00089   case kATxB:
00090     NI = A.GetNcols(); fNrows = NI;
00091     NJ = A.GetNrows();
00092     assert (NJ == B.GetNrows());
00093     NK = B.GetNcols(); fNcols = NK;
00094     Set(NI*NK);
00095     TCL::mxmpy2(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00096     break;
00097   case kATxBT:
00098     NI = A.GetNcols(); fNrows = NI;
00099     NJ = A.GetNrows();
00100     assert (NJ == B.GetNcols());
00101     NK = B.GetNrows(); fNcols = NK;
00102     Set(NI*NK);
00103     TCL::mxmpy3(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00104     break;
00105 #if 1
00106   case kAxBxAT: //[NI,NJ]x[NJ,NJ]x[NI,NJ]T => [NI,NI]
00107     NI = A.GetNcols();
00108     NJ = A.GetNrows();
00109     assert (NJ == B.GetNcols());
00110     Set(NI*NI);
00111     TCL::mxmlrt(A.GetArray(),B.GetArray(),fArray,NI,NJ);
00112     break;
00113   case kATxBxA: //[NI,NJ]Tx[NI,NI]x[NI,NJ] => [NJ,NJ]
00114     NI = A.GetNcols(); 
00115     NJ = A.GetNrows(); fNrows = NJ; fNcols = NJ;
00116     assert (NI == B.GetNcols());
00117     Set(NJ*NJ);
00118     TCL::mxmltr(A.GetArray(),B.GetArray(),fArray,NJ,NI);
00119     break;
00120 #endif
00121   default:
00122     Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
00123     
00124   }
00125 }
00126 //________________________________________________________________________________
00127 void TRMatrix::Add(const TRMatrix& A, ETRMatrixCreatorsOp kop,const TRMatrix& B) {
00128   Int_t NI, NJ, NK;
00129   switch (kop) {
00130   case kAxB:
00131     NI = A.GetNrows();
00132     NJ = A.GetNcols();
00133     assert (NJ == B.GetNrows());
00134     NK = B.GetNcols(); 
00135     assert(NI == fNrows && NK == fNcols);
00136     TCL::mxmad(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00137     break;
00138   case kAxBT:
00139     NI = A.GetNrows();
00140     NJ = A.GetNcols();
00141     assert (NJ == B.GetNcols());
00142     NK = B.GetNrows(); 
00143     assert(NI == fNrows && NK == fNcols);
00144     TCL::mxmad1(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00145     break;
00146   case kATxB:
00147     NI = A.GetNcols();
00148     NJ = A.GetNrows();
00149     assert (NJ == B.GetNrows());
00150     NK = B.GetNcols(); 
00151     assert(NI == fNrows && NK == fNcols);
00152     TCL::mxmad2(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00153     break;
00154   case kATxBT:
00155     NI = A.GetNcols();
00156     NJ = A.GetNrows();
00157     assert (NJ == B.GetNcols());
00158     NK = B.GetNrows(); 
00159     assert(NI == fNrows && NK == fNcols);
00160     TCL::mxmad3(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00161     break;
00162     
00163   default:
00164     Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
00165   }
00166 }
00167 //________________________________________________________________________________
00168 void TRMatrix::Substruct(const TRMatrix& A, ETRMatrixCreatorsOp kop,const TRMatrix& B) {
00169   Int_t NI, NJ, NK;
00170   switch (kop) {
00171   case kAxB:
00172     NI = A.GetNrows();
00173     NJ = A.GetNcols();
00174     assert (NJ == B.GetNrows());
00175     NK = B.GetNcols(); 
00176     assert(NI == fNrows && NK == fNcols);
00177     TCL::mxmub(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00178     break;
00179   case kAxBT:
00180     NI = A.GetNrows();
00181     NJ = A.GetNcols();
00182     assert (NJ == B.GetNcols());
00183     NK = B.GetNrows(); 
00184     assert(NI == fNrows && NK == fNcols);
00185     TCL::mxmub1(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00186     break;
00187   case kATxB:
00188     NI = A.GetNcols();
00189     NJ = A.GetNrows();
00190     assert (NJ == B.GetNrows());
00191     NK = B.GetNcols(); 
00192     assert(NI == fNrows && NK == fNcols);
00193     TCL::mxmub2(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00194     break;
00195   case kATxBT: //[NJ,NI]T x [NK,NJ]T => [NI,NK]
00196     NI = A.GetNcols();
00197     NJ = A.GetNrows();
00198     assert (NJ == B.GetNcols());
00199     NK = B.GetNrows(); 
00200     assert(NI == fNrows && NK == fNcols);
00201     TCL::mxmub3(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
00202     break;
00203     
00204   default:
00205     Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
00206   }
00207 }
00208 //________________________________________________________________________________
00209 TRMatrix::TRMatrix(const TRSymMatrix &S, ETRMatrixCreatorsOp kop,const TRMatrix& A){
00210   Int_t M,N; // A[M,N], B[N,M]], C[M,N], S[M,M], R[N,N]
00211   switch (kop) { 
00212   case kSxA: // S[M,M]*A[M,N] => C[M,N]
00213     M = A.GetNrows();
00214     N = A.GetNcols();
00215     assert(M == S.GetNrows());
00216     fNrows = M;
00217     fNcols = N;
00218     Set(fNrows*fNcols);
00219     TCL::trsa(S.GetArray(),A.GetArray(),fArray,M,N);
00220     break;
00221   case kSxAT: // S[M,M]*BT[N,M] => C[M,N]
00222     N = A.GetNrows();
00223     M = A.GetNcols();
00224     assert(M == S.GetNrows());
00225     fNrows = M;
00226     fNcols = N;
00227     Set(fNrows*fNcols);
00228     TCL::trsat(S.GetArray(),A.GetArray(),fArray,M,N);
00229     break;
00230   default:
00231     Error("TRMatrix SxA (ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
00232   }
00233 }
00234 //________________________________________________________________________________
00235 TRMatrix::TRMatrix(const TRMatrix& A, ETRMatrixCreatorsOp kop,const TRSymMatrix &S){
00236   Int_t M,N; // A[M,N], B[N,M]], C[M,N], S[M,M], R[N,N]
00237   switch (kop) { 
00238   case kAxS: // A[M,N]*R[N,N] => C[M,N]
00239     M = A.GetNrows();
00240     N = A.GetNcols();
00241     assert(N == S.GetNrows());
00242     fNrows = M;
00243     fNcols = N;
00244     Set(fNrows*fNcols);
00245     TCL::tras(A.GetArray(),S.GetArray(),fArray,M,N);
00246     break;
00247   case kATxS: // BT[N,M]*R[N,N] => C[M,N]
00248     M = A.GetNcols();
00249     N = A.GetNrows();
00250     assert(N == S.GetNrows());
00251     fNrows = M;
00252     fNcols = N;
00253     Set(fNrows*fNcols);
00254     TCL::trats(A.GetArray(),S.GetArray(),fArray,M,N);
00255     break;
00256   default:
00257     Error("TRMatrix AxS (ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
00258   }
00259 }
00260 //________________________________________________________________________________
00261 TRMatrix::TRMatrix(const TRSymMatrix &S) {
00262   Int_t M = S.GetNcols();
00263   fNrows = M;
00264   fNcols = M;
00265   Set(fNrows*fNcols);
00266   TCL::trupck(S.GetArray(),fArray,M);
00267 }
00268 //________________________________________________________________________________
00269 ostream& operator<<(ostream& s,const TRMatrix &target) {
00270   Int_t Nrows = target.GetNrows();
00271   Int_t Ncols = target.GetNcols();
00272   const Double_t *Array = target.GetArray();
00273   s << "Rectangular Matrix Size \t[" << Nrows << "," << Ncols << "]" << endl;
00274   if (Array) { 
00275     s.setf(ios::fixed,ios::scientific);
00276     s.setf(ios::showpos);
00277     for (int i = 0; i< Nrows; i++) {
00278       Int_t Nzeros = 0;
00279       for (int j = Ncols-1; j >= 0; j--) if (Array[j + i*Ncols] == 0.0) {Nzeros++;} else break;
00280       if (Nzeros == 1) Nzeros = 0;
00281       for (int j = 0; j < Ncols-Nzeros; j++) s << Form("%10.3f",Array[j + i*Ncols]);
00282       if (Nzeros) s << Form("%8i*0",Nzeros);
00283       s << endl;
00284     }
00285     s.unsetf(ios::showpos);
00286   }
00287   else s << " Empty";
00288   return s;
00289 }
00290 //________________________________________________________________________________
00291 void TRMatrix::Print(Option_t *opt) const {if (opt) {}; cout << *this << endl;}