Logo Search packages:      
Sourcecode: gambit version File versions  Download package

cmatrix.h

/* Copyright 2002 Ben Blum, Christian Shelton
 * 
 * This file is part of GameTracer.
 *
 * GameTracer is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * GameTracer is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with GameTracer; if not, write to the Free Software Foundation, 
 * Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#ifndef _CMATRIX_H_
#define _CMATRIX_H_

#include <math.h>
#ifdef SOLARIS
#include <ieeefp.h>
#endif
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <assert.h>
#include <strings.h>
#include <string>
#include <iomanip>

using namespace std;
class cmatrix;

class cmatrixrow {
public:
      double operator[](int) const {
            return 0.0;
      }
};

class cvector {
friend class cmatrix;
public:
 static int num_vec_cons; 
      inline cvector() {
        cvector::num_vec_cons++;
            m = 1;
            x = new double[1];
      }
      inline cvector(int m) {
        cvector::num_vec_cons++;
            this->m = m;
            x = new double[m];
      }
      ~cvector(); 
      inline cvector(const cvector &v) {
        cvector::num_vec_cons++;
            m = v.m;
            x = new double[m];
            //for(int i=0;i<m;i++) x[i] = v.x[i];
            memcpy(x,v.x,m*sizeof(double));
      }
      inline cvector(int m, const double &a) {
        cvector::num_vec_cons++;
            this->m = m;
            x = new double[m];
            for(int i=0;i<m;i++) x[i] = a;
      }
      inline cvector(double *v, int m, bool keep=false) {
        cvector::num_vec_cons++;
            this->m = m;
            if (keep) x = v;
            else {
                  x = new double[m];
                  //for(int i=0;i<m;i++) x[i] = v[i];
                  memcpy(x,v,m*sizeof(double));
            }
      }
      inline cvector operator-() const {
            cvector ret(m);
            for(int i=0;i<m;i++) ret.x[i] = -x[i];
            return ret;
      }
      inline cvector& operator=(double a) {
            for(int i=0;i<m;i++) x[i] = a;
            return *this;
      }
      inline cvector& operator=(const cvector &v) {
            if (&v==this) return *this;
            if (v.m != m) {
                  delete []x;
                  m = v.m;
                  x = new double[m];
            }
            //for(int i=0;i<m;i++) x[i] = v.x[i];
            memcpy(x,v.x,m*sizeof(double));
            return *this;
      }
      inline bool isvalid() const {
            for(int i=0;i<m;i++) if (!finite(x[i])) return false;
            return true;
      }
      inline double operator*(const cvector &v) const {
            if (m!=v.m) {
                  cerr << "invalid cvector dot product" << endl;
                  assert(0);
            }
            double ret = 0.0;
            for(int i=0;i<m;i++) ret += x[i]*v.x[i];
            return ret;
      }
      inline double operator*(const double *v) const {
            double ret = 0.0;
            for(int i=0;i<m;i++) ret += x[i]*v[i];
            return ret;
      }
      inline cvector outer(const cvector &v) const {
            cvector ret(m*v.m);
            for(int i=0,c=0;i<m;i++)
                  for(int j=0;j<v.m;j++,c++)
                        ret.x[c] = x[i]*v.x[j];
            return ret;
      }
      inline double operator[](int i) const {
            return x[i];
      }
      inline double &operator[](int i) {
            return x[i];
      }
      inline cvector &operator+=(const cvector &v) {
            if (v.m!=m) {
                  cerr << "invalid cvector addition" << endl;
                  assert(0);
            }
            for(int i=0;i<m;i++) x[i] += v.x[i];
            return *this;
      }
      inline cvector &operator-=(const cvector &v) {
            if (v.m!=m) {
                  cerr << "invalid cvector subtraction" << endl;
                  assert(0);
            }
            for(int i=0;i<m;i++) x[i] -= v.x[i];
            return *this;
      }
      inline cvector &operator*=(const double &a) {
            for(int i=0;i<m;i++) x[i] *= a;
            return *this;
      }
      inline cvector &operator+=(const double &a) {
            for(int i=0;i<m;i++) x[i] += a;
            return *this;
      }
      inline cvector &operator-=(const double &a) {
            for(int i=0;i<m;i++) x[i] -= a;
            return *this;
      }
      inline cvector &operator/=(const double &a) {
            for(int i=0;i<m;i++) x[i] /= a;
            return *this;
      }
      inline double max() const {
            double t,ma = x[0];
            for(int i=1;i<m;i++)
                  if((t=x[i])>ma) ma = t;
            return ma;
      }
      inline double min() const {
            double t,mi = x[0];
            for(int i=1;i<m;i++)
                  if ((t=x[i])<mi) mi = t;
            return mi;
      }
      inline double absmax() const {
            double t,ma = x[0]>0?x[0]:-x[0];
            for(int i=1;i<m;i++)
                  if ((t=(x[i]>0?x[i]:-x[i]))>ma) ma = t;
            return ma;
      }
      inline double absmin() const {
            double t,mi = x[0]>0?x[0]:-x[0];
            for(int i=1;i<m;i++)
                  if ((t=(x[i]>0?x[i]:-x[i]))<mi) mi = t;
            return mi;
      }
      inline double normalize() {
            double norm = 0.0;
            for(int i = 0; i < m; i++)
                  norm += x[i] * x[i];
            norm = sqrt(norm);
            for(int i = 0; i < m; i++)
                  x[i] /= norm;
            return norm;
      }
      inline bool operator==(const cvector &v) const {
            if (m!=v.m) return false;
#if !defined(HAVE_BCMP)
            for(int i=0;i<m;i++) if (v.x[i]!=x[i]) return false;
            return true;
#else
            return bcmp(x,v.x,m*sizeof(double))==0;
#endif // HAVE_BCMP
      }
      inline bool IsEqual(cvector *v) const {
            if (m!=v->m) return false;
#if !defined(HAVE_BCMP)
            for(int i=0;i<m;i++) if (v->x[i]!=x[i]) return false;
            return true;
#else
            return bcmp(x,v->x,m*sizeof(double))==0;
#endif // HAVE_BCMP
      }     
      inline bool operator==(const double &a) const {
            for(int i=0;i<m;i++) if (a!=x[i]) return false;
            return true;
      }
      inline bool operator!=(const cvector &v) const {
            if (m!=v.m) return true;
#if !defined(HAVE_BCMP)
            for(int i=0;i<m;i++) if (v.x[i]!=x[i]) return true;
            return false;
#else
            return bcmp(x,v.x,m*sizeof(double))!=0;
#endif // HAVE_BCMP
      }
      inline bool operator!=(const double &a) const {
            for(int i=0;i<m;i++) if (a!=x[i]) return true;
            return false;
      }
      inline double norm2() const {
            double ret=x[0]*x[0];
            for(int i=1;i<m;i++) ret += x[i]*x[i];
            return ret;
      }
      inline double norm() const {
            return sqrt(norm2());
      }
      friend ostream& operator<<(ostream &s, const cvector &v);
      friend istream& operator>>(istream &s, cvector &v);

      inline double *values() {
            return x;
      }
      
      inline int getm() const { return m; }

      inline ostream &niceprint(ostream &s) {
            for(int i=0;i<m;i++)
                  s << x[i] << ' ';
            s << endl;
            return s;
      }

      inline void unfuzz(double fuzz) {
        for(int i=0; i < m; i++)
          if(x[i] < fuzz) x[i] = 0.0;
      }

      inline double sum() {
        double total = 0.0;
        for(int i = 0; i < m; i++)
          total += x[i];
        return total;
      }

      inline void support(int *s) {
        for(int i = 0; i < m; i++)
          if(!s[i])
            x[i] = 0.0;
      }

      inline void negate() {
        for(int i = 0; i < m; i++) {
          x[i] = -x[i];
        }
      }

private:
      int m;
      double *x;
};

inline double max(double f1, double f2) {
        return ((f1 > f2) ? f1 : f2);
}
inline cvector operator+(const cvector &a, const cvector &b) {
      return cvector(a)+=b;
}
inline cvector operator-(const cvector &a, const cvector &b) {
      return cvector(a)-=b;
}
inline cvector operator+(const cvector &a, const double &b) {
      return cvector(a)+=b;
}
inline cvector operator-(const cvector &a, const double &b) {
      return cvector(a)-=b;
}
inline cvector operator+(const double &a, const cvector &b) {
      return cvector(b)+=a;
}
inline cvector operator-(const double &a, const cvector &b) {
      return cvector(b.getm(),a)-=b;
}
inline cvector operator*(const cvector &a, const double &b) {
      return cvector(a)*=b;
}
inline cvector operator*(const double &a, const cvector &b) {
      return cvector(b)*=a;
}
inline cvector operator/(const cvector &a, const double &b) {
      return cvector(a)/=b;
}

inline ostream &operator<<(ostream &s, const cvector& v) {
//    s << v.m << ' ';
      for(int i=0;i<v.m;i++) { s << v.x[i]; if (i!=v.m) s << ' '; }
      return s;
}

inline istream &operator>>(istream &s, cvector& v) {
      int tm;
      s >> tm;
      if (tm!=v.m) {
            delete []v.x;
            v.m = tm;
            v.x = new double[tm];
      }
      for(int i=0;i<tm;i++) s >> v.x[i];
      return s;
}

class cmatrix {
public:
      inline cmatrix(int m=1, int n=1) {
            this->m = m; this->n = n;
            s = m*n;
            x = new double[s];
      }
      ~cmatrix();
      inline cmatrix(const cmatrix &ma, bool transpose=false) {
            s = ma.m*ma.n;
            x = new double[s];
            if (transpose) {
                  int i,j,c;
                  n = ma.m; m = ma.n;
                  c = 0;
                  for(i=0;i<m;i++) for(j=0;j<n;j++,c++)
                        x[c] = ma.x[i+j*m];
            } else {
                  n = ma.n; m = ma.m;
                  int i;
                  for(i=0;i<s;i++) x[i] = ma.x[i];
            }
      }
      inline cmatrix(int m, int n,const double &a,bool diaonly=false) {
            this->m = m;
            this->n = n;
            s = m*n;
            x = new double[s];
            if (diaonly) {
                  int i;
                  //for(i=0;i<s;i++) x[i] = 0;
                  memset(x,0,s*sizeof(double));
                  if (n>=m)
                        for(i=0;i<m;i++) x[i*n+i] = a;
                  else for(i=0;i<n;i++) x[i*n+i] = a;
            } else {
                  int i;
                  if (a==0.0) memset(x,0,s*sizeof(double));
                  else for(i=0;i<s;i++) x[i] = a;
            }
      }
      // put v on the diagonal
      inline cmatrix(int m, int n,const cvector &v) {
            this->m = m;
            this->n = n;
            s = m*n;
            x = new double[s];
            //for(int i=0;i<s;i++) x[i] = 0;
            memset(x,0,s*sizeof(double));
            int l = m;
            if (n<l) l = n;
            if (v.m<l) l = v.m;
            for(int i=0,c=0;i<l;i++,c+=n+1) x[c] = v.x[i]; 
      }
      inline cmatrix(const cvector &v) {
            m = v.m;
            n = 1;
            s = m;
            x = new double[s];
            //for(int i=0;i<s;i++) x[i] = v.x[i];
            memcpy(x,v.x,s*sizeof(double));
      }
            
      inline cmatrix(double *v,int m, int n) {
            this->m = m;
            this->n = n;
            s = m*n;
            //int i;
            x = new double[s];
            //for(i=0;i<s;i++) x[i] = v[i];
            memcpy(x,v,s*sizeof(double));
      }
      // forms a cmatrix of the outer product (ie v1*v2') -- v2 is
      // "transposed" temporarily for this operation
      inline cmatrix(const cmatrix &v1, const cmatrix &v2) {
            if (v1.n!=v2.n) {
                  s = 1;
                  m=1; n=1; x = new double[1];
                  //x[0] = NaN;
                  //x[0] = 0.0/0.0;
                  x[0] = 0;
            } else {
                  n = v2.m; m = v1.m;
                  s = n*m;
                  x = new double[s];
                  int i,j,k,c=0;
                  for(i=0;i<m;i++) for(j=0;j<n;j++,c++) {
                        x[c] = 0;
                        for(k=0;k<v1.n;k++)
                              x[c] += v1.x[i*v1.n+k]*
                                    v2.x[j*v1.n+k];
                  }
            }
      }
      inline cmatrix(const cvector &v1, const cvector &v2) {
            n = v2.m; m = v1.m;
            s = n*m;
            x = new double[s];
            int i,j,c=0;
            for(i=0;i<m;i++) for(j=0;j<n;j++,c++)
                  x[c] = v1.x[i]*v2.x[j];
      }

      inline cmatrix operator-() const {
            cmatrix ret(m,n);
            for(int i=0;i<s;i++) ret.x[i] = -x[i];
            return ret;
      }
      inline cmatrix& operator=(double a) {
            if (a==0) memset(x,0,s*sizeof(double));
            else for(int i=0;i<s;i++) x[i] = a;
            return *this;
      }
      inline cmatrix& operator=(const cmatrix &ma) {
            if (&ma==this) return *this;
            if (ma.n != n || ma.m != m) {
                  s = ma.s; m = ma.m; n = ma.n;
                  delete []x;
                  x = new double[s];
            }
            //for(int i=0;i<s;i++) x[i] = ma.x[i];
            memcpy(x,ma.x,s*sizeof(double));
            return *this;
      }

      inline bool isvalid() const {
            for(int i=0;i<s;i++) if(!finite(x[i])) return false;
            return true;
      }

      inline cmatrix operator*(const cmatrix &ma) const {
            if (n!=ma.m) {
                  cerr << "invalid cmatrix multiply" << endl;
                  assert(0);
            }
            cmatrix ret(m,ma.n);
            int c=0;
            for(int i=0;i<m;i++) for(int j=0;j<ma.n;j++,c++) {
                  ret.x[c] = 0;
                  for(int k=0;k<n;k++)
                        ret.x[c] += x[i*n+k] * ma.x[k*ma.n+j];
            }
            return ret;
      }
      inline cvector operator*(const cvector &v) const {
            if (n!=v.m) {
                  cerr << "invalid cvector-cmatrix multiply" << endl;
                  assert(0);
            }
            cvector ret(m);
            int c = 0;
            for(int i=0;i<m;i++,c+=n) {
                  ret.x[i] = 0;
                  for(int j=0;j<n;j++)
                        ret.x[i] += x[c+j] * v.x[j];
            }
            return ret;
      }

      inline double dot(const cmatrix &ma) const {
            if (n!=ma.n || m!=ma.m) {
                  cerr << "invalid cmatrix dot-product" << endl;
                  assert(0);
            }
            int c = 0;
            double ret = 0.0;
            for(int i=0;i<m;i++,c+=n)
                  for(int j=0;j<n;j++)
                        ret += x[c+j]*ma.x[c+j];
            return ret;
      }
            
      inline void outer(const cmatrix &ma, cmatrix &ret) const {
            if (n!=ma.n || ret.m!=m || ret.n!=ma.m) {
                  cerr << "invalid cmatrix outer multiply" << endl;
                  assert(0);
            }
            int c=0;
            for(int i=0;i<m;i++) for(int j=0;j<ma.m;j++,c++) {
                  ret.x[c] = 0;
                  for(int k=0;k<n;k++)
                        ret.x[c] += x[i*n+k] * ma.x[j*ma.m+k];
            }
      }
      inline void inner(const cmatrix &ma, cmatrix &ret) const {
            if (m!=ma.m || ret.m!=n || ret.n!=ma.n) {
                  cerr << "invalid cmatrix inner multiply" << endl;
                  assert(0);
            }
            int c=0;
            for(int i=0;i<n;i++) for(int j=0;j<ma.n;j++,c++) {
                  ret.x[c] = 0;
                  for(int k=0;k<m;k++)
                        ret.x[c] += x[k*m+i] * ma.x[k*m+j];
            }
      }

      inline double rowmult(int r, const cvector &v, int exclude) const {
            if (n!=v.m) {
                  cerr << "invalid matrix-vector multiply" << endl;
                  assert(0);
            }
            double ret = 0.0;
            int c=n*r;
            for(int j=0;j<n;j++,c++)
                  if (j!=exclude) ret += x[c]*v[j];
            return ret;
      }
      inline double rowmult(int r, const cvector &v) const {
            if (n!=v.m) {
                  cerr << "invalid matrix-vector multiply" << endl;
                  assert(0);
            }
            double ret = 0.0;
            int c=n*r;
            for(int j=0;j<n;j++,c++)
                  ret += x[c]*v[j];
            return ret;
      }

      inline cmatrix &multbyrow(const double *v) {
            int c = 0;
            for(int i=0;i<m;i++)
                  for(int j=0;j<n;j++,c++)
                        x[c] *= v[j];
            return *this;
      }
      inline cmatrix &multbycol(const double *v) {
            int c = 0;
            for(int i=0;i<m;i++)
                  for(int j=0;j<n;j++,c++)
                        x[c] *= v[i];
            return *this;
      }
      inline cmatrix &multbyrow(const cvector &v) {
            if (n!=v.m) {
                  cerr << "invalid multbycol" << endl;
                  assert(0);
            }
            int c = 0;
            for(int i=0;i<m;i++)
                  for(int j=0;j<n;j++,c++)
                        x[c] *= v[j];
            return *this;
      }
      inline cmatrix &multbycol(const cvector &v) {
            if (m!=v.m) {
                  cerr << "invalid multbycol" << endl;
                  assert(0);
            }
            int c = 0;
            for(int i=0;i<m;i++)
                  for(int j=0;j<n;j++,c++)
                        x[c] *= v[i];
            return *this;
      }
      inline cmatrix &dividebyrow(const double *v) {
            int c = 0;
            for(int i=0;i<m;i++)
                  for(int j=0;j<n;j++,c++)
                        x[c] /= v[j];
            return *this;
      }
      inline cmatrix &dividebycol(const double *v) {
            int c = 0;
            for(int i=0;i<m;i++)
                  for(int j=0;j<n;j++,c++)
                        x[c] /= v[i];
            return *this;
      }
      inline cmatrix &dividebyrow(const cvector &v) {
            if (n!=v.m) {
                  cerr << "invalid dividebycol" << endl;
                  assert(0);
            }
            int c = 0;
            for(int i=0;i<m;i++)
                  for(int j=0;j<n;j++,c++)
                        x[c] /= v[j];
            return *this;
      }
      inline cmatrix &dividebycol(const cvector &v) {
            if (m!=v.m) {
                  cerr << "invalid dividebycol" << endl;
                  assert(0);
            }
            int c = 0;
            for(int i=0;i<m;i++)
                  for(int j=0;j<n;j++,c++)
                        x[c] /= v[i];
            return *this;
      }
      
      inline double operator()(int i, int j) const {
            return x[i*n+j];
      }

      inline const double *operator[](int i) const {
            return x+(i*n);
      }
      inline double *operator[](int i) {
            return x+(i*n);
      }
      inline cmatrix t() const {
            return cmatrix(*this,true);
      }

      inline cmatrix &operator+=(const cmatrix &ma) {
            if (m!=ma.m||n!=ma.n) {
                  cerr << "invalid cmatrix addition" << endl;
                  assert(0);
            }
            for(int i=0;i<s;i++) x[i] += ma.x[i];
            return *this;
      }

      inline cmatrix &operator-=(const cmatrix &ma) {
            if (m!=ma.m||n!=ma.n) {
                  cerr << "invalid cmatrix addition" << endl;
                  assert(0);
            }
            for(int i=0;i<s;i++) x[i] -= ma.x[i];
            return *this;
      }

      inline cmatrix &operator*=(const cmatrix &ma) {
            if (n!=ma.m || n != ma.n) {
                  cerr << "invalid cmatrix multiply" << endl;
                  assert(0);
            }
            int i,j,k,c=0;
            double newrow[n];
            for(i=0;i<m;i++) {
              for(j=0;j<n;j++) {
                newrow[j] = 0;
                for(k=0;k<n;k++)
                  newrow[j] += x[c+k] * ma.x[k*n+j];
              }
              for(j=0; j < n; j++,c++)
                x[c] = newrow[j];
            }
            return *this;
      }

      inline cmatrix &operator+=(const double &a) {
            for(int i=0;i<s;i++) x[i] += a;
            return *this;
      }

      inline cmatrix &operator-=(const double &a) {
            for(int i=0;i<s;i++) x[i] -= a;
            return *this;
      }

      inline cmatrix &operator*=(const double &a) {
            for(int i=0;i<s;i++) x[i] *= a;
            return *this;
      }

      inline cmatrix &operator/=(const double &a) {
            for(int i=0;i<s;i++) x[i] /= a;
            return *this;
      }

      inline double max() const {
            double t,ma = x[0];
            for(int i=1;i<s;i++) 
                  if ((t=x[i])>ma) ma=t;
            return ma;
      }

      inline double min() const {
            double t,mi = x[0];
            for(int i=1;i<s;i++) 
                  if ((t=x[i])<mi) mi=t;
            return mi;
      }

      inline double absmin() const {
            double t,mi = x[0]>0?x[0]:-x[0];
            for(int i=1;i<s;i++)
                  if ((t=(x[i]>0?x[i]:-x[i]))<mi) mi=t;
            return mi;
      }

      inline double absmax() const {
            double t,ma = x[0]>0?x[0]:-x[0];
            for(int i=1;i<s;i++)
                  if ((t=(x[i]>0?x[i]:-x[i]))>ma) ma=t;
            return ma;
      }

      inline bool operator==(const cmatrix &ma) const {
            if (ma.n!=n||ma.m!=m) return false;
#if !defined(HAVE_BCMP)
            for(int i=0;i<s;i++) if (ma.x[i]!=x[i]) return false;
            return true;
#else
            return bcmp(ma.x,x,s*sizeof(double))==0.0;
#endif // HAVE_BCMP
      }
      inline bool operator==(const double &a) const {
            for(int i=0;i<s;i++) if (x[i]!=a) return false;
            return true;
      }
      inline bool operator!=(const cmatrix &ma) const {
            if (ma.n!=n||ma.m!=m) return true;
#if !defined(HAVE_BCMP)
            for(int i=0;i<s;i++) if (ma.x[i]!=x[i]) return true;
            return false;
#else
            return bcmp(ma.x,x,s*sizeof(double))!=0.0;
#endif // HAVE_BCMP
      }
      inline bool operator!=(const double &a) const {
            for(int i=0;i<s;i++) if (x[i]!=a) return true;
            return false;
      }

      inline double norm2() const { // returns the square of the frobenius norm
            double ret=x[0]*x[0];
            for(int i=1;i<s;i++) ret += x[i]*x[i];
            return ret;
      }
      inline double norm() const { // returns the frobenius norm
            return sqrt(norm2());
      }

      friend ostream& operator<<(ostream& s, const cmatrix& ma);
      friend istream& operator>>(istream& s, cmatrix& ma);

      // LU decomposition -- ix is the row permutations
      int LUdecomp(cmatrix &LU, int *ix) const;
      // LU back substitution --
      //    ix from above fn call (this should be an LU combination)
      void LUbacksub(int *ix, double *col) const;

      // solves equation Ax=b (A is this, x is the returned value)
bool solve(cvector &b, cvector &dest);
      double *solve(const double *b, bool &worked) const;
      inline double *solve(const double *b) const { bool w; return solve(b,w); }
      
      inline void negate() { for(int i = 0; i < s; i++) x[i] = -x[i]; }
      cmatrix inv(bool &worked) const;
      inline cmatrix inv() const { bool w; return inv(w); }
      double adjoint();
      inline double trace();
      double testAdjoint();
      inline void multiply(const cvector &source, cvector &dest) {
        assert(n == source.m && m == dest.m);
        int i,j,c=0;
        for(i = 0; i < m; i++) {
          dest[i] = 0;
          for(j = 0; j < n; j++,c++)
            dest[i] += x[c] * source[j];
        }
      }

      inline double *values() { return x; }



      // w needs to points to an array of n fp numbers
      void svd(cmatrix &u, cmatrix &v, double *w);

       // makes the matrix tri-diagonal (if symmetric)
      // The matrix is replaced by Q, and the vectors d and e hold
      // the diagonal and off-diagonal elements respectively
      void tridiag(cvector &d, cvector &e);

      // converts the Q matrix from tridiag into the eigenvectors
      // d and e are the diagonal and off-diagonal elements from
      // tridiag and the eigenvalues are left in d
      void tridiagev(cvector &d, cvector &e);

      // returns the inverse of a real, symmetric matrix
      // via eigenvalue decomposition.  If any of the e.v. are
      // very close to zero, it uses the psuedo inverse
      // (ie 1/ev = 0)
      cmatrix syminv(double tol=1e-10) const;

      // uses the above two to compute eigen decomp of symmetric
      // matrix (destroys matrix and replaces with eigenvectors)
      cvector eigensym() { cvector d,e; tridiag(d,e); tridiagev(d,e); return d; }


      inline ostream &niceprint(ostream& s) {
            s << m << ' ' << n << endl;
            for(int i=0;i<m;i++) {
                  for(int j=0;j<n;j++)
                        s << x[i*n+j] << ' ';
                  s << endl;
            }
            return s;
      }

      inline int getm() const { return m; }
      inline int getn() const { return n; }

      inline void compact() { }

private:
      static double pythag(double a, double b);

      int m,n,s;
      double *x;
};

inline cmatrix operator+(const cmatrix &a, const cmatrix &b) {
      return cmatrix(a)+=b;
}
inline cmatrix operator-(const cmatrix &a, const cmatrix &b) {
      return cmatrix(a)-=b;
}
inline cmatrix operator+(const cmatrix &a, const double &b) {
      return cmatrix(a)+=b;
}
inline cmatrix operator-(const cmatrix &a, const double &b) {
      return cmatrix(a)-=b;
}
inline cmatrix operator+(const double &a, const cmatrix &b) {
      return cmatrix(b)+=a;
}
inline cmatrix operator-(const double &a, const cmatrix &b) {
      return cmatrix(b.getn(),b.getm(),a)-=b;
}
inline cmatrix operator*(const cmatrix &a, const double &b) {
      return cmatrix(a)*=b;
}
inline cmatrix operator*(const double &b, const cmatrix &a) {
      return cmatrix(a)*=b;
}
inline cmatrix operator/(const cmatrix &a, const double &b) {
      return cmatrix(a)/=b;
}

inline ostream& operator<<(ostream& s, const cmatrix& ma) {
//    s << ma.m << ' ' << ma.n << ' ';
      for(int i=0;i<ma.s;i++) { 
if (i%ma.n==0) s << endl; s << ma.x[i]; if (i!=ma.s) s << ' '; }
      return s;
}

inline istream& operator>>(istream& s, cmatrix& ma) {
      int tn,tm;
      s >> tm >> tn;
      if (tm!=ma.m || tn!=ma.n) {
            delete []ma.x;
            ma.s = tm*tn;
            ma.x = new double[ma.s];
            ma.m = tm; ma.n = tn;
      }
      for(int i=0;i<ma.s;i++) { s >> ma.x[i]; }
      return s;
}
#endif

Generated by  Doxygen 1.6.0   Back to index