00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "TMatrixT.h"
00025 #include "TMatrixTSym.h"
00026 #include "TMatrixTLazy.h"
00027 #include "TMath.h"
00028
00029 #ifndef R__ALPHA
00030 templateClassImp(TMatrixTLazy)
00031 templateClassImp(TMatrixTSymLazy)
00032 templateClassImp(THaarMatrixT)
00033 templateClassImp(THilbertMatrixT)
00034 templateClassImp(THilbertMatrixTSym)
00035 #endif
00036
00037
00038 template<class Element>
00039 THaarMatrixT<Element>::THaarMatrixT(Int_t order,Int_t no_cols)
00040 : TMatrixTLazy<Element>(1<<order, no_cols == 0 ? 1<<order : no_cols)
00041 {
00042 if (order <= 0)
00043 Error("THaarMatrixT","Haar order(%d) should be > 0",order);
00044 if (no_cols < 0)
00045 Error("THaarMatrixT","#cols(%d) in Haar should be >= 0",no_cols);
00046 }
00047
00048
00049 template<class Element>
00050 void MakeHaarMat(TMatrixT<Element> &m)
00051 {
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 R__ASSERT(m.IsValid());
00062 const Int_t no_rows = m.GetNrows();
00063 const Int_t no_cols = m.GetNcols();
00064
00065 if (no_rows < no_cols) {
00066 Error("MakeHaarMat","#rows(%d) should be >= #cols(%d)",no_rows,no_cols);
00067 return;
00068 }
00069 if (no_cols <= 0) {
00070 Error("MakeHaarMat","#cols(%d) should be > 0",no_cols);
00071 return;
00072 }
00073
00074
00075
00076
00077 TMatrixT<Element> mtr(no_cols,no_rows);
00078 Element *cp = mtr.GetMatrixArray();
00079 const Element *m_end = mtr.GetMatrixArray()+no_rows*no_cols;
00080
00081 Element norm_factor = 1/TMath::Sqrt((Element)no_rows);
00082
00083
00084 Int_t j;
00085 for (j = 0; j < no_rows; j++)
00086 *cp++ = norm_factor;
00087
00088
00089
00090
00091
00092 Int_t step_length = no_rows/2;
00093 while (cp < m_end && step_length > 0) {
00094 for (Int_t step_position = 0; cp < m_end && step_position < no_rows;
00095 step_position += 2*step_length, cp += no_rows) {
00096 Element *ccp = cp+step_position;
00097 for (j = 0; j < step_length; j++)
00098 *ccp++ = norm_factor;
00099 for (j = 0; j < step_length; j++)
00100 *ccp++ = -norm_factor;
00101 }
00102 step_length /= 2;
00103 norm_factor *= TMath::Sqrt(2.0);
00104 }
00105
00106 R__ASSERT(step_length != 0 || cp == m_end);
00107 R__ASSERT(no_rows != no_cols || step_length == 0);
00108
00109 m.Transpose(mtr);
00110 }
00111
00112
00113 template<class Element>
00114 void THaarMatrixT<Element>::FillIn(TMatrixT<Element> &m) const
00115 {
00116 MakeHaarMat(m);
00117 }
00118
00119
00120 template<class Element>
00121 THilbertMatrixT<Element>::THilbertMatrixT(Int_t no_rows,Int_t no_cols)
00122 : TMatrixTLazy<Element>(no_rows,no_cols)
00123 {
00124 if (no_rows <= 0)
00125 Error("THilbertMatrixT","#rows(%d) in Hilbert should be > 0",no_rows);
00126 if (no_cols <= 0)
00127 Error("THilbertMatrixT","#cols(%d) in Hilbert should be > 0",no_cols);
00128 }
00129
00130
00131 template<class Element>
00132 THilbertMatrixT<Element>::THilbertMatrixT(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb)
00133 : TMatrixTLazy<Element>(row_lwb,row_upb,col_lwb,col_upb)
00134 {
00135 if (row_upb < row_lwb)
00136 Error("THilbertMatrixT","row_upb(%d) in Hilbert should be >= row_lwb(%d)",row_upb,row_lwb);
00137 if (col_upb < col_lwb)
00138 Error("THilbertMatrixT","col_upb(%d) in Hilbert should be >= col_lwb(%d)",col_upb,col_lwb);
00139 }
00140
00141
00142 template<class Element>
00143 void MakeHilbertMat(TMatrixT<Element> &m)
00144 {
00145
00146
00147
00148 R__ASSERT(m.IsValid());
00149 const Int_t no_rows = m.GetNrows();
00150 const Int_t no_cols = m.GetNcols();
00151
00152 if (no_rows <= 0) {
00153 Error("MakeHilbertMat","#rows(%d) should be > 0",no_rows);
00154 return;
00155 }
00156 if (no_cols <= 0) {
00157 Error("MakeHilbertMat","#cols(%d) should be > 0",no_cols);
00158 return;
00159 }
00160
00161 Element *cp = m.GetMatrixArray();
00162 for (Int_t i = 0; i < no_rows; i++)
00163 for (Int_t j = 0; j < no_cols; j++)
00164 *cp++ = 1.0/(i+j+1.0);
00165 }
00166
00167
00168 template<class Element>
00169 void THilbertMatrixT<Element>::FillIn(TMatrixT<Element> &m) const
00170 {
00171 MakeHilbertMat(m);
00172 }
00173
00174
00175 template<class Element>
00176 THilbertMatrixTSym<Element>::THilbertMatrixTSym(Int_t no_rows)
00177 : TMatrixTSymLazy<Element>(no_rows)
00178 {
00179 if (no_rows <= 0)
00180 Error("THilbertMatrixTSym","#rows(%d) in Hilbert should be > 0",no_rows);
00181 }
00182
00183
00184 template<class Element>
00185 THilbertMatrixTSym<Element>::THilbertMatrixTSym(Int_t row_lwb,Int_t row_upb)
00186 : TMatrixTSymLazy<Element>(row_lwb,row_upb)
00187 {
00188 if (row_upb < row_lwb)
00189 Error("THilbertMatrixTSym","row_upb(%d) in Hilbert should be >= row_lwb(%d)",row_upb,row_lwb);
00190 }
00191
00192
00193 template<class Element>
00194 void MakeHilbertMat(TMatrixTSym<Element> &m)
00195 {
00196
00197
00198
00199 R__ASSERT(m.IsValid());
00200 const Int_t no_rows = m.GetNrows();
00201 if (no_rows <= 0) {
00202 Error("MakeHilbertMat","#rows(%d) should be > 0",no_rows);
00203 return;
00204 }
00205
00206 Element *cp = m.GetMatrixArray();
00207 for (Int_t i = 0; i < no_rows; i++)
00208 for (Int_t j = 0; j < no_rows; j++)
00209 *cp++ = 1.0/(i+j+1.0);
00210 }
00211
00212
00213 template<class Element>
00214 void THilbertMatrixTSym<Element>::FillIn(TMatrixTSym<Element> &m) const
00215 {
00216 MakeHilbertMat(m);
00217 }
00218
00219 template class TMatrixTLazy <Float_t>;
00220 template class TMatrixTSymLazy <Float_t>;
00221 template class THaarMatrixT <Float_t>;
00222 template class THilbertMatrixT <Float_t>;
00223 template class THilbertMatrixTSym<Float_t>;
00224
00225 template class TMatrixTLazy <Double_t>;
00226 template class TMatrixTSymLazy <Double_t>;
00227 template class THaarMatrixT <Double_t>;
00228 template class THilbertMatrixT <Double_t>;
00229 template class THilbertMatrixTSym<Double_t>;