#include "MAT.h"
#define EPS 0.000000001
#define ITERATIONS 20000
static int freadv(FILE *fp, int *value) { return fscanf(fp,"%d", value); }
static int freadv(FILE *fp, double *value) { return fscanf(fp,"%lf", value); }
static int fprintv(int v){ return printf("%d ", v); }
static int fprintv(double v){ return printf("%lf ", v); }
int dluav(double a[],int m,int n,double u[],double v[],double eps); //svd
int eejcb(double a[],int n,double v[],double eps,int jt); //eig
template <class T>
void VEC<T>::Print(char *msg)
{
int i, t;
printf(msg);
printf("T= %d\n", (t=(int)(*this)));
for (i=0; i < t; i++)
fprintv((*this)[i]);
printf("\n");
}
template <class T>
void VEC<T>::Read(char *fl)
{
int i, t;
FILE *fp;
if((fp = fopen(fl, "r"))==0) throw new exception("File open error!\n");
fscanf(fp, "T= %d\n", &t);
VEC<T> s(t);
for (i=0; i < t; i++)
freadv(fp, &s[i]);
fclose(fp);
*this=s;
}
template <class T>
VEC<T>::VEC():r(0), vec(0){ }
template <class T>
VEC<T>::VEC(int r): vec(new T[r]), r(vec?r:0){ }
template <class T>
VEC<T>::VEC(const VEC<T>& m)
{
if((vec=new T[m.r])==0) throw new exception("Memory allocate error!\n");
for(int x=(r=m.r)-1; x>=0; x--)
vec[x]=m.vec[x];
}
template <class T>
VEC<T>& VEC<T>::operator=(const VEC<T>& m)
{
VEC<T>::~VEC();
if((vec=new T[m.r])==0) throw new exception("Memory allocate error!\n");
for(int x=(r=m.r)-1; x>=0; x--)
vec[x]=m.vec[x];
return *this;
}
template <class T>
VEC<T>::operator int()const { return r; }
template <class T>
T & VEC<T>::operator[](int ind)const
{
if(ind>=r) throw new exception("Index error!\n");
return vec[ind];
}
template <class T>
VEC<T>::~VEC(){ if(vec){ delete []vec; r=0; vec=0; } }
template <class T>
void MAT<T>::Print(char *msg) const
{
int i, j;
printf(msg);
printf("r= %d, c= %d\n", r, c);
for (i=0; i < r; i++)
{
for(j=0; j< c; j++)
fprintv(mat[i*c+j]);
printf("\n");
}
printf("\n");
}
template <class T>
MAT<T>::MAT(): r(0), c(0), mat(0) { }
template <class T>
MAT<T>::MAT(int r, int c): mat(new T[r*c]), r(mat?r:0), c(mat?c:0) { }
template <class T>
MAT<T>::MAT(const MAT<T>& m)
{
if((mat=new T[m.r*m.c])==0) throw new exception("Memory allocate error!\n");
r=m.r; c=m.c;
for(int x=r*c-1; x>=0; x--)
mat[x]=m.mat[x];
}
template <class T>
MAT<T>& MAT<T>::operator=(const MAT<T> &m)
{
MAT<T>::~MAT( );
if((mat=new T[m.r*m.c])==0) throw new exception("Memory not enough1\n");
r=m.r; c=m.c;
for(int x=r*c-1; x>=0; x--)
mat[x]=m.mat[x];
return *this;
}
template <class T>
T *const MAT<T>::operator[](int ind) const
{
if(ind>=r) throw new exception("Matrix index bound error!\n");
return mat+ind*c;
}
template <class T>
MAT<T> MAT<T>::operator+(const MAT<T> &m)const //for +
{
if(r!=m.r || c!=m.c) throw new exception("Matrix index bound error!\n");
MAT<T> s(*this);
for(int h=0; h<r; h++)
for(int k=0; k<c; k++){
s[h][k]+=m[h][k];
}
return s;
}
template <class T>
MAT<T> MAT<T>::operator+(T m)const //for +
{
MAT<T> s(*this);
for(int h=0; h<r; h++)
for(int k=0; k<c; k++){
s[h][k]+=m;
}
return s;
}
template <class T>
MAT<T> MAT<T>::operator-(const MAT<T> &m)const //for -
{
if(r!=m.r || c!=m.c) throw new exception("Matrix index bound error!\n");
MAT<T> s(*this);
for(int h=0; h<r; h++)
for(int k=0; k<c; k++){
s[h][k]-=m[h][k];
}
return s;
}
template <class T>
MAT<T> MAT<T>::operator-(T m)const //for -
{
MAT<T> s(*this);
for(int h=0; h<r; h++)
for(int k=0; k<c; k++){
s[h][k]-=m;
}
return s;
}
template <class T>
MAT<T> operator-(T m, const MAT<T> n)//for -
{
MAT<T> s(n.r, n.c);
for(int h=0; h<r; h++)
for(int k=0; k<c; k++){
s[h][k]=m-n[h][k];
}
return s;
}
template <class T>
MAT<T> MAT<T>::operator*(const MAT<T> &m)const //for *
{
if(c!=m.r) throw new exception("Matrix index bound error!\n");
MAT<T> s(r, m.c);
for(int h=0; h<r; h++)
for(int j=0; j<m.c; j++){
s.mat[h*s.c+j]=0;
for(int k=0; k<c; k++)
s.mat[h*s.c+j]+=mat[h*c+k]*m.mat[k*m.c+j];
}
return s;
}
template <class T>
MAT<T> MAT<T>::dot_multiply(const MAT<T> &m)const //for .*
{
if((c!=m.c)||(r!=m.r)) throw new exception("Matrix index bound error!\n");
MAT<T> s(r, c);
for(int h=r*c-1; h>=0; h--)
s.mat[h]=mat[h]*m.mat[h];
return s;
}
template <class T>
MAT<T> MAT<T>::dot_multiply(T m)const //for .*
{
MAT<T> s(*this);
for(int h=r*c-1; h>=0; h--) s.mat[h]*=m;
return s;
}
template <class T>
MAT<T> MAT<T>::operator^(int p)const //for power
{
if(r!=c) throw new exception("Must be a square Matrix!\n");
MAT<T> s(r,c);
for(int h=0; h<r; h++) s[h][h]=1;
for(int h=0; h<p; h++) s=s*(*this);
return s;
}
template <class T>
MAT<T> MAT<T>::dot_power(double p)const //for dot_power
{
MAT<T> s(r,c);
for(int h=0; h<r; h++)
for(int k=0; k<c; k++)
s[h][k]=pow((*this)[h][k],p);
return s;
}
template <class T>
MAT<T> MAT<T>::dot_divide(T m)const //for ./
{
if(m==0) throw new exception("Divide by Zero!\n");
MAT<T> s(*this);
for(int h=r*c-1; h>=0; h--)
s.mat[h]/=m;
return s;
}
template <class T>
MAT<T> MAT<T>::dot_divide(const MAT<T> &m)const //for ./
{
if((c!=m.c)||(r!=m.r)) throw new exception("Matrix index bound error!\n");
MAT<T> s(r, c);
for(int h=r*c-1; h>=0; h--)
s.mat[h]=mat[h]/m.mat[h];
return s;
}
template <class T>
MAT<T> MAT<T>::dot_divide(T x, const MAT<T> &m) //for ./
{
for( int h=m.r*m.c-1; h>=0; h--)
if(m.mat[h]==0) throw new exception("Matrix index bound error!\n");
MAT<T> s(m.r,m.c);
for(int h=m.r*m.c-1; h>=0; h--)
s.mat[h]=x/m.mat[h];
return s;
}
template <class T>
MAT<T> MAT<T>::transpose()const //for Transpose
{
MAT<T> s(c, r);
for(int i=0; i<r; i++)
for(int j=0; j<c; j++)
s[j][i]=(*this)[i][j];
return s;
}
template <class T>
void MAT<T>::svd(MAT<T> &U,MAT<T> &S,MAT<T> &V)const //for Singular value decompose
{
MAT<T> a(*this);
MAT<T> u(r, r), v(c,c);
dluav(a.mat, r, c, u.mat, v.mat, EPS);
U=u;
S=a;
V=*v;
}
template <class T>
int MAT<T>::eig(MAT<T> &V,MAT<T> &D)const //for Singular value decompose
{
if(r!=c) throw new exception("Must be a square matrix!\n");
MAT<T> D(*this);
MAT<T> V(r, r);
int k, h;
k=eejcb(D.mat, r, V.mat, EPS, ITERATIONS);
if(k>=0) {
for(k=h=0; h<r; h++) if(fabs(D[h][h])>=EPS*0.1) k++;
}
return k;
}
template <class T>
int MAT<T>::rank()const
{
int i, j, l, is, js, nn, Rank;
T d, q, t;
MAT<T> M(*this);
nn=r>c?c:r;
Rank=0;
for(l=0;l<nn;l++)
{
q=0;
for(i=l;i<r;i++)
for(j=l;j<c;j++)
{
d=fabs(M[i][j]);
if(d>q) {q=d;is=i;js=j;}
}
if(q==0) break;
Rank=Rank+1;
if(js!=l) for(i=l;i<r;i++) {t=M[i][l];M[i][l]=M[i][js];M[i][js]=t;}
if(is!=l) for(j=l;j<c;j++) {t=M[l][j];M[l][j]=M[is][j];M[is][j]=t;}
for(i=l+1;i<r;i++)
for(j=l+1;j<c;j++)
{
d=M[i][l]*M[l][j]/M[l][l];
M[i][j]=M[i][j]-d;
}
}
return Rank;
}
template <class T>
T MAT<T>::det( ) const
{
if(r!=c) throw new exception("Must be square matrix!\n");