// @(#)root/matrix:$Id$
// Authors: Fons Rademakers, Eddy Offermann  Nov 2003

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// Templates of Lazy Matrix classes.                                    //
//                                                                      //
//   TMatrixTLazy                                                       //
//   TMatrixTSymLazy                                                    //
//   THaarMatrixT                                                       //
//   THilbertMatrixT                                                    //
//   THilbertMatrixTSym                                                 //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "TMatrixT.h"
#include "TMatrixTSym.h"
#include "TMatrixTLazy.h"
#include "TMath.h"

templateClassImp(TMatrixTLazy)
templateClassImp(TMatrixTSymLazy)
templateClassImp(THaarMatrixT)
templateClassImp(THilbertMatrixT)
templateClassImp(THilbertMatrixTSym)


//______________________________________________________________________________
template<class Element>
THaarMatrixT<Element>::THaarMatrixT(Int_t order,Int_t no_cols)
    : TMatrixTLazy<Element>(1<<order, no_cols == 0 ? 1<<order : no_cols)
{
   if (order <= 0)
      Error("THaarMatrixT","Haar order(%d) should be > 0",order);
   if (no_cols < 0)
      Error("THaarMatrixT","#cols(%d) in Haar should be >= 0",no_cols);
}

//______________________________________________________________________________
template<class Element>
void MakeHaarMat(TMatrixT<Element> &m)
{
   // Create an orthonormal (2^n)*(no_cols) Haar (sub)matrix, whose columns
   // are Haar functions. If no_cols is 0, create the complete matrix with
   // 2^n columns. Example, the complete Haar matrix of the second order is:
   // column 1: [ 1  1  1  1]/2
   // column 2: [ 1  1 -1 -1]/2
   // column 3: [ 1 -1  0  0]/sqrt(2)
   // column 4: [ 0  0  1 -1]/sqrt(2)
   // Matrix m is assumed to be zero originally.

   R__ASSERT(m.IsValid());
   const Int_t no_rows = m.GetNrows();
   const Int_t no_cols = m.GetNcols();

   if (no_rows < no_cols) {
      Error("MakeHaarMat","#rows(%d) should be >= #cols(%d)",no_rows,no_cols);
      return;
   }
   if (no_cols <= 0) {
      Error("MakeHaarMat","#cols(%d) should be > 0",no_cols);
      return;
   }

   // It is easier to calculate a Haar matrix when the elements are stored
   // column-wise . Since we are row-wise, the transposed Haar is calculated

   TMatrixT<Element> mtr(no_cols,no_rows);
         Element *cp    = mtr.GetMatrixArray();
   const Element *m_end = mtr.GetMatrixArray()+no_rows*no_cols;

   Element norm_factor = 1/TMath::Sqrt((Element)no_rows);

   // First row is always 1 (up to normalization)
   Int_t j;
   for (j = 0; j < no_rows; j++)
      *cp++ = norm_factor;

   // The other functions are kind of steps: stretch of 1 followed by the
   // equally long stretch of -1. The functions can be grouped in families
   // according to their order (step size), differing only in the location
   // of the step
   Int_t step_length = no_rows/2;
   while (cp < m_end && step_length > 0) {
      for (Int_t step_position = 0; cp < m_end && step_position < no_rows;
              step_position += 2*step_length, cp += no_rows) {
         Element *ccp = cp+step_position;
         for (j = 0; j < step_length; j++)
            *ccp++ = norm_factor;
         for (j = 0; j < step_length; j++)
            *ccp++ = -norm_factor;
      }
      step_length /= 2;
      norm_factor *= TMath::Sqrt(2.0);
   }

   R__ASSERT(step_length != 0       || cp == m_end);
   R__ASSERT(no_rows     != no_cols || step_length == 0);

   m.Transpose(mtr);
}

//______________________________________________________________________________
template<class Element>
void THaarMatrixT<Element>::FillIn(TMatrixT<Element> &m) const
{
   MakeHaarMat(m);
}

//______________________________________________________________________________
template<class Element>
THilbertMatrixT<Element>::THilbertMatrixT(Int_t no_rows,Int_t no_cols)
    : TMatrixTLazy<Element>(no_rows,no_cols)
{
   if (no_rows <= 0)
      Error("THilbertMatrixT","#rows(%d) in Hilbert should be > 0",no_rows);
   if (no_cols <= 0)
      Error("THilbertMatrixT","#cols(%d) in Hilbert should be > 0",no_cols);
}

//______________________________________________________________________________
template<class Element>
THilbertMatrixT<Element>::THilbertMatrixT(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb)
    : TMatrixTLazy<Element>(row_lwb,row_upb,col_lwb,col_upb)
{
   if (row_upb < row_lwb)
      Error("THilbertMatrixT","row_upb(%d) in Hilbert should be >= row_lwb(%d)",row_upb,row_lwb);
   if (col_upb < col_lwb)
      Error("THilbertMatrixT","col_upb(%d) in Hilbert should be >= col_lwb(%d)",col_upb,col_lwb);
}

//______________________________________________________________________________
template<class Element>
void MakeHilbertMat(TMatrixT<Element> &m)
{
   // Make a Hilbert matrix. Hilb[i,j] = 1/(i+j+1),
   // i,j=0...max-1 (matrix need not be a square one).

   R__ASSERT(m.IsValid());
   const Int_t no_rows = m.GetNrows();
   const Int_t no_cols = m.GetNcols();

   if (no_rows <= 0) {
      Error("MakeHilbertMat","#rows(%d) should be > 0",no_rows);
      return;
   }
   if (no_cols <= 0) {
      Error("MakeHilbertMat","#cols(%d) should be > 0",no_cols);
      return;
   }

   Element *cp = m.GetMatrixArray();
   for (Int_t i = 0; i < no_rows; i++)
      for (Int_t j = 0; j < no_cols; j++)
         *cp++ = 1.0/(i+j+1.0);
}

//______________________________________________________________________________
template<class Element>
void THilbertMatrixT<Element>::FillIn(TMatrixT<Element> &m) const
{
   MakeHilbertMat(m);
}

//______________________________________________________________________________
template<class Element>
THilbertMatrixTSym<Element>::THilbertMatrixTSym(Int_t no_rows)
    : TMatrixTSymLazy<Element>(no_rows)
{
   if (no_rows <= 0)
      Error("THilbertMatrixTSym","#rows(%d) in Hilbert should be > 0",no_rows);
}

//______________________________________________________________________________
template<class Element>
THilbertMatrixTSym<Element>::THilbertMatrixTSym(Int_t row_lwb,Int_t row_upb)
    : TMatrixTSymLazy<Element>(row_lwb,row_upb)
{
   if (row_upb < row_lwb)
      Error("THilbertMatrixTSym","row_upb(%d) in Hilbert should be >= row_lwb(%d)",row_upb,row_lwb);
}

//______________________________________________________________________________
template<class Element>
void MakeHilbertMat(TMatrixTSym<Element> &m)
{
   // Make a Hilbert matrix. Hilb[i,j] = 1/(i+j+1),
   // i,j=0...max-1 (matrix must be square).

   R__ASSERT(m.IsValid());
   const Int_t no_rows = m.GetNrows();
   if (no_rows <= 0) {
      Error("MakeHilbertMat","#rows(%d) should be > 0",no_rows);
      return;
   }

   Element *cp = m.GetMatrixArray();
   for (Int_t i = 0; i < no_rows; i++)
      for (Int_t j = 0; j < no_rows; j++)
         *cp++ = 1.0/(i+j+1.0);
}

//______________________________________________________________________________
template<class Element>
void THilbertMatrixTSym<Element>::FillIn(TMatrixTSym<Element> &m) const
{
   MakeHilbertMat(m);
}

template class TMatrixTLazy      <Float_t>;
template class TMatrixTSymLazy   <Float_t>;
template class THaarMatrixT      <Float_t>;
template class THilbertMatrixT   <Float_t>;
template class THilbertMatrixTSym<Float_t>;

template class TMatrixTLazy      <Double_t>;
template class TMatrixTSymLazy   <Double_t>;
template class THaarMatrixT      <Double_t>;
template class THilbertMatrixT   <Double_t>;
template class THilbertMatrixTSym<Double_t>;
 TMatrixTLazy.cxx:1
 TMatrixTLazy.cxx:2
 TMatrixTLazy.cxx:3
 TMatrixTLazy.cxx:4
 TMatrixTLazy.cxx:5
 TMatrixTLazy.cxx:6
 TMatrixTLazy.cxx:7
 TMatrixTLazy.cxx:8
 TMatrixTLazy.cxx:9
 TMatrixTLazy.cxx:10
 TMatrixTLazy.cxx:11
 TMatrixTLazy.cxx:12
 TMatrixTLazy.cxx:13
 TMatrixTLazy.cxx:14
 TMatrixTLazy.cxx:15
 TMatrixTLazy.cxx:16
 TMatrixTLazy.cxx:17
 TMatrixTLazy.cxx:18
 TMatrixTLazy.cxx:19
 TMatrixTLazy.cxx:20
 TMatrixTLazy.cxx:21
 TMatrixTLazy.cxx:22
 TMatrixTLazy.cxx:23
 TMatrixTLazy.cxx:24
 TMatrixTLazy.cxx:25
 TMatrixTLazy.cxx:26
 TMatrixTLazy.cxx:27
 TMatrixTLazy.cxx:28
 TMatrixTLazy.cxx:29
 TMatrixTLazy.cxx:30
 TMatrixTLazy.cxx:31
 TMatrixTLazy.cxx:32
 TMatrixTLazy.cxx:33
 TMatrixTLazy.cxx:34
 TMatrixTLazy.cxx:35
 TMatrixTLazy.cxx:36
 TMatrixTLazy.cxx:37
 TMatrixTLazy.cxx:38
 TMatrixTLazy.cxx:39
 TMatrixTLazy.cxx:40
 TMatrixTLazy.cxx:41
 TMatrixTLazy.cxx:42
 TMatrixTLazy.cxx:43
 TMatrixTLazy.cxx:44
 TMatrixTLazy.cxx:45
 TMatrixTLazy.cxx:46
 TMatrixTLazy.cxx:47
 TMatrixTLazy.cxx:48
 TMatrixTLazy.cxx:49
 TMatrixTLazy.cxx:50
 TMatrixTLazy.cxx:51
 TMatrixTLazy.cxx:52
 TMatrixTLazy.cxx:53
 TMatrixTLazy.cxx:54
 TMatrixTLazy.cxx:55
 TMatrixTLazy.cxx:56
 TMatrixTLazy.cxx:57
 TMatrixTLazy.cxx:58
 TMatrixTLazy.cxx:59
 TMatrixTLazy.cxx:60
 TMatrixTLazy.cxx:61
 TMatrixTLazy.cxx:62
 TMatrixTLazy.cxx:63
 TMatrixTLazy.cxx:64
 TMatrixTLazy.cxx:65
 TMatrixTLazy.cxx:66
 TMatrixTLazy.cxx:67
 TMatrixTLazy.cxx:68
 TMatrixTLazy.cxx:69
 TMatrixTLazy.cxx:70
 TMatrixTLazy.cxx:71
 TMatrixTLazy.cxx:72
 TMatrixTLazy.cxx:73
 TMatrixTLazy.cxx:74
 TMatrixTLazy.cxx:75
 TMatrixTLazy.cxx:76
 TMatrixTLazy.cxx:77
 TMatrixTLazy.cxx:78
 TMatrixTLazy.cxx:79
 TMatrixTLazy.cxx:80
 TMatrixTLazy.cxx:81
 TMatrixTLazy.cxx:82
 TMatrixTLazy.cxx:83
 TMatrixTLazy.cxx:84
 TMatrixTLazy.cxx:85
 TMatrixTLazy.cxx:86
 TMatrixTLazy.cxx:87
 TMatrixTLazy.cxx:88
 TMatrixTLazy.cxx:89
 TMatrixTLazy.cxx:90
 TMatrixTLazy.cxx:91
 TMatrixTLazy.cxx:92
 TMatrixTLazy.cxx:93
 TMatrixTLazy.cxx:94
 TMatrixTLazy.cxx:95
 TMatrixTLazy.cxx:96
 TMatrixTLazy.cxx:97
 TMatrixTLazy.cxx:98
 TMatrixTLazy.cxx:99
 TMatrixTLazy.cxx:100
 TMatrixTLazy.cxx:101
 TMatrixTLazy.cxx:102
 TMatrixTLazy.cxx:103
 TMatrixTLazy.cxx:104
 TMatrixTLazy.cxx:105
 TMatrixTLazy.cxx:106
 TMatrixTLazy.cxx:107
 TMatrixTLazy.cxx:108
 TMatrixTLazy.cxx:109
 TMatrixTLazy.cxx:110
 TMatrixTLazy.cxx:111
 TMatrixTLazy.cxx:112
 TMatrixTLazy.cxx:113
 TMatrixTLazy.cxx:114
 TMatrixTLazy.cxx:115
 TMatrixTLazy.cxx:116
 TMatrixTLazy.cxx:117
 TMatrixTLazy.cxx:118
 TMatrixTLazy.cxx:119
 TMatrixTLazy.cxx:120
 TMatrixTLazy.cxx:121
 TMatrixTLazy.cxx:122
 TMatrixTLazy.cxx:123
 TMatrixTLazy.cxx:124
 TMatrixTLazy.cxx:125
 TMatrixTLazy.cxx:126
 TMatrixTLazy.cxx:127
 TMatrixTLazy.cxx:128
 TMatrixTLazy.cxx:129
 TMatrixTLazy.cxx:130
 TMatrixTLazy.cxx:131
 TMatrixTLazy.cxx:132
 TMatrixTLazy.cxx:133
 TMatrixTLazy.cxx:134
 TMatrixTLazy.cxx:135
 TMatrixTLazy.cxx:136
 TMatrixTLazy.cxx:137
 TMatrixTLazy.cxx:138
 TMatrixTLazy.cxx:139
 TMatrixTLazy.cxx:140
 TMatrixTLazy.cxx:141
 TMatrixTLazy.cxx:142
 TMatrixTLazy.cxx:143
 TMatrixTLazy.cxx:144
 TMatrixTLazy.cxx:145
 TMatrixTLazy.cxx:146
 TMatrixTLazy.cxx:147
 TMatrixTLazy.cxx:148
 TMatrixTLazy.cxx:149
 TMatrixTLazy.cxx:150
 TMatrixTLazy.cxx:151
 TMatrixTLazy.cxx:152
 TMatrixTLazy.cxx:153
 TMatrixTLazy.cxx:154
 TMatrixTLazy.cxx:155
 TMatrixTLazy.cxx:156
 TMatrixTLazy.cxx:157
 TMatrixTLazy.cxx:158
 TMatrixTLazy.cxx:159
 TMatrixTLazy.cxx:160
 TMatrixTLazy.cxx:161
 TMatrixTLazy.cxx:162
 TMatrixTLazy.cxx:163
 TMatrixTLazy.cxx:164
 TMatrixTLazy.cxx:165
 TMatrixTLazy.cxx:166
 TMatrixTLazy.cxx:167
 TMatrixTLazy.cxx:168
 TMatrixTLazy.cxx:169
 TMatrixTLazy.cxx:170
 TMatrixTLazy.cxx:171
 TMatrixTLazy.cxx:172
 TMatrixTLazy.cxx:173
 TMatrixTLazy.cxx:174
 TMatrixTLazy.cxx:175
 TMatrixTLazy.cxx:176
 TMatrixTLazy.cxx:177
 TMatrixTLazy.cxx:178
 TMatrixTLazy.cxx:179
 TMatrixTLazy.cxx:180
 TMatrixTLazy.cxx:181
 TMatrixTLazy.cxx:182
 TMatrixTLazy.cxx:183
 TMatrixTLazy.cxx:184
 TMatrixTLazy.cxx:185
 TMatrixTLazy.cxx:186
 TMatrixTLazy.cxx:187
 TMatrixTLazy.cxx:188
 TMatrixTLazy.cxx:189
 TMatrixTLazy.cxx:190
 TMatrixTLazy.cxx:191
 TMatrixTLazy.cxx:192
 TMatrixTLazy.cxx:193
 TMatrixTLazy.cxx:194
 TMatrixTLazy.cxx:195
 TMatrixTLazy.cxx:196
 TMatrixTLazy.cxx:197
 TMatrixTLazy.cxx:198
 TMatrixTLazy.cxx:199
 TMatrixTLazy.cxx:200
 TMatrixTLazy.cxx:201
 TMatrixTLazy.cxx:202
 TMatrixTLazy.cxx:203
 TMatrixTLazy.cxx:204
 TMatrixTLazy.cxx:205
 TMatrixTLazy.cxx:206
 TMatrixTLazy.cxx:207
 TMatrixTLazy.cxx:208
 TMatrixTLazy.cxx:209
 TMatrixTLazy.cxx:210
 TMatrixTLazy.cxx:211
 TMatrixTLazy.cxx:212
 TMatrixTLazy.cxx:213
 TMatrixTLazy.cxx:214
 TMatrixTLazy.cxx:215
 TMatrixTLazy.cxx:216
 TMatrixTLazy.cxx:217
 TMatrixTLazy.cxx:218
 TMatrixTLazy.cxx:219
 TMatrixTLazy.cxx:220
 TMatrixTLazy.cxx:221
 TMatrixTLazy.cxx:222
 TMatrixTLazy.cxx:223
 TMatrixTLazy.cxx:224
 TMatrixTLazy.cxx:225
 TMatrixTLazy.cxx:226
 TMatrixTLazy.cxx:227
 TMatrixTLazy.cxx:228