---

DenseMatrix.h

/****************************************************************************
 *                                                                          *
 *  Program: FDDlib                                                         *
 *  Version: 1.1                                                            *
 *                                                                          *
 *  Copyright (C) 2002 - 2004 by Eric Miller and Dana Brooks                *
 *  All rights reserved.                                                    *
 *                                                                          *
 *  This software is Version 1.1 of the fddlib tomography toolbox.          *
 *  It is not to be redistributed or used for any commercial purpose        *
 *  without the prior written consent of the authors and Northeastern       *
 *  University.                                                             *
 *                                                                          *
 *  This software is provided as is, and any express or implied warranties, *
 *  including but not limited to the implied warranty of merchantability    *
 *  and the implied warranty of fitness for a particular purpose, are dis-  *
 *  claimed.  In no event shall the authors or Northeastern University be   *
 *  liable for any direct, indirect, incidental, special, exemplary, or     *
 *  consequential damages (including but not limited to procurement of      * 
 *  substitute goods or services; loss of use, data, or profits; or busi-   *
 *  ness interruption) however caused and on any theory of liability,       *
 *  whether in contract, strict liability, or tort (including negligence or *
 *  otherwise) arising in any way from the use of this software, even if    *
 *  advised of the possibility of such damage.                              *
 *                                                                          *
 *  Portions of this code benefit from ideas from Kyle Guilbert, Greg       *
 *  Boverman, Derek Uluski, David Kaeli, and Jennifer Black                 *
 *                                                                          *
 ****************************************************************************/

#ifndef _DENSE_MATRIX_H_
#define _DENSE_MATRIX_H_

#include "RealVector.h"
#include "Matrix.h"
#include <string>

// Minimal tolerance for maintaining linear independence
#define QR_MIN_NORM 1.0e-20

namespace FDDlib {

template<class T>
class DenseMatrix : public Matrix<T>
{
protected:
  int rows_;
  int cols_;
  T*  valVec_;
  
public:
  DenseMatrix(int rows, int cols);
  
  DenseMatrix();

  DenseMatrix(const DenseMatrix<T>& dm);

  virtual ~DenseMatrix();
  
  void deallocate();

  void init(int rows, int cols);

  DenseMatrix<T>& operator=(const DenseMatrix<T>& dm);

  DenseMatrix<T>& operator*=(const T v);

  template<class Vector>
  Vector operator*(const Vector& m) const;

  int getNumRows() const;

  int getNumCols() const;

  T val(int row, int col) const;

  void set(int row, int col, T v) throw(std::string);

  void QRDecomp(DenseMatrix<T>& Q, DenseMatrix<T>& R) const throw(std::string);

  template<class Vector>
  void LSSolve(Vector& x, const Vector& b) const throw(std::string);

  DenseMatrix<T> operator*(const DenseMatrix<T>& a) const;

  DenseMatrix<T> operator+(const DenseMatrix<T>& a) const;

  DenseMatrix<T> operator-(const DenseMatrix<T>& a) const;

  DenseMatrix<T> transpose() const;

  RealVector<T> getRow(int row) const;

  RealVector<T> getColumn(int col) const;

  T& operator()(int row, int col);

  T operator()(int row, int col) const;

};

template<class T>
DenseMatrix<T>::DenseMatrix(int rows, int cols)
{
  int i;

  rows_ = rows;
  cols_ = cols;
  valVec_ = new T[rows_*cols_];
  for (i = 0; i < rows_*cols_; i++)
  {
    valVec_[i] = (T) 0;
  }
}

template<class T>
DenseMatrix<T>::DenseMatrix()
{
  rows_ = 0;
  cols_ = 0;
  valVec_ = (T*) NULL;
}

template<class T>
DenseMatrix<T>::DenseMatrix(const DenseMatrix<T>& dm)
{
  int i, j;

  rows_ = dm.getNumRows();
  cols_ = dm.getNumCols();
  valVec_ = new T[rows_*cols_];
  for (i = 0; i <  rows_; i++)
  {
    for (j = 0; j < cols_; j++)
    {
      set(i,j, dm.val(i,j));
    }
  }
}

template<class T>
void DenseMatrix<T>::deallocate()
{
  delete[] valVec_;
}


template<class T>
DenseMatrix<T>::~DenseMatrix()
{
  delete[] valVec_;
}

template<class T>
void DenseMatrix<T>::init(int rows, int cols)
{
  int i;

  rows_ = rows;
  cols_ = cols;
  if (valVec_ != (T*) NULL)
  {
    delete[] valVec_;
  }

  valVec_ = new T[rows*cols];
  for (i = 0; i < rows_*cols_; i++)
  {
    valVec_[i] = (T) 0;
  }
}

template<class T>
DenseMatrix<T>& DenseMatrix<T>::operator=(const DenseMatrix<T>& dm)
{
  int i, j;

  if (valVec_ != (T*) NULL)
    delete[] valVec_;

  rows_ = dm.getNumRows();
  cols_ = dm.getNumCols();
  valVec_ = new T[rows_*cols_];
  for (i = 0; i < rows_; i++)
  {
    for (j =0; j < cols_; j++)
    {
      set(i,j,dm.val(i,j));
    }
  }
  return *this;
}

template<class T>
inline DenseMatrix<T>& DenseMatrix<T>::operator*=(const T v)
{
  int i, j;

  for (i = 0; i < rows_; i++)
  {
    for (j = 0; j < cols_; j++)
    {
      set(i,j, val(i,j) * v);
    }
  }

  return *this;
}

template<class T>
template<class Vector>
inline Vector DenseMatrix<T>::operator*(const Vector& m) const
{
  int    i, j;
  T      v;
  Vector ret(getNumRows());

  for (i = 0; i < getNumRows(); i++)
  {
    v = (T) 0;
    for (j = 0; j < getNumCols(); j++)
    {
     v += (T) (val(i,j) * m(j));
    }
    ret(i) = v;
  }

  return ret;
}

template<class T>
inline int DenseMatrix<T>::getNumRows() const
{
  return rows_;
}

template<class T>
inline int DenseMatrix<T>::getNumCols() const
{
  return cols_;
}

template<class T>
inline T DenseMatrix<T>::val(int row, int col) const
{
  return valVec_[row*cols_ + col];
}

template<class T>
inline void DenseMatrix<T>::set(int row, int col, T v) throw(std::string)
{
  valVec_[row*cols_ + col] = v;
}

template<class T>
void DenseMatrix<T>::QRDecomp(DenseMatrix<T>& Q, DenseMatrix<T>& R) const throw(std::string)
{
  int        i, j;
  T          v;
  double     n;
  T          proj;

  // What will be an orthonormal basis
  // of the columns of A
  RealVector<T> * rv = new RealVector<T>[getNumCols()];
  for (i = 0; i < getNumCols(); i++)
  {
    rv[i].init(getNumRows());
  }

  for (i = 0; i < getNumRows(); i++)
  {
    for (j = 0; j < getNumCols(); j++)
    {
      v = val(i,j);
      rv[j].set(i, v);
    }
  }

  // Now we assume that the columns of the matrix are linearly
  // independent

  Q.init(getNumRows(), getNumCols());
  R.init(getNumCols(), getNumCols());

  for (i = 0; i < getNumCols(); i++)
  {
    n = norm(rv[i]);
    if (n < QR_MIN_NORM)
      throw std::string("DenseMatrix::QRDecomp: norm of basis is too small");
    rv[i] *= (T) (1/n);
    R(i,i) = (T) n;
    for (j = i+1; j < getNumCols(); j++)
    {
      proj = (T) dot(rv[i], rv[j]);
      R(i,j) = proj;
      rv[j] -= proj*rv[i];
    }
  }

  // Now copy the linearly independent
  // basis back to Q

  for (i = 0; i < getNumRows(); i++)
  {
    for (j = 0; j < getNumCols(); j++)
    {
      v = rv[j].val(i);
      Q(i,j) = v;
    }
  }

}


template<class T>
template<class Vector>
void DenseMatrix<T>::LSSolve(Vector& x, const Vector& b) const throw(std::string)
{
  int            i, j;
  DenseMatrix<T> Q,R;
  Vector         rhs;
  T              v;

  try {
    QRDecomp(Q,R);
  } catch (std::string s) { 
    throw std::string("DenseMatrix::LSSolve: QR Decomposition failed.");
  }

  // Solution is Rx = Q'*b
  rhs = (Q.transpose())*b;

  Vector sol(R.getNumCols());
  for (i = 0;  i < R.getNumCols(); i++)
  {
    sol(i) = 0;
  }

  // Now do a backsolve of R, which is upper-triangular
  for (i = R.getNumCols()-1; i >= 0; i--)
  {
    v = rhs(i);
    for (j = i+1; j < R.getNumCols(); j++)
    {
      v -= (T) sol(j)*R(i,j);
    }
    v = v/( (T) R(i,i));
    sol(i) = v;
  }

  x = sol;
}

template<class T>
inline DenseMatrix<T> DenseMatrix<T>::operator*(const DenseMatrix<T>& a) const
{

  int i,j,k;
  T   v;
  DenseMatrix<T> dm(getNumRows(), a.getNumCols());

  for (i = 0; i < getNumRows(); i++)
  {
    for (j = 0; j < a.getNumCols(); j++)
    {
      v = (T) 0;
      for (k = 0; k < getNumCols(); k++)
      {
        v += val(i,k) * a.val(k,j);
      }
      dm(i,j) = v;
    }
  }

  return dm;
}

template<class T>
inline DenseMatrix<T> DenseMatrix<T>::operator+(const DenseMatrix<T>& a) const
{
  int i,j;

  DenseMatrix<T> dm(getNumRows(), getNumCols());

  for (i = 0; i < getNumRows(); i++)
  {
    for (j = 0; j < getNumCols(); i++)
    {
      dm(i,j) = val(i,j) + b.val(i,j);
    }
  }

  return dm;

}

template<class T>
inline DenseMatrix<T> DenseMatrix<T>::operator-(const DenseMatrix<T>& a) const
{
  int i,j;

  DenseMatrix<T> dm(getNumRows(), getNumCols());

  for (i = 0; i < getNumRows(); i++)
  {
    for (j = 0; j < getNumCols(); i++)
    {
      dm(i,j) = val(i,j) - a.val(i,j);
    }
  }

  return dm;

}

template<class T>
inline DenseMatrix<T> DenseMatrix<T>::transpose() const
{
  int i,j;

  DenseMatrix<T> dm(getNumCols(), getNumRows());

  for (i = 0; i < getNumRows(); i++)
  {
    for (j = 0; j < getNumCols(); j++)
    {
      dm(j,i) = val(i,j);
    }

  }

  return dm;

}

template<class T>
RealVector<T> DenseMatrix<T>::getRow(int row) const
{
  RealVector<T> ret(cols_);
  int i, j;

  for(i=(row*cols_), j=0; i < ((row+1)*cols_); i++, j++)
    ret.set(j, valVec_[i]);

  return ret;
}

template <class T>
RealVector<T> DenseMatrix<T>::getColumn(int col) const
{
  RealVector<T> ret(rows_);
  int i, j;

  for(i=col, j=0; i<rows_*cols_; i+=cols_, j++)
    ret.set(j, valVec_[i]);

  return ret;
}

template <class T>
inline T& DenseMatrix<T>::operator()(int row, int col)
{
  return valVec_[row*cols_ + col];
}

template <class T>
inline T DenseMatrix<T>::operator()(int row, int col) const
{
  return valVec_[row*cols_ + col];
}

} // end namespace
#endif

---