gausselim.h

Go to the documentation of this file.

#ifndef GAUSSELIM_H
#define GAUSSELIM_H

#include <sstream>
using namespace std;

#include <typedefs.h>
#include <zero.h>

template< typename T >
class gausselim
{
  gausselim() { assert(false); }
public:

  T * mat;
  T * row;
  uintc dim;  

  gausselim( T * mat_, uintc dim_ )
    : mat(mat_), dim(dim_)
    { row = new T[dim+1]; }
  ~gausselim()
    { delete[] row; }

  boolc eval();

  void columnC( T * C ) const;


  // Possible Individual Operations on the state.

  void rowget( uintc i );
  void rowput( uintc i );
  void mult(T const alpha);

  void swap( uintc i, uintc k );
  void norm( uintc row, uintc col );

  boolc pivot(uintc k);

  void substitute( uintc i );

  operator stringc () const;
};


//----------------------------------------------------------
//  Implementation


template< typename T >
void gausselim<T>::columnC(T * C) const
{
  for (uint i=0; i<dim; ++i)
    C[i] = mat[3+i*4];
}

template< typename T >
boolc gausselim<T>::eval()
{
  for (uint i=0; i<dim; ++i)
  {
    if (pivot(i)==false)
      return false;

    substitute(i);
  }

  return true;
}

template< typename T >
void gausselim<T>::substitute( uintc k )
{
  norm(k,k);

  rowget(k);
  for (uint i=0; i<dim+1; ++i)
    row[i] *= -1;

  uint b=0;
  T r;

  for (uint i=0; i<dim; ++i)
  {
    if (i==k)
    {
      b += (dim+1);
      continue;
    }

    r = mat[b+k];

    for (uint m=0; m<dim+1; ++m)
      mat[b++] += row[m]*r;
  }

}

template< typename T >
boolc gausselim<T>::pivot( uintc k )
{
  assert(k<dim);

  T lgnum = mat[ k*(dim+1) + k ];
  T lgnum2 = lgnum*lgnum;
  uint lgi = k;
  T x;
  T x2;
  for (uint i=k+1; i<dim; ++i)
  {
    x = mat[ i*(dim+1) + k ];
    x2 = x*x;
    if (x2>lgnum2)
    {
      lgnum2 = x2;
      lgnum = x;
      lgi = i;
    }
  }

  if (k!=lgi)
    swap(k,lgi);

  if (lgnum2 < zero<T>::val )
    return false;

  return true;
}

template< typename T >
void gausselim<T>::rowget( uintc i )
{
  assert(i<dim);

  uint b = i*(dim+1);
  for (uint k=0; k<dim+1; ++k)
    row[k] = mat[b+k];
} 

template< typename T >
void gausselim<T>::rowput( uintc i )
{
  assert(i<dim);

  uint b = i*(dim+1);
  for (uint k=0; k<dim+1; ++k)
    mat[b+k] = row[k];
} 

template< typename T >
void gausselim<T>::mult(T const alpha)
{
  for (uint i=0; i<dim+1; ++i)
    row[i] *= alpha;
}

template< typename T >
void gausselim<T>::norm( uintc row, uintc col )
{
  assert(row<dim);
  assert(col<dim);

  uint b = row*(dim+1);

  T alpha = mat[b+col];
  if (alpha==0)
    return;

  alpha = 1.0/alpha;

  for (uint i=0; i<dim+1; ++i)
    mat[b++] *= alpha;
}

template< typename T >
void gausselim<T>::swap( uintc i, uintc k )
{
  assert(i<dim);
  assert(k<dim);

  T t;
  uint bi = i*(dim+1);
  uint bk = k*(dim+1);
  for (uint j=0; j<dim+1; ++j)
  {
    t = mat[bi];
    mat[bi] = mat[bk];
    mat[bk] = t;
    ++bi;
    ++bk;
  }
}

template< typename T >
gausselim<T>::operator stringc () const
{
  stringstream ss;
  
  for (uint i=0; i<dim; ++i)
  {
    for (uint k=0; k<dim+1; ++k)
    {
      ss << mat[i*(dim+1)+k] << " ";
    }
    if (i!=dim-1)
      ss << "\n";
  }

  return ss.str();
}
 


#endif