Matrix.h

Go to the documentation of this file.

#ifndef Matrix_HEADER
#define Matrix_HEADER

#ifdef HAVE_CONFIG_H
# include <config.h>
#endif

#include <boost/array.hpp>
#include <utilities/foreach.h>
#include "Point.h"

#if 0
#define dprintf(...) printf(__VA_ARGS__);
#else
#define dprintf(...)
#endif

namespace Geometry {

template <unsigned rows, unsigned cols>
class Matrix {
  // Stored in column-major order, for interoperability with LAPACK
  static const unsigned n = rows * cols;
  boost::array<double, n> data_;

  public:
  Matrix() { }
  Matrix(double value) {
    data_.assign(value);
  }
  // Outer product pq^T constructor
  Matrix(const Point<rows>& p, const Point<cols>& q) {
    for(unsigned r = 0; r < rows; ++r) {
      for(unsigned c = 0; c < cols; ++c) {
        get(r, c) = p[r] * q[c];
      }
    }
  }

  const double *column_major_array() const {
    return data_.c_array();
  }
  double *column_major_array() {
    return data_.c_array();
  }

  double operator() (unsigned r, unsigned c) const {
    return data_[r * cols + c];
  }
  double& operator() (unsigned r, unsigned c) {
    return data_[r * cols + c];
  }
  double get(unsigned r, unsigned c) const {
    return data_[r * cols + c];
  }
  double& get(unsigned r, unsigned c) {
    return data_[r * cols + c];
  }

  void operator+= (const Matrix& other) {
    for(unsigned i = 0; i < n; ++i) {
      data_[i] += other.data_[i];
    }
  }
  void operator-= (const Matrix& other) {
    for(unsigned i = 0; i < n; ++i) {
      data_[i] -= other.data_[i];
    }
  }
  void operator*= (double x) {
    for(unsigned i = 0; i < n; ++i) {
      data_[i] *= x;
    }
  }
  void operator/= (double x) {
    for(unsigned i = 0; i < n; ++i) {
      data_[i] /= x;
    }
  }


  Point<rows> operator * (const Point<cols>& p) const {
    Point<rows> q;
    for(unsigned i = 0; i < rows; ++i) {
      q[i] = 0.0;
      for(unsigned j = 0; j < cols; ++j) {
        q[i] += (*this)(i,j) * p[j];
      }
    }
    return q;
  }

  void qd_print()
    {
      for(unsigned i = 0; i < rows; ++i) {
        for(unsigned j = 0; j < cols; ++j) {
          dprintf("%f   ",get(i,j));
        }
        dprintf("\n");
      }
      return;
    }

  /***************************************************************************
    Functions that assume square matrices
   ***************************************************************************/

  template <unsigned newcols>
  Matrix<rows, newcols>
  operator * (const Matrix<cols, newcols>& other) {
    Matrix<rows, newcols> m;
    for(unsigned i = 0; i < newcols; ++i) {
      for(unsigned j = 0; j < rows; ++j) {
        double d = 0.0;
        for(unsigned k = 0; k < cols; ++k) {
          d += get(j, k) * other.get(i, k);
        }
        m(i, j) = d;
      }
    }
    return m;
  }

  private:
  void largest_eigenvector(Matrix<rows, rows> M, Point<rows>& v) const {
    // get the largest eigenvector by repeated squaring
    enum { NSQUARE = 20 };

    for (unsigned i = 0 ; i < NSQUARE; ++i) {
      M = M * M;
      double maxabs = 0.0;
      BOOST_FOREACH(double d, M.data_) {
        d = fabs(d);
        maxabs = maxabs > d ? maxabs : d;
      }
      if (maxabs != 0.0) {
        M /= maxabs;
      } 
    }

    // Look for a non-zero column
    unsigned use_col = 0;
    for (unsigned c = 0; c < cols; ++c) {
      for(unsigned r = 0; r < rows; ++r) {
        if (M(r,c) != 0.0) { 
          use_col = c;
          break;
        }
      }
    }

    // compute the eigenvector, and normalize it
    for (unsigned r = 0; r < rows; ++r) {
      v[r] = M(r,use_col);
    }
    if (v.norm2() != 0.0) {
      v /= sqrt(v.norm2());
    }
  }


  public:
 template<unsigned k> void eigenvectors(boost::array<Point<rows>, k>& vectors) const {
    BOOST_STATIC_ASSERT(rows == cols);

    Matrix<rows, rows> M = *this;
    for (unsigned eigen = 0; eigen < k; ++eigen) {
      largest_eigenvector(M, vectors[eigen]);

      dprintf("Big EIV is %f %f %f\n",vectors[eigen][0],vectors[eigen][1],vectors[eigen][2]);

      // Find the associated eigenvalue by computing the Rayleigh quotient.
      const Point<rows>& v = vectors[eigen];
      double lambda = v.dot(M * v);  // vTv is 1.0 by construction
      assert(lambda != 0); // don't run out of eigenvalues

      dprintf("L is %f\n",lambda);

      // subtract out this eigenvector and move on to the next
      Matrix<rows,rows> vv(v, v);
      vv *= lambda;
      M -= vv;

      dprintf("Lvv then M-Lvv\n");
      vv.qd_print();
      dprintf("----\n");
      M.qd_print();
    }
  }

};

}

#undef dprintf
#endif

---