Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TMatrixTLazy.cxx
Go to the documentation of this file.
1 // @(#)root/matrix:$Id$
2 // Authors: Fons Rademakers, Eddy Offermann Nov 2003
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class TMatrixTLazy
13  \ingroup Matrix
14 
15  Templates of Lazy Matrix classes.
16 ~~~
17  TMatrixTLazy
18  TMatrixTSymLazy
19  THaarMatrixT
20  THilbertMatrixT
21  THilbertMatrixTSym
22 ~~~
23 */
24 
25 #include "TMatrixT.h"
26 #include "TMatrixTSym.h"
27 #include "TMatrixTLazy.h"
28 #include "TMath.h"
29 
30 templateClassImp(TMatrixTLazy);
31 templateClassImp(TMatrixTSymLazy);
32 templateClassImp(THaarMatrixT);
33 templateClassImp(THilbertMatrixT);
34 templateClassImp(THilbertMatrixTSym);
35 
36 
37 ////////////////////////////////////////////////////////////////////////////////
38 
39 template<class Element>
40 THaarMatrixT<Element>::THaarMatrixT(Int_t order,Int_t no_cols)
41  : TMatrixTLazy<Element>(1<<order, no_cols == 0 ? 1<<order : no_cols)
42 {
43  if (order <= 0)
44  Error("THaarMatrixT","Haar order(%d) should be > 0",order);
45  if (no_cols < 0)
46  Error("THaarMatrixT","#cols(%d) in Haar should be >= 0",no_cols);
47 }
48 
49 ////////////////////////////////////////////////////////////////////////////////
50 /// Create an orthonormal (2^n)*(no_cols) Haar (sub)matrix, whose columns
51 /// are Haar functions. If no_cols is 0, create the complete matrix with
52 /// 2^n columns. Example, the complete Haar matrix of the second order is:
53 /// column 1: [ 1 1 1 1]/2
54 /// column 2: [ 1 1 -1 -1]/2
55 /// column 3: [ 1 -1 0 0]/sqrt(2)
56 /// column 4: [ 0 0 1 -1]/sqrt(2)
57 /// Matrix m is assumed to be zero originally.
58 
59 template<class Element>
60 void MakeHaarMat(TMatrixT<Element> &m)
61 {
62  R__ASSERT(m.IsValid());
63  const Int_t no_rows = m.GetNrows();
64  const Int_t no_cols = m.GetNcols();
65 
66  if (no_rows < no_cols) {
67  Error("MakeHaarMat","#rows(%d) should be >= #cols(%d)",no_rows,no_cols);
68  return;
69  }
70  if (no_cols <= 0) {
71  Error("MakeHaarMat","#cols(%d) should be > 0",no_cols);
72  return;
73  }
74 
75  // It is easier to calculate a Haar matrix when the elements are stored
76  // column-wise . Since we are row-wise, the transposed Haar is calculated
77 
78  TMatrixT<Element> mtr(no_cols,no_rows);
79  Element *cp = mtr.GetMatrixArray();
80  const Element *m_end = mtr.GetMatrixArray()+no_rows*no_cols;
81 
82  Element norm_factor = 1/TMath::Sqrt((Element)no_rows);
83 
84  // First row is always 1 (up to normalization)
85  Int_t j;
86  for (j = 0; j < no_rows; j++)
87  *cp++ = norm_factor;
88 
89  // The other functions are kind of steps: stretch of 1 followed by the
90  // equally long stretch of -1. The functions can be grouped in families
91  // according to their order (step size), differing only in the location
92  // of the step
93  Int_t step_length = no_rows/2;
94  while (cp < m_end && step_length > 0) {
95  for (Int_t step_position = 0; cp < m_end && step_position < no_rows;
96  step_position += 2*step_length, cp += no_rows) {
97  Element *ccp = cp+step_position;
98  for (j = 0; j < step_length; j++)
99  *ccp++ = norm_factor;
100  for (j = 0; j < step_length; j++)
101  *ccp++ = -norm_factor;
102  }
103  step_length /= 2;
104  norm_factor *= TMath::Sqrt(2.0);
105  }
106 
107  R__ASSERT(step_length != 0 || cp == m_end);
108  R__ASSERT(no_rows != no_cols || step_length == 0);
109 
110  m.Transpose(mtr);
111 }
112 
113 ////////////////////////////////////////////////////////////////////////////////
114 
115 template<class Element>
116 void THaarMatrixT<Element>::FillIn(TMatrixT<Element> &m) const
117 {
118  MakeHaarMat(m);
119 }
120 
121 ////////////////////////////////////////////////////////////////////////////////
122 
123 template<class Element>
124 THilbertMatrixT<Element>::THilbertMatrixT(Int_t no_rows,Int_t no_cols)
125  : TMatrixTLazy<Element>(no_rows,no_cols)
126 {
127  if (no_rows <= 0)
128  Error("THilbertMatrixT","#rows(%d) in Hilbert should be > 0",no_rows);
129  if (no_cols <= 0)
130  Error("THilbertMatrixT","#cols(%d) in Hilbert should be > 0",no_cols);
131 }
132 
133 ////////////////////////////////////////////////////////////////////////////////
134 
135 template<class Element>
136 THilbertMatrixT<Element>::THilbertMatrixT(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb)
137  : TMatrixTLazy<Element>(row_lwb,row_upb,col_lwb,col_upb)
138 {
139  if (row_upb < row_lwb)
140  Error("THilbertMatrixT","row_upb(%d) in Hilbert should be >= row_lwb(%d)",row_upb,row_lwb);
141  if (col_upb < col_lwb)
142  Error("THilbertMatrixT","col_upb(%d) in Hilbert should be >= col_lwb(%d)",col_upb,col_lwb);
143 }
144 
145 ////////////////////////////////////////////////////////////////////////////////
146 /// Make a Hilbert matrix. Hilb[i,j] = 1/(i+j+1),
147 /// i,j=0...max-1 (matrix need not be a square one).
148 
149 template<class Element>
150 void MakeHilbertMat(TMatrixT<Element> &m)
151 {
152  R__ASSERT(m.IsValid());
153  const Int_t no_rows = m.GetNrows();
154  const Int_t no_cols = m.GetNcols();
155 
156  if (no_rows <= 0) {
157  Error("MakeHilbertMat","#rows(%d) should be > 0",no_rows);
158  return;
159  }
160  if (no_cols <= 0) {
161  Error("MakeHilbertMat","#cols(%d) should be > 0",no_cols);
162  return;
163  }
164 
165  Element *cp = m.GetMatrixArray();
166  for (Int_t i = 0; i < no_rows; i++)
167  for (Int_t j = 0; j < no_cols; j++)
168  *cp++ = 1.0/(i+j+1.0);
169 }
170 
171 ////////////////////////////////////////////////////////////////////////////////
172 
173 template<class Element>
174 void THilbertMatrixT<Element>::FillIn(TMatrixT<Element> &m) const
175 {
176  MakeHilbertMat(m);
177 }
178 
179 ////////////////////////////////////////////////////////////////////////////////
180 
181 template<class Element>
182 THilbertMatrixTSym<Element>::THilbertMatrixTSym(Int_t no_rows)
183  : TMatrixTSymLazy<Element>(no_rows)
184 {
185  if (no_rows <= 0)
186  Error("THilbertMatrixTSym","#rows(%d) in Hilbert should be > 0",no_rows);
187 }
188 
189 ////////////////////////////////////////////////////////////////////////////////
190 
191 template<class Element>
192 THilbertMatrixTSym<Element>::THilbertMatrixTSym(Int_t row_lwb,Int_t row_upb)
193  : TMatrixTSymLazy<Element>(row_lwb,row_upb)
194 {
195  if (row_upb < row_lwb)
196  Error("THilbertMatrixTSym","row_upb(%d) in Hilbert should be >= row_lwb(%d)",row_upb,row_lwb);
197 }
198 
199 ////////////////////////////////////////////////////////////////////////////////
200 /// Make a Hilbert matrix. Hilb[i,j] = 1/(i+j+1),
201 /// i,j=0...max-1 (matrix must be square).
202 
203 template<class Element>
204 void MakeHilbertMat(TMatrixTSym<Element> &m)
205 {
206  R__ASSERT(m.IsValid());
207  const Int_t no_rows = m.GetNrows();
208  if (no_rows <= 0) {
209  Error("MakeHilbertMat","#rows(%d) should be > 0",no_rows);
210  return;
211  }
212 
213  Element *cp = m.GetMatrixArray();
214  for (Int_t i = 0; i < no_rows; i++)
215  for (Int_t j = 0; j < no_rows; j++)
216  *cp++ = 1.0/(i+j+1.0);
217 }
218 
219 ////////////////////////////////////////////////////////////////////////////////
220 
221 template<class Element>
222 void THilbertMatrixTSym<Element>::FillIn(TMatrixTSym<Element> &m) const
223 {
224  MakeHilbertMat(m);
225 }
226 
227 template class TMatrixTLazy <Float_t>;
228 template class TMatrixTSymLazy <Float_t>;
229 template class THaarMatrixT <Float_t>;
230 template class THilbertMatrixT <Float_t>;
231 template class THilbertMatrixTSym<Float_t>;
232 
233 template class TMatrixTLazy <Double_t>;
234 template class TMatrixTSymLazy <Double_t>;
235 template class THaarMatrixT <Double_t>;
236 template class THilbertMatrixT <Double_t>;
237 template class THilbertMatrixTSym<Double_t>;