#ifndef FILE_BANDMATRIX
#define FILE_BANDMATRIX

/****************************************************************************/
/* File:   bandmatrix.hpp                                                   */
/* Author: Joachim Schoeberl                                                */
/* Date:   14. Aug. 2002                                                    */
/****************************************************************************/

/// A symmetric band-matrix
template <class T>
class FlatSymBandMatrix
{
protected:
  ///
  int n;
  ///
  int bw;
  ///
  T *data;
public:
  typedef typename mat_traits<T>::TV_COL TV;

  ///
  FlatSymBandMatrix (int an, int abw, T * adata)
    : n(an), bw(abw), data(adata)
  { ; }

  ///
  void Mult (const FlatVector<TV> & x, FlatVector<TV> & y) const
  {
    for (int i = 0; i < n; i++)
      y(i) = (*this)(i,i) * x(i);
    for (int i = 0; i < n; i++)
      for (int j = max2(i-bw+1, 0); j < i; j++)
	{
	  y(i) += (*this)(i,j) * x(j);
	  y(j) += Trans((*this)(i,j)) * x(i);
	}
  }

  ///
  ostream & Print (ostream & ost) const
  {
    for (int i = 0; i < n; i++)
      {
	for (int j = 0; j < n; j++)
	  if (Used (i, j))
	    ost << setw(8) << (*this)(i,j) << " ";
	  else if (Used (j,i))
	    ost << setw(8) << "sym" << " ";
	  else
	    ost << setw(8) << 0;
	ost << endl;
      }
    return *this;
  }

  ///
  int Height() const { return n; }

  ///
  int BandWidth() const { return bw; }
  
  ///
  const T & operator() (int i, int j) const
  { return data[i * bw + j - i + bw-1]; }

  ///
  T & operator() (int i, int j) 
  { return data[i * bw + j - i + bw-1]; }

  ///
  bool Used (int i, int j) const
  {
    return (n > i && i >= j && j >= 0 && i-j < bw);
  }
  
  ///
  FlatSymBandMatrix & operator= (const T & val)
  {
    for (int i = 0; i < bw * n; i++)
      data[i] = val;
    return *this;
  }
};


template<typename T>
inline std::ostream & operator<< (std::ostream & s, const FlatSymBandMatrix<T> & m)
{
  m.Print (s);
  return s;
}



///  A symmetric band-matrix with memory management
template <class T>
class SymBandMatrix : public FlatSymBandMatrix<T>
{
public:
  typedef typename mat_traits<T>::TV_COL TV;

  ///
  SymBandMatrix (int an, int abw)
    : FlatSymBandMatrix<T> (an, abw, new T[an*abw])
  { ; }

  ///
  ~SymBandMatrix ()
  { delete [] this->data; }

  ///
  SymBandMatrix & operator= (const T & val)
  {
    for (int i = 0; i < this->bw * this->n; i++)
      this->data[i] = val;
    return *this;
  }
};












/*
  lfact (bw = 3)

  d0 
  0   d1
  1   2   d2
      3   4   d3
*/

/**
   Object for Cholesky factors of a band matrix.
   This class does not provide memory management.
 */
template <class T>
class FlatBandCholeskyFactors
{
protected:
  ///
  int n;
  ///
  int bw;
  ///
  T * mem; // first diag, than lfact
public:
  typedef typename mat_traits<T>::TV_COL TV;

  /// assigne dimension, bandwith and memory
  FlatBandCholeskyFactors (int an, int abw, T * amem)
  { n = an, bw = abw, mem = amem; }

  ///
  FlatBandCholeskyFactors ()
  { n = bw = 0; mem = 0; }

  /// factor bandmatrix a
  void Factor (const FlatSymBandMatrix<T> & a);

  /// solve with factored matrices
  void Mult (const FlatVector<TV> & x, FlatVector<TV> & y) const;

  ///
  ostream & Print (ostream & ost) const;

  /// compute linear position of matrix element (i,j)
  int Index (int i, int j) const
  {
    if (i < bw)
      return n + (i * (i-1)) / 2 + j;
    else
      return n + i * (bw-2) + j - ((bw-1)*(bw-2))/2;   
  }
  
  /// matrix element (i,j), (i,j) must be a valid position
  const T & operator() (int i, int j) const
  { 
    if (i < bw)
      return mem[n + (i * (i-1)) / 2 + j];
    else
      return mem[n + i * (bw-2) + j - ((bw-1)*(bw-2))/2];   
  }

  /// matrix element (i,j), (i,j) must be a valid position
  T & operator() (int i, int j) 
  {
    if (i < bw)
      return mem[n + (i * (i-1)) / 2 + j];
    else
      return mem[n + i * (bw-2) + j - ((bw-1)*(bw-2))/2];   
  }

  ///
  int Size() const { return n; }
  ///
  int BandWidth() const { return bw; }
  ///
  static int RequiredMem (int n, int bw)
  { return n*bw - (bw * (bw-1)) / 2 + n; }
};
  

///  
template<typename T>
inline std::ostream & operator<< (std::ostream & s, const FlatBandCholeskyFactors<T> & m)
{
  m.Print (s);
  return s;
}





/**
   Cholesky factors of a band matrix.
 */
template <class T>
class BandCholeskyFactors : public FlatBandCholeskyFactors<T>
{
public:
  /// allocate memory and factor a
  BandCholeskyFactors (const SymBandMatrix<T> & a)
    : FlatBandCholeskyFactors<T> (a.Height(), 
				  a.BandWidth(),
				  new T[RequiredMem (a.Height(), a.BandWidth())])
  {
    Factor (a);
  }

  ///
  ~BandCholeskyFactors ()
  {
    delete [] this->mem;
  }
};


#endif