TMatrixTLazy.cxx

Go to the documentation of this file.
00001 // @(#)root/matrix:$Id: TMatrixTLazy.cxx 20882 2007-11-19 11:31:26Z rdm $
00002 // Authors: Fons Rademakers, Eddy Offermann  Nov 2003
00003 
00004 /*************************************************************************
00005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
00006  * All rights reserved.                                                  *
00007  *                                                                       *
00008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
00009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
00010  *************************************************************************/
00011 
00012 //////////////////////////////////////////////////////////////////////////
00013 //                                                                      //
00014 // Templates of Lazy Matrix classes.                                    //
00015 //                                                                      //
00016 //   TMatrixTLazy                                                       //
00017 //   TMatrixTSymLazy                                                    //
00018 //   THaarMatrixT                                                       //
00019 //   THilbertMatrixT                                                    //
00020 //   THilbertMatrixTSym                                                 //
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    // Create an orthonormal (2^n)*(no_cols) Haar (sub)matrix, whose columns
00053    // are Haar functions. If no_cols is 0, create the complete matrix with
00054    // 2^n columns. Example, the complete Haar matrix of the second order is:
00055    // column 1: [ 1  1  1  1]/2
00056    // column 2: [ 1  1 -1 -1]/2
00057    // column 3: [ 1 -1  0  0]/sqrt(2)
00058    // column 4: [ 0  0  1 -1]/sqrt(2)
00059    // Matrix m is assumed to be zero originally.
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    // It is easier to calculate a Haar matrix when the elements are stored
00075    // column-wise . Since we are row-wise, the transposed Haar is calculated
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    // First row is always 1 (up to normalization)
00084    Int_t j;
00085    for (j = 0; j < no_rows; j++)
00086       *cp++ = norm_factor;
00087 
00088    // The other functions are kind of steps: stretch of 1 followed by the
00089    // equally long stretch of -1. The functions can be grouped in families
00090    // according to their order (step size), differing only in the location
00091    // of the step
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    // Make a Hilbert matrix. Hilb[i,j] = 1/(i+j+1),
00146    // i,j=0...max-1 (matrix need not be a square one).
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    // Make a Hilbert matrix. Hilb[i,j] = 1/(i+j+1),
00197    // i,j=0...max-1 (matrix must be square).
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>;

Generated on Tue Jul 5 14:36:34 2011 for ROOT_528-00b_version by  doxygen 1.5.1