Logo ROOT  
Reference Guide
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
35
36
37////////////////////////////////////////////////////////////////////////////////
38
39template<class Element>
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
59template<class Element>
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
115template<class Element>
117{
118 MakeHaarMat(m);
119}
120
121////////////////////////////////////////////////////////////////////////////////
122
123template<class Element>
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
135template<class Element>
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
149template<class Element>
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
173template<class Element>
175{
177}
178
179////////////////////////////////////////////////////////////////////////////////
180
181template<class Element>
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
191template<class Element>
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
203template<class Element>
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
221template<class Element>
223{
225}
226
227template class TMatrixTLazy <Float_t>;
228template class TMatrixTSymLazy <Float_t>;
229template class THaarMatrixT <Float_t>;
230template class THilbertMatrixT <Float_t>;
231template class THilbertMatrixTSym<Float_t>;
232
233template class TMatrixTLazy <Double_t>;
234template class TMatrixTSymLazy <Double_t>;
235template class THaarMatrixT <Double_t>;
236template class THilbertMatrixT <Double_t>;
237template class THilbertMatrixTSym<Double_t>;
int Int_t
Definition: RtypesCore.h:45
#define templateClassImp(name)
Definition: Rtypes.h:408
#define R__ASSERT(e)
Definition: TError.h:118
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
void MakeHaarMat(TMatrixT< Element > &m)
Create an orthonormal (2^n)*(no_cols) Haar (sub)matrix, whose columns are Haar functions.
void MakeHilbertMat(TMatrixT< Element > &m)
Make a Hilbert matrix.
void FillIn(TMatrixT< Element > &m) const
void FillIn(TMatrixTSym< Element > &m) const
void FillIn(TMatrixT< Element > &m) const
Templates of Lazy Matrix classes.
Definition: TMatrixTLazy.h:43
TMatrixTSym.
Definition: TMatrixTSym.h:34
TMatrixT.
Definition: TMatrixT.h:39
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
Double_t Sqrt(Double_t x)
Definition: TMath.h:641
auto * m
Definition: textangle.C:8