linalg.h

Go to the documentation of this file.

 
#ifndef LINALG_H
#define LINALG_H
 
#include <vector>
using namespace std;
 
#include "mycomplex.h"
#include "xpr-vector.h"
#include "xpr-matrix.h"
#include "def.h"
#include "utils.h"
#include <iosfwd>
#include <valarray>
#include <sstream>
using namespace softsusy;
 
/************************************
 * 
 *   *** DOUBLE VECTOR CLASS ***
 *
 ************************************/
 
class DoubleVector : public Indexable<double, DoubleVector> {
private:
  int start, end; 
  std::valarray<double> x; 
  DoubleVector(const std::valarray<double> & valarr, int s, int e)
    : Indexable<double,DoubleVector>(), start(s), end(e), x(valarr) {}
  friend class DoubleMatrix;
 
public:
 
  /*
   *  CONSTRUCTORS
   */
 
 
  explicit DoubleVector(int e) 
    : Indexable<double,DoubleVector>(), start(1), end(e), x(0.0,e) {}
  DoubleVector(int s, int e) 
    : Indexable<double,DoubleVector>(), start(s), end(e), x(0.0,e-s+1) {}
  DoubleVector(double a, double b)
    : Indexable<double,DoubleVector>(), start(1), end(2), x(0.0,2) {
    x[0] = a; x[1] = b; }
  ~DoubleVector() {}
 
  template<class E>
  DoubleVector & operator=(const Xpr<double,E> & x) {
    *this = copy_from(x);
    return *this;
  }
 
  DoubleVector& operator=(const DoubleVector& other) {
    x.resize(other.x.size());
    x = other.x;
    start = other.start;
    end = other.end;
    return *this;
  }
 
  template <typename E>
  DoubleVector(const Xpr<double,E> & v)
    : Indexable<double,DoubleVector>(), 
      start(v.displayStart()), end(v.displayEnd()), 
      x(0.0,v.displayEnd()-v.displayStart()+1)
  {
    assign_from(v);
  }
 
  /*
   *  ELEMENT ACCESS
   */
 
  double & operator() (int i); 
  double operator() (int i) const { return x[i-start]; }
  double display(int i) const { return x[i-start]; }
  void set(int i, double f);
 
  int displayStart() const { return start; }
  int displayEnd() const { return end; }
  const DoubleVector & display() const { return *this; }
  double compare(const DoubleVector & a) const; 
 
  void append(double);
 
  void setEnd(int e);
  DoubleVector abs() const { return DoubleVector(std::abs(x), start, end); }
 
  double norm() const;
  DoubleVector apply(double (*fn)(double)) const { 
    return DoubleVector(x.apply(fn),start,end); 
  }
  
  double nmin(int & p) const;
  double max() const { return x.max(); }
 
  double max(int & p) const {
    double m = double(0);
    int i; 
    for (i=displayStart(); i<=displayEnd();  i++)
      if (display(i) > m) {
        m = display(i);
        p = i;
      }
    return m;
  }
 
  double min(int & p) const;
 
  double absmax(int & p) const { return abs().max(p); }
  double sumElements() const { return std::abs(x).sum(); }  
  std::size_t size() const { return x.size(); }
  
  void swap(int i, int j); 
  
  double dot(const DoubleVector & v) const { 
    return (x * v.x).sum();
  }
 
  DoubleVector sort() const {
    DoubleVector temp(*this), saveIt(*this);
    
    int j, pos;
    for (j=start; j<=end; j++) {
      double d = temp.min(pos);
      if (d < 1.0e66) saveIt(j) = d;
      
      temp(pos) = 6.66e66; // take that element out of the min calculation
    }
    return saveIt;
  }
  
  DoubleVector divide(DoubleVector const &v) const;
  double average() const; 
  int closest(double a) const; 
  void fillArray(double* array, unsigned offset = 0) const;
};
 
 
/*
 *  STANDARD INPUT / OUTPUT
 */
 
std::ostream & operator<<(std::ostream &left, const DoubleVector &V);
std::istream & operator>>(std::istream & left, DoubleVector &V);
 
 
/*
 *  INLINE FUNCTION DEFINITIONS
 */
 
inline double & DoubleVector::operator() (int i) {
#ifdef ARRAY_BOUNDS_CHECKING
  if (i>end || i<start) {
      std::ostringstream ii;
      ii << "Trying to access " << i << "th element of DoubleVector\n";
      ii << "start " << start << " end " << end;
      ii << *this;
      throw ii.str();
  }
#endif
  return x[i-start];
}
 
inline void DoubleVector::set(int i, double f) {
#ifdef ARRAY_BOUNDS_CHECKING
  if (i < start || i > end) 
    {
      std::ostringstream ii;
      ii << "Cannot access " << i << "th element of vector " << *this;
      throw ii.str();
    }
#endif
  x[i-start] = f;
}
 
 
/************************************
 * 
 *   *** DOUBLE MATRIX CLASS ***
 *
 ************************************/
 
 
class ComplexVector; class ComplexMatrix;
class DoubleMatrix : public MatIndexable<double, DoubleMatrix> {
private:
  int rows, cols;
  std::valarray<double> x; 
  
  // for usage see e.g. http://www.roguewave.com/support/docs/sourcepro/stdlibref/2-10.html
  std::slice_array<double>  row(int i) { return x[std::slice((i-1)*cols,cols,1)]; }
  std::slice_array<double>  col(int i) { return x[std::slice(i-1,rows,cols)]; }
  const std::valarray<double>  row(int i) const {
    return std::valarray<double>(x[std::slice((i-1)*cols,cols,1)]); 
  }
  const std::valarray<double>  col(int i) const {
    return std::valarray<double>(x[std::slice(i-1,rows,cols)]);
  }
 
  std::slice_array<double>  diag() { return x[std::slice(0,rows,cols+1)]; }
  const std::valarray<double>  diag() const {
    return std::valarray<double>(x[std::slice(0,rows,cols+1)]);
  }
 
  double & elmt(int i, int j) { return x[(i-1) * cols + (j-1)]; }
  double elmt(int i, int j) const { return x[(i-1) * cols + (j-1)]; }
 
  DoubleMatrix(const std::valarray<double> & valarr, int r, int c)
    : MatIndexable<double,DoubleMatrix>(), rows(r), cols(c), x(valarr) {}
 
  friend class ComplexMatrix;
  
public:
 
  /*
   *  CONSTRUCTORS
   */
 
  DoubleMatrix(int r, int c) 
    : MatIndexable<double,DoubleMatrix>(), rows(r), cols(c), x(0.0, r*c) {}
  explicit DoubleMatrix(const DoubleVector &v);  
  ~DoubleMatrix() {}
 
  template<class E>
  DoubleMatrix & operator=(const MatXpr<double,E> & x) {
    *this = copy_from(x);
    return *this;
  }
 
  DoubleMatrix& operator=(const DoubleMatrix& other) {
    x.resize(other.x.size());
    x = other.x;
    rows = other.rows;
    cols = other.cols;
    return *this;
  }
 
  template <typename E>
  DoubleMatrix(const MatXpr<double, E> & m)
    : MatIndexable<double,DoubleMatrix>(), 
      rows(m.displayRows()), cols(m.displayCols()), x(0.0, rows*cols) 
  {
    assign_from(m);
  }
 
  /*
   *  ELEMENT ACCESS
   */
  double operator() (int i, int j) const; 
  double & operator() (int i, int j); 
  double display(int i, int j) const;
  double sumElements() const { return abs(x).sum(); }  
  int displayRows() const { return rows; };
  int displayCols() const { return cols; };
  const DoubleMatrix & display() const { return *this; }
 
  DoubleVector displayRow(int i) const {
    DoubleVector temp(displayCols());
    int j; for (j=1; j<=displayCols(); j++) temp(j) = display(i, j);
    return temp;
  }
 
  DoubleVector displayCol(int i) const {
    DoubleVector temp(displayRows());
    int j; for (j=1; j<=displayRows(); j++) temp(j) = display(j, i);
    return temp;
  }
 
 
  void operator+=(DoubleMatrix & f) {
    x += f.x;
  }
 
  double nmin(int & k, int & l) const;
 
  DoubleMatrix apply(double (*fn)(double)) const { 
    return DoubleMatrix(x.apply(fn),rows,cols); 
  }
 
  const DoubleMatrix & operator=(double v);  
 
  double min(int & k, int & l) const; 
  double max(int & k, int & l) const; 
  void swaprows(int i, int j);
  void swapcols(int i,int j);
  std::size_t size() const { return x.size(); }
  void setCols(int);
  void setRows(int);
  void resize(int, int);
 
  double trace() const;
  DoubleMatrix transpose() const; 
  bool testNan() const; 
 
  /*
   *  NUMERICAL DIAGONALIZATION ROUTINES ETC.
   */
 
  void associateOrderAbs(DoubleVector &v);
  // Perform on asymmetric matrices M such that
  void associateOrderAbs(DoubleMatrix & u, DoubleMatrix & v, DoubleVector &w) 
    const;
  void symmetrise(); 
  double compare(const DoubleMatrix & a) const;
  double compare(const ComplexMatrix & a) const;
  DoubleVector diagVals() const;
  DoubleMatrix inverse() const; 
 
  double diagonalise(DoubleMatrix & u, DoubleMatrix & v, DoubleVector & w) const;
  double diagonaliseSym(DoubleMatrix & v, DoubleVector & w) const;
  double diagonaliseSym(ComplexMatrix & v, DoubleVector & w) const;
 
  // -[ cos theta    sin theta ]   A   [ cos theta -sin theta ] = diagonal
  // -[ -sin theta   cos theta ]       [ sin theta  cos theta ]
  DoubleVector sym2by2(double & theta) const; 
  // -[ cos thetaL    sin thetaL ]   A   [ cos thetaR -sin thetaR ] = diagonal
  // -[ -sin thetaL   cos thetaL ]       [ sin thetaR  cos thetaR ]
  DoubleVector asy2by2(double & thetaL, double & thetaR) const;
  DoubleVector vectorfy() const;
  double determinant() const; 
 
  DoubleMatrix ludcmp(double & d) const;
  void fillArray(double* array, unsigned offset = 0) const;
};
 
/*
 *  STANDARD INPUT / OUTPUT
 */
 
std::istream & operator>>(std::istream & left, DoubleMatrix &M);
std::ostream & operator <<(std::ostream &left, const DoubleMatrix &V);
 
/*
 *  NON-MEMBER DIAGONALIZATION ROUTINES ETC.
 */
 
// [ cos theta    sin theta ]
// [-sin theta    cos theta ]
DoubleMatrix rot2d(double theta);
// [ -sin(theta)  cos(theta) ]
// [  cos(theta)  sin(theta) ] --
DoubleMatrix rot2dTwist(double theta);
 
// [ -cos theta     sin theta            0 ]
// [ sin theta                  cos theta                0 ]
// [ 0                                                  0                                                1 ]
DoubleMatrix rot3d(double theta);
 
// [ cos thetaL    sin thetaL ]   A   [ cos thetaR -sin thetaR ]  = diag
// [ -sin thetaL   cos thetaL ]       [ sin thetaR  cos thetaR ]
// as given by asy2by2!
void positivise(double thetaL, double thetaR, const DoubleVector & diag,
                  ComplexMatrix & u, ComplexMatrix & v);
void diagonaliseSvd(DoubleMatrix & a, DoubleVector & w, DoubleMatrix & v);
double pythagoras(double a, double b);
void diagonaliseJac(DoubleMatrix & a,  int n,  DoubleVector & d, DoubleMatrix
                    & v,  int *nrot);
 
void fillArray(const std::valarray<double>&, double*, unsigned offset = 0);
 
/*
 *  INLINE FUNCTION DEFINITIONS
 */
 
inline double DoubleMatrix::operator() (int i, int j) const { 
#ifdef ARRAY_BOUNDS_CHECKING
  if (i > rows || j > cols || i < 1 || j < 1) { 
    std::ostringstream ii;
    ii << "Trying to access " << i << "," << j << 
      "th element of DoubleMatrix\n" << *this;
    throw ii.str();
  }
#endif
  return elmt(i,j); 
}  
 
inline double &  DoubleMatrix::operator() (int i, int j) { 
#ifdef ARRAY_BOUNDS_CHECKING
  if (i > rows || j > cols || i < 1 || j < 1) { 
    std::ostringstream ii;
    ii << "Trying to access " << i << "," << j << 
      "th element of DoubleMatrix\n" << *this;
    throw ii.str();
  }
#endif
  return elmt(i,j); 
}  
 
inline double DoubleMatrix::display(int i, int j) const { 
#ifdef ARRAY_BOUNDS_CHECKING
  if (i < 1 || i > displayRows() || j < 1 || j > displayCols()) { 
    std::ostringstream ii;
    ii << "Error: Requested display (" << i << "," << j << 
      ")th element of matrix " << *this; 
    throw ii.str();
  }
#endif
  return elmt(i,j); 
}  
 
/************************************
 * 
 *   *** COMPLEX VECTOR CLASS ***
 *
 ************************************/
 
class ComplexVector : public Indexable<Complex, ComplexVector> {
private:
  int start, end; 
  std::valarray<Complex> x; 
  ComplexVector(const std::valarray<Complex> & valarr, int s, int e)
    : Indexable<Complex, ComplexVector>(), start(s), end(e), x(valarr) {}  
  friend class ComplexMatrix;
 
 
public:
 
  /*
   *  CONSTRUCTORS
   */
 
  explicit ComplexVector(int e)
    : Indexable<Complex, ComplexVector>(), start(1), end(e), x(0.0,e) {}
  ComplexVector(int s, int e)
    : Indexable<Complex, ComplexVector>(), start(s), end(e), x(0.0,e-s+1) {}
  ~ComplexVector() {}
 
  template<class E>
  ComplexVector & operator=(const Xpr<Complex,E> & x) {
    *this = copy_from(x);
    return *this;
  }
 
  ComplexVector& operator=(const ComplexVector& other) {
    x.resize(other.x.size());
    x = other.x;
    start = other.start;
    end = other.end;
    return *this;
  }
 
  template <typename E>
  ComplexVector(const Xpr<Complex,E> & v)
    : Indexable<Complex,ComplexVector>(), 
      start(v.displayStart()), end(v.displayEnd()), 
      x(0.0,v.displayEnd()-v.displayStart()+1)
  {
    assign_from(v);
  }
 
 
  /*
   *  ELEMENT ACCESS
   */
 
  Complex & operator() (int i);
  Complex operator() (int i) const { return x[i-start]; }
  Complex display(int i) const { return x[i-start]; }
  void set(int i, Complex f);
 
  int displayStart() const { return start; }
  int displayEnd() const { return end; }
  const ComplexVector & display() const { return *this; }
 
 
 
  void setEnd(int e);
 
  ComplexVector apply(Complex (*fn)(Complex)) const { 
    return ComplexVector(x.apply(fn),start,end);
  }
  
  Complex max() const { return x.max(); }
  Complex min(int & p) const;
  
  void swap(int i, int j); 
 
};
 
/*
 *  STANDARD INPUT / OUTPUT
 */
 
std::ostream & operator<<(std::ostream &left, const ComplexVector &V);
std::istream & operator>>(std::istream & left, ComplexVector &V);
 
/*
 *  INLINE FUNCTION DEFINITIONS
 */
 
inline Complex & ComplexVector::operator() (int i) {
#ifdef ARRAY_BOUNDS_CHECKING
  if (i>end || i<start) {
    std::ostringstream ii;
    ii << "Trying to access " << i << "th element of DoubleVector\n";
    ii << "start " << start << " end " << end;
    ii << *this;
    throw ii.str();
  }
#endif
  return x[i-start];
}
 
inline void ComplexVector::set(int i, Complex f) {
#ifdef ARRAY_BOUNDS_CHECKING
  if (i < start || i > end) {
    std::ostringstream ii;
    ii << "Cannot access " << i << "th element of vector " << *this;
    throw ii.str();
  }
#endif
  x[i-start] = f;
}
 
/************************************
 * 
 *   *** COMPLEX MATRIX CLASS ***
 *
 ************************************/
 
class ComplexMatrix : public MatIndexable<Complex, ComplexMatrix> {
private:
  int rows, cols;
  std::valarray<Complex> x;
  
  // for usage see e.g. http://www.roguewave.com/support/docs/sourcepro/stdlibref/2-10.html
  std::slice_array<Complex>  row(int i) { return x[std::slice((i-1)*cols,cols,1)]; }
  std::slice_array<Complex>  col(int i) { return x[std::slice(i-1,rows,cols)]; }
  const std::valarray<Complex>  row(int i) const {
    return std::valarray<Complex>(x[std::slice((i-1)*cols,cols,1)]); 
  }
  const std::valarray<Complex>  col(int i) const {
    return std::valarray<Complex>(x[std::slice(i-1,rows,cols)]);
  }
 
  std::slice_array<Complex>  diag() { return x[std::slice(0,rows,cols+1)]; }
  const std::valarray<Complex>  diag() const {
    return std::valarray<Complex>(x[std::slice(0,rows,cols+1)]);
  }
 
  Complex & elmt(int i, int j) { return x[(i-1)*cols+(j-1)]; }
  const Complex elmt(int i, int j) const { return x[(i-1)*cols+(j-1)]; }
 
  ComplexMatrix(const std::valarray<Complex> & valarr, int r, int c)
    : MatIndexable<Complex,ComplexMatrix>(), rows(r), cols(c), x(valarr) {}
 
  friend class DoubleMatrix;
  
public:
 
  /*
   *  CONSTRUCTORS
   */
 
  ComplexMatrix(int r, int c) 
    : MatIndexable<Complex,ComplexMatrix>(), rows(r), cols(c), x(0.0, r*c) {}
  explicit ComplexMatrix(const DoubleMatrix & m);
  explicit ComplexMatrix(const ComplexVector &v);  
  ~ComplexMatrix() {}
 
  template<class E>
  ComplexMatrix & operator=(const MatXpr<Complex,E> & x) {
    *this = copy_from(x);
    return *this;
  }
 
  ComplexMatrix& operator=(const ComplexMatrix& other) {
    x.resize(other.x.size());
    x = other.x;
    rows = other.rows;
    cols = other.cols;
    return *this;
  }
 
  template <typename E>
  ComplexMatrix(const MatXpr<Complex, E> & m)
    : MatIndexable<Complex,ComplexMatrix>(), 
      rows(m.displayRows()), cols(m.displayCols()), x(0.0, rows*cols) 
  {
    assign_from(m);
  }
  
  /*
   *  ELEMENT ACCESS
   */
 
  Complex & operator() (int i, int j); 
  Complex operator() (int i, int j) const;
  Complex display(int i, int j) const;
 
  int displayRows() const { return rows; }; 
  int displayCols() const { return cols; };
  const ComplexMatrix & display() const { return *this; };
  std::size_t size() const { return x.size(); }
 
  const ComplexMatrix & operator=(const Complex &v);  
 
  Complex min(int & k, int & l) const; 
  void swaprows(int i, int j);
  void swapcols(int i,int j);
 
  void setCols(int);
  void setRows(int);
  void resize(int, int);
 
  Complex trace() const;  
  ComplexMatrix transpose() const;
  ComplexMatrix hermitianConjugate() const;
  ComplexMatrix complexConjugate() const;
  DoubleMatrix real() const;
  DoubleMatrix imag() const;
  
  /*
   *  NUMERICAL DIAGONALIZATION ROUTINES ETC.
   */ 
  double nonHermiticity() const;
  void symmetrise();
  double compare(const ComplexMatrix & a) const;
  DoubleMatrix makeHermitianRealForDiag() const;
  double diagonaliseHerm(ComplexMatrix & v, DoubleVector & w) const;
  double takagi(ComplexMatrix & v, ComplexVector & w) const;
  double diagonaliseSym2by2(ComplexMatrix & v, ComplexVector & w) const;
  double diagonalise(ComplexMatrix & u, ComplexMatrix & v, DoubleVector & w) 
    const;
  ComplexMatrix apply(Complex (*fn)(Complex)) const { 
    ComplexMatrix a(rows, cols);
    for (int i=1; i<=rows; i++)
      for (int j=1; j<=rows; j++) 
        a(i, j) = fn(display(i, j));
      
    return a;
  }
};
 
/*
 *  STANDARD INPUT / OUTPUT
 */
 
std::istream & operator>>(std::istream & left, ComplexMatrix &M);
std::ostream & operator<<(std::ostream &left, const ComplexMatrix &V);
 
/*
 *  INLINE FUNCTION DEFINITIONS
 */
 
inline Complex & ComplexMatrix::operator() (int i, int j) { 
#ifdef ARRAY_BOUNDS_CHECKING
  if (i > rows || j > cols || i < 0 || j < 0) { 
    std::ostringstream ii;
    ii << "Trying to access " << i << "," << j << 
      "th element of ComplexMatrix\n" << *this;
    throw ii.str();
  }
#endif
  return elmt(i,j); 
}  
 
inline Complex ComplexMatrix::operator() (int i, int j) const { 
#ifdef ARRAY_BOUNDS_CHECKING
  if (i > rows || j > cols || i < 0 || j < 0) { 
    std::ostringstream ii;
    ii << "Trying to access " << i << "," << j << 
      "th element of ComplexMatrix\n" << *this;
    throw ii.str();
  }
#endif
  return elmt(i,j); 
}  
 
inline Complex ComplexMatrix::display(int i, int j) const { 
#ifdef ARRAY_BOUNDS_CHECKING
  if (i < 1 || i > displayRows() || j < 1 || j > displayCols()) { 
    std::ostringstream ii;
    ii << "Error: Requested display (" << i << "," << j << 
        ")th element of matrix " << *this; 
    throw ii.str();
    }
#endif
  return elmt(i,j); 
}  
 
 
#endif