StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TRMatrix.cxx
1 #include <iomanip>
2 #include "TRMatrix.h"
3 #include "TMath.h"
4 #include "TError.h"
5 #include "TString.h"
6 ClassImp(TRMatrix);
7 //________________________________________________________________________________
8 TRMatrix::TRMatrix(const TRMatrix &S, Int_t NI, Int_t NJ, Int_t I, Int_t J){
9  if (NI == 0) NI = S.NI();
10  if (NJ == 0) NJ = S.NJ();
11  if (NI > S.NI()) NI = S.NI();
12  if (NJ > S.NJ()) NJ = S.NJ();
13  if (I == 0) {::Error("TRMatrix::TRMatrix(const TRMatrix &)", "index i %d out of bounds (size: %d, this: %p)",
14  I, S.NI(), this); I = 1;}
15  if (J == 0) {::Error("TRMatrix::TRMatrix(const TRMatrix &)", "index j %d out of bounds (size: %d, this: %p)",
16  J, S.NJ(), this); J = 1;}
17  if (I+NI-1 > S.NI()) {::Error("TRMatrix::TRMatrix(const TRMatrix &)", "index i %d out of bounds (size: %d, this: %p)",
18  I+NI-1, S.NI(), this); I = 1;}
19  if (J+NJ-1 > S.NJ()) {::Error("TRMatrix::TRMatrix(const TRMatrix &)", "index j %d out of bounds (size: %d, this: %p)",
20  J+NJ-1, S.NJ(), this); J = 1;}
21  fNrows = NI;
22  fNcols = NJ;
23  Set(fNrows*fNcols);
24  const Double_t *Array = S.GetArray();
25  for (int i = 0; i < fNrows; i++)
26  for (int j = 0; j < fNcols; j++)
27  fArray[j + i*fNcols] = Array[J+j-1 + (I+i-1)*S.NJ()];
28 }
29 //________________________________________________________________________________
30 TRMatrix &TRMatrix::operator=(const TRMatrix &rhs) {
31  if (this != &rhs) SetMatrix(rhs.GetNrows(),rhs.GetNcols(),rhs.GetArray());
32  return *this;
33  }
34 //________________________________________________________________________________
35 TRMatrix::TRMatrix(Int_t nrows,Int_t ncols,Double_t a0, ...) :
36  TRArray(nrows*ncols), fNrows(nrows), fNcols(ncols) {
37  __VA_LIST__(a0);
38 }
39 //________________________________________________________________________________
40 void TRMatrix::SetMatrix(Int_t nrows,Int_t ncols,const Double_t *array) {
41  fNrows = nrows; fNcols = ncols; Set(fNrows*fNcols,array);
42 }
43 //________________________________________________________________________________
44 TRMatrix::TRMatrix(ETRMatrixCreatorsOp kop,Int_t nrows):
45 TRArray(nrows*nrows), fNrows(nrows), fNcols(nrows) {
46  switch (kop) {
47  case kZero:
48  break;
49  case kUnit:
50  for (int i=0; i<fNrows; i++) fArray[i+fNrows*i] = 1;
51  break;
52  default:
53  Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
54  }
55 }
56 //________________________________________________________________________________
57 TRMatrix::TRMatrix(const TRMatrix& A, ETRMatrixCreatorsOp kop) {
58  switch (kop) {
59  case kTransposed:
60  fNrows = A.GetNcols();
61  fNcols = A.GetNrows();
62  Set(fNrows*fNcols);
63  TCL::mxtrp(A.GetArray(),fArray,fNcols,fNrows);
64  break;
65  default:
66  Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
67  }
68 }
69 //________________________________________________________________________________
70 TRMatrix::TRMatrix(const TRMatrix& A, ETRMatrixCreatorsOp kop,const TRMatrix& B): TRArray(0){
71  Int_t NI, NJ, NK;
72  switch (kop) {
73  case kAxB:
74  NI = A.GetNrows(); fNrows = NI;
75  NJ = A.GetNcols();
76  assert (NJ == B.GetNrows());
77  NK = B.GetNcols(); fNcols = NK;
78  Set(NI*NK);
79  TCL::mxmpy(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
80  break;
81  case kAxBT:
82  NI = A.GetNrows(); fNrows = NI;
83  NJ = A.GetNcols();
84  assert (NJ == B.GetNcols());
85  NK = B.GetNrows(); fNcols = NK;
86  Set(NI*NK);
87  TCL::mxmpy1(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
88  break;
89  case kATxB:
90  NI = A.GetNcols(); fNrows = NI;
91  NJ = A.GetNrows();
92  assert (NJ == B.GetNrows());
93  NK = B.GetNcols(); fNcols = NK;
94  Set(NI*NK);
95  TCL::mxmpy2(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
96  break;
97  case kATxBT:
98  NI = A.GetNcols(); fNrows = NI;
99  NJ = A.GetNrows();
100  assert (NJ == B.GetNcols());
101  NK = B.GetNrows(); fNcols = NK;
102  Set(NI*NK);
103  TCL::mxmpy3(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
104  break;
105 #if 1
106  case kAxBxAT: //[NI,NJ]x[NJ,NJ]x[NI,NJ]T => [NI,NI]
107  NI = A.GetNcols();
108  NJ = A.GetNrows();
109  assert (NJ == B.GetNcols());
110  Set(NI*NI);
111  TCL::mxmlrt(A.GetArray(),B.GetArray(),fArray,NI,NJ);
112  break;
113  case kATxBxA: //[NI,NJ]Tx[NI,NI]x[NI,NJ] => [NJ,NJ]
114  NI = A.GetNcols();
115  NJ = A.GetNrows(); fNrows = NJ; fNcols = NJ;
116  assert (NI == B.GetNcols());
117  Set(NJ*NJ);
118  TCL::mxmltr(A.GetArray(),B.GetArray(),fArray,NJ,NI);
119  break;
120 #endif
121  default:
122  Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
123 
124  }
125 }
126 //________________________________________________________________________________
127 void TRMatrix::Add(const TRMatrix& A, ETRMatrixCreatorsOp kop,const TRMatrix& B) {
128  Int_t NI, NJ, NK;
129  switch (kop) {
130  case kAxB:
131  NI = A.GetNrows();
132  NJ = A.GetNcols();
133  assert (NJ == B.GetNrows());
134  NK = B.GetNcols();
135  assert(NI == fNrows && NK == fNcols);
136  TCL::mxmad(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
137  break;
138  case kAxBT:
139  NI = A.GetNrows();
140  NJ = A.GetNcols();
141  assert (NJ == B.GetNcols());
142  NK = B.GetNrows();
143  assert(NI == fNrows && NK == fNcols);
144  TCL::mxmad1(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
145  break;
146  case kATxB:
147  NI = A.GetNcols();
148  NJ = A.GetNrows();
149  assert (NJ == B.GetNrows());
150  NK = B.GetNcols();
151  assert(NI == fNrows && NK == fNcols);
152  TCL::mxmad2(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
153  break;
154  case kATxBT:
155  NI = A.GetNcols();
156  NJ = A.GetNrows();
157  assert (NJ == B.GetNcols());
158  NK = B.GetNrows();
159  assert(NI == fNrows && NK == fNcols);
160  TCL::mxmad3(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
161  break;
162 
163  default:
164  Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
165  }
166 }
167 //________________________________________________________________________________
168 void TRMatrix::Substruct(const TRMatrix& A, ETRMatrixCreatorsOp kop,const TRMatrix& B) {
169  Int_t NI, NJ, NK;
170  switch (kop) {
171  case kAxB:
172  NI = A.GetNrows();
173  NJ = A.GetNcols();
174  assert (NJ == B.GetNrows());
175  NK = B.GetNcols();
176  assert(NI == fNrows && NK == fNcols);
177  TCL::mxmub(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
178  break;
179  case kAxBT:
180  NI = A.GetNrows();
181  NJ = A.GetNcols();
182  assert (NJ == B.GetNcols());
183  NK = B.GetNrows();
184  assert(NI == fNrows && NK == fNcols);
185  TCL::mxmub1(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
186  break;
187  case kATxB:
188  NI = A.GetNcols();
189  NJ = A.GetNrows();
190  assert (NJ == B.GetNrows());
191  NK = B.GetNcols();
192  assert(NI == fNrows && NK == fNcols);
193  TCL::mxmub2(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
194  break;
195  case kATxBT: //[NJ,NI]T x [NK,NJ]T => [NI,NK]
196  NI = A.GetNcols();
197  NJ = A.GetNrows();
198  assert (NJ == B.GetNcols());
199  NK = B.GetNrows();
200  assert(NI == fNrows && NK == fNcols);
201  TCL::mxmub3(A.GetArray(),B.GetArray(),fArray,NI,NJ,NK);
202  break;
203 
204  default:
205  Error("TRMatrix(ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
206  }
207 }
208 //________________________________________________________________________________
209 TRMatrix::TRMatrix(const TRSymMatrix &S, ETRMatrixCreatorsOp kop,const TRMatrix& A){
210  Int_t M,N; // A[M,N], B[N,M]], C[M,N], S[M,M], R[N,N]
211  switch (kop) {
212  case kSxA: // S[M,M]*A[M,N] => C[M,N]
213  M = A.GetNrows();
214  N = A.GetNcols();
215  assert(M == S.GetNrows());
216  fNrows = M;
217  fNcols = N;
218  Set(fNrows*fNcols);
219  TCL::trsa(S.GetArray(),A.GetArray(),fArray,M,N);
220  break;
221  case kSxAT: // S[M,M]*BT[N,M] => C[M,N]
222  N = A.GetNrows();
223  M = A.GetNcols();
224  assert(M == S.GetNrows());
225  fNrows = M;
226  fNcols = N;
227  Set(fNrows*fNcols);
228  TCL::trsat(S.GetArray(),A.GetArray(),fArray,M,N);
229  break;
230  default:
231  Error("TRMatrix SxA (ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
232  }
233 }
234 //________________________________________________________________________________
235 TRMatrix::TRMatrix(const TRMatrix& A, ETRMatrixCreatorsOp kop,const TRSymMatrix &S){
236  Int_t M,N; // A[M,N], B[N,M]], C[M,N], S[M,M], R[N,N]
237  switch (kop) {
238  case kAxS: // A[M,N]*R[N,N] => C[M,N]
239  M = A.GetNrows();
240  N = A.GetNcols();
241  assert(N == S.GetNrows());
242  fNrows = M;
243  fNcols = N;
244  Set(fNrows*fNcols);
245  TCL::tras(A.GetArray(),S.GetArray(),fArray,M,N);
246  break;
247  case kATxS: // BT[N,M]*R[N,N] => C[M,N]
248  M = A.GetNcols();
249  N = A.GetNrows();
250  assert(N == S.GetNrows());
251  fNrows = M;
252  fNcols = N;
253  Set(fNrows*fNcols);
254  TCL::trats(A.GetArray(),S.GetArray(),fArray,M,N);
255  break;
256  default:
257  Error("TRMatrix AxS (ETRMatrixCreatorsOp)", "operation %d not yet implemented", kop);
258  }
259 }
260 //________________________________________________________________________________
261 TRMatrix::TRMatrix(const TRSymMatrix &S) {
262  Int_t M = S.GetNcols();
263  fNrows = M;
264  fNcols = M;
265  Set(fNrows*fNcols);
266  TCL::trupck(S.GetArray(),fArray,M);
267 }
268 //________________________________________________________________________________
269 std::ostream& operator<<(std::ostream& s,const TRMatrix &target) {
270  Int_t Nrows = target.GetNrows();
271  Int_t Ncols = target.GetNcols();
272  const Double_t *Array = target.GetArray();
273  s << "Rectangular Matrix Size \t[" << Nrows << "," << Ncols << "]" << std::endl;
274  if (Array) {
275  s.setf(std::ios::fixed,std::ios::scientific);
276  s.setf(std::ios::showpos);
277  for (int i = 0; i< Nrows; i++) {
278  Int_t Nzeros = 0;
279  for (int j = Ncols-1; j >= 0; j--) if (Array[j + i*Ncols] == 0.0) {Nzeros++;} else break;
280  if (Nzeros == 1) Nzeros = 0;
281  for (int j = 0; j < Ncols-Nzeros; j++) s << Form("%10.3f",Array[j + i*Ncols]);
282  if (Nzeros) s << Form("%8i*0",Nzeros);
283  s << std::endl;
284  }
285  s.unsetf(std::ios::showpos);
286  }
287  else s << " Empty";
288  return s;
289 }
290 //________________________________________________________________________________
291 void TRMatrix::Print(Option_t *opt) const {if (opt) {}; std::cout << *this << std::endl;}
Error
Definition: FJcore.h:367
TCL::mxtrp
static float * mxtrp(const float *a, float *b, int i, int j)
Definition: TCernLib.cxx:230
TCL::trsat
static float * trsat(const float *s, const float *a, float *b, int m, int n)
Definition: TCernLib.cxx:1298
TCL::trsa
static float * trsa(const float *s, const float *a, float *b, int m, int n)
Definition: TCernLib.cxx:1111
TCL::trats
static float * trats(const float *a, const float *s, float *b, int m, int n)
Definition: TCernLib.cxx:629
TCL::tras
static float * tras(const float *a, const float *s, float *b, int m, int n)
Definition: TCernLib.cxx:470
TCL::trupck
static float * trupck(const float *u, float *s, int m)
Definition: TCernLib.cxx:1248
Generated by   doxygen 1.8.5