C++类写的矩阵计算

#pragma once
#include "DibImage.h"
#include "Exception.h"
//矩阵类-声明

class Matrix
{
private:
 //矩阵宽度
 long width;
 //矩阵高度
 long height;
protected:
 //矩阵的首地址
 double * p;
public:
 static Matrix One(long n);
 void info();//显示信息
 long getHeight();
 long getWidth();
 bool isVector();//如果是false，那就是一个数
 double Deter();//求行列式determinant
 bool isPtv();//是否正定
 Matrix Eigen();//特征值
 Matrix Adjoint();//伴随矩阵
 Matrix  T();//转置
 Matrix  Inv();//逆矩阵
 void Unity();//归一化
 void Store(const char * filename);//把矩阵存储到文件里
 Matrix  BiCubic(long n);//构造如下类似矩阵
        /*{{2,0.5,0,0,0,0},
       {0.5,2,0.5,0,0,0},
       {0,0.5,2,0.5,0,0},
       {0,0,0.5,2,0.5,0},
       {0,0,0,0.5,2,0.5,},
       {0,0,0,0,0.5,2,}
       };*/
 Matrix subMatrix(long offset) ;//throw (Exception &);
 Matrix subMatrix(long x,long y,long WidthOffset,long HeightOffset);//取整数位置子区
 //Matrix subMatrix(double x,double y,long WidthOffset,long HeightOffset);//取小数位置子区
 void Test();//测试函数
    /***********重载部分-Overloaded Part*****************/
 Matrix operator+(Matrix &);
 Matrix operator-(Matrix &);
 Matrix operator*(Matrix &);
 friend Matrix operator*(double alpha,Matrix &);//实数与矩阵相乘
 Matrix operator*(double alpha);//矩阵与实数相乘
 Matrix operator/(Matrix &);//实际是实数相除或矩阵和实数相除
 Matrix operator/(double sub);
 Matrix operator+=(Matrix &);
 Matrix operator-=(Matrix &);
 Matrix operator*=(Matrix &);//矩阵与实数相乘
 Matrix operator*=(double alpha);//矩阵与实数相乘
 Matrix &operator=(Matrix &);//赋值
 double * operator[](long heightPos);//用于实现用[][]操作矩阵元素
 friend Matrix  sqrt(Matrix m);//开方
 friend double abs(Matrix &);//取绝对值(泛数)
 friend double sum(Matrix &);//求和
 friend Matrix multiply(Matrix &m1,Matrix & m2);//按元素相乘
 friend double operator+(double dbl,Matrix &);
 friend double operator+(Matrix &,double dbl);
 friend double operator-(double dbl,Matrix &);
 friend double operator-(Matrix &,double dbl);
 friend ostream & operator<<(ostream &,Matrix &);//输出
 Matrix(void);//默认构造函数
 Matrix(long n);//单位矩阵
 Matrix(DibImage & dibImage);//类型转换(把图像文件转换成矩阵，用于图像计算处理)
 Matrix(double * arrAddress,long arrWidth);//构造一维矩阵(一行，arrWidth列)
 Matrix(double * arrAddress,long arrHeight,long arrWidth);//构造二维矩阵
 Matrix(long Height,long Width);//构造空矩阵
 Matrix(const Matrix &);//复制构造函数
 virtual ~Matrix(void);//默认析构函数
 /***********重载部分-Overloaded Part*****************/
};
#include "StdAfx.h"
#include "./matrix.h"
//矩阵类-实现
#define SIGN(a,b) ((b) > 0 ? fabs(a) : -fabs(a))
Matrix::Matrix(void):width(1),height(1)
{
 p=new double[1];
 //width=1;
 //height=1;
 //cout<<"Default Constructor Invoked."<<endl;
}
Matrix::Matrix(long n)
{
 height=width=n;
 p=new double[n*n];
 for(long i=0;i<n;i++){
  for(long j=0;j<n;j++){
   if(i==j) *(p+n*i+j)=1;
   else *(p+n*i+j)=0;
  }
 }
}

Matrix::Matrix(DibImage & dibImage)
{
 
 width=dibImage.getWidth();height=dibImage.getHeight();
 p=new double[width*height];
 if(p==NULL) assert(0);//内存分配失败
 long lLineBytes = WIDTHBYTES(width * 8);
 LPSTR Src = ::FindDIBBits((LPSTR)::GlobalLock((HGLOBAL) dibImage.GetHDIB()));
 for(long i=0;i<height;i++){
  for(long j=0;j<width;j++){
   *(p+width*i+j)=double(*((unsigned char *)Src + lLineBytes * (height-1-i) + j));
  // *(p+width*i+j)=double(*((unsigned char *)Src + lLineBytes * (i) + j));
  }
 }
 ::GlobalUnlock((HGLOBAL) dibImage.GetHDIB());
}
Matrix::Matrix(double * arrAddress,long arrWidth)
{
 long arrHeight=1;
 p=new double[arrWidth*arrHeight];
 for(long i=0;i<arrHeight;i++){
  for(long j=0;j<arrWidth;j++){
   *(p+arrWidth*i+j)=*(arrAddress+arrWidth*i+j);
  }
 }
 width=arrWidth;
 height=arrHeight;
}
Matrix::Matrix(double * arrAddress,long arrHeight,long arrWidth)
{
 p=new double[arrWidth*arrHeight];
 for(long i=0;i<arrHeight;i++){
  for(long j=0;j<arrWidth;j++){
   *(p+arrWidth*i+j)=*(arrAddress+arrWidth*i+j);
  }
 }
 width=arrWidth;
 height=arrHeight;
}
Matrix::Matrix(long Height,long Width)
{
 p=new double[Height*Width];
 width=Width;
 height=Height;
}
Matrix::Matrix(const Matrix & m)//copy constructor
{
 height=m.height;
 width=m.width;
 p=new double[height*width];
 for(long i=0;i<height;i++){
   for(long j=0;j<width;j++){
    *(p+width*i+j)=*(m.p+width*i+j);
   }
 }
}
long Matrix::getWidth(){
 return width;
}
long Matrix::getHeight(){
 return height;
}
Matrix Matrix::One(long n)
{
 Matrix m(1,n);
 for(long i=0;i<n;i++)
  m[0][i]=i;
 return m;
}
bool Matrix::isVector()
{
 return !(width==1 && height==1);
}
Matrix Matrix::subMatrix(long offset) //throw (Exception &)
{
 if(!(height==width && offset<=width && offset>=0))
  throw Exception("Error at Matrix Matrix::subMatrix(long offset)./n NOT height==width && offset<=width && offset>=0");
 double * t=new double[offset*offset];
 for(long i=0;i<offset;i++){
  for(long j=0;j<offset;j++){
   *(t+offset*i+j)=*(p+width*i+j);
  }
 }
 Matrix m(t,offset,offset);
 delete [] t;
 return m;
}
Matrix Matrix::subMatrix(long x,long y,long WidthOffset,long HeightOffset)
{
    //cout<<"x="<<x<<" y="<<y<<" w="<<WidthOffset<<"  h="<<HeightOffset<<"  in subMatrix 1/n";
 bool log=true;
 if(x<0 || y<0){
  cout<<"x<0 or y<0 in Matrix::subMatrix()!/n";
  log=false;
 }else if(WidthOffset<0 || HeightOffset<0){  
  cout<<"WidthOffset<0 || HeightOffset<0 in Matrix::subMatrix()!/n";
  log=false;
 }else if(WidthOffset*2+1>width || HeightOffset*2+1>height){
  cout<<">Subarea is bigger than Matrix in Matrix::subMatrix()!/n";
  log=false;
 }else if((x-WidthOffset)<0 || (x+WidthOffset)>width){
  cout<<"Width error in Matrix::subMatrix()!/n";
  log=false;
 }else if((y-HeightOffset)<0 || (y+HeightOffset)>height){
  cout<<"Height error in Matrix::subMatrix()!/n";
  log=false;
 }
 assert(log);
 Matrix m(HeightOffset*2+1,WidthOffset*2+1);
 for(long i=0;i<HeightOffset*2+1;i++){
  for(long j=0;j<WidthOffset*2+1;j++){
   m[i][j]=(*this)[i+y-HeightOffset][j+x-WidthOffset];
  }
 }
 return m; 
}

double Matrix::Deter()//矩阵的行列式
{
 assert(width==height);
 double result=1;
 double k;
 Matrix m=*this;
 for(long i=0;i<height-1;i++){//Gauss消去法，变换成上三角矩阵
  for(long j=i+1;j<height;j++){
   k=m[j][i]/m[i][i];
   m[j][i]=0;
   for(long n=i+1;n<width;n++){
    m[j][n]=m[j][n]-k*m[i][n];
   }
  }
 }
 for(long i=0;i<height;i++){
  for(long j=0;j<width;j++){
   if(i==j) result*=m[i][j];
  }
 }
 return result;
}

bool Matrix::isPtv()
{
 assert(width==height);//是方阵才可以计算
 bool result=true;
 Matrix m;
 for(long i=1;i<=height;i++){
   m=this->subMatrix(i);
   if(m.Deter()<=0){
    result=false;
    break;
   }
 }
 return result;
}

Matrix  Matrix::T()
{
 double * t=new double[width*height];
 for(long i=0;i<height;i++){
  for(long j=0;j<width;j++){
   *(t+height*j+i)=*(p+width*i+j);
  }
 }
 Matrix m(t,width,height);
 delete [] t;
 return m;
}
Matrix Matrix::Eigen()
{
 Matrix matrix=*this;
 assert(matrix.getHeight()==matrix.getWidth());

 double **a,*wr,*wi;
 long nn,m,l,k,j,its,i,mmin;
 double z,y,x,w,v,u,t,s,r,q,p,anorm;
 
 
    long  n=matrix.getWidth();
 a=new double*[n+1];
 wr=new double[n+1];
 wi=new double[n+1];
 for(i=0;i<n+1;i++)
  a[i]=new double[n+1];

 for(i=0;i<n;i++){
  for(j=0;j<n;j++){
   a[i+1][j+1]=matrix[i][j];
  }
 }

 anorm=fabs(a[1][1]);
 for (i=2;i<=n;i++)
  for (j=(i-1);j<=n;j++)
   anorm += fabs(a[i][j]);
 nn=n;
 t=0.0;
 while (nn >= 1) {
  its=0;
  do {
   for (l=nn;l>=2;l--) {
    s=fabs(a[l-1][l-1])+fabs(a[l][l]);
    if (s == 0.0) s=anorm;
    if (fabs(a[l][l-1]) + s == s) break;
   }
   x=a[nn][nn];
   if (l == nn) {
    wr[nn]=x+t;
    wi[nn--]=0.0;
   } else {
    y=a[nn-1][nn-1];
    w=a[nn][nn-1]*a[nn-1][nn];
    if (l == (nn-1)) {
     p=0.5*(y-x);
     q=p*p+w;
     z=sqrt(fabs(q));
     x += t;
     if (q >= 0.0) {
      z=p+SIGN(z,p);
      wr[nn-1]=wr[nn]=x+z;
      if (z) wr[nn]=x-w/z;
      wi[nn-1]=wi[nn]=0.0;
     } else {
      wr[nn-1]=wr[nn]=x+p;
      wi[nn-1]= -(wi[nn]=z);
     }
     nn -= 2;
    } else {
     if (its == 30) printf("Too many iterations in HQR");
     if (its == 10 || its == 20) {
      t += x;
      for (i=1;i<=nn;i++) a[i][i] -= x;
      s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);
      y=x=0.75*s;
      w = -0.4375*s*s;
     }
     ++its;
     for (m=(nn-2);m>=l;m--) {
      z=a[m][m];
      r=x-z;
      s=y-z;
      p=(r*s-w)/a[m+1][m]+a[m][m+1];
      q=a[m+1][m+1]-z-r-s;
      r=a[m+2][m+1];
      s=fabs(p)+fabs(q)+fabs(r);
      p /= s;
      q /= s;
      r /= s;
      if (m == l) break;
      u=fabs(a[m][m-1])*(fabs(q)+fabs(r));
      v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
      if (u+v == v) break;
     }
     for (i=m+2;i<=nn;i++) {
      a[i][i-2]=0.0;
      if  (i != (m+2)) a[i][i-3]=0.0;
     }
     for (k=m;k<=nn-1;k++) {
      if (k != m) {
       p=a[k][k-1];
       q=a[k+1][k-1];
       r=0.0;
       if (k != (nn-1)) r=a[k+2][k-1];
       if (x=fabs(p)+fabs(q)+fabs(r)) {
        p /= x;
        q /= x;
        r /= x;
       }
      }
      if (s=SIGN(sqrt(p*p+q*q+r*r),p)) {
       if (k == m) {
        if (l != m)
        a[k][k-1] = -a[k][k-1];
       } else
        a[k][k-1] = -s*x;
       p += s;
       x=p/s;
       y=q/s;
       z=r/s;
       q /= p;
       r /= p;
       for (j=k;j<=nn;j++) {
        p=a[k][j]+q*a[k+1][j];
        if (k != (nn-1)) {
         p += r*a[k+2][j];
         a[k+2][j] -= p*z;
        }
        a[k+1][j] -= p*y;
        a[k][j] -= p*x;
       }
       mmin = nn<k+3 ? nn : k+3;
       for (i=l;i<=mmin;i++) {
        p=x*a[i][k]+y*a[i][k+1];
        if (k != (nn-1)) {
         p += z*a[i][k+2];
         a[i][k+2] -= p*r;
        }
        a[i][k+1] -= p*q;
        a[i][k] -= p;
       }
      }
     }
    }
   }
  } while (l < nn-1);
 }
 Matrix eigenvalue(1,n);
 for(i=0;i<n;i++){
  eigenvalue[0][i]=wr[i+1];
 }
 return eigenvalue;
}
Matrix Matrix::Adjoint()
{
 return this->Inv()*this->Deter();
}
void Matrix::Unity()
{
 double total;
 total=sqrt(sum(multiply(*this,*this)));
 for(long i=0;i<this->height;i++){
  for(long j=0;j<this->width;j++){
   (*this)[i][j]/=total;
  }
 }
    return;
}
void Matrix::Store(const char * filename){
 ofstream fout(filename);
 int w=5;
 for(long i=0;i<height;i++){
  for(long j=0;j<width;j++){
   fout.width(w);
   fout<<i;
   fout.width(w);
   fout<<j;
   fout.width(15);
   fout<<(*this)[i][j]<<endl;/*
   cout.width(width);
   cout<<i;
   cout.width(width);
   cout<<j;
   cout.width(width);
   cout<<(*this)[i][j]<<endl;*/
  }
 }
 fout.close();
}
Matrix Matrix::Inv()
{
 assert(width==height && width>0);
 Matrix m=*this;
 double * pDbSrc=m.p;
 long nLen=m.width;
 int *is,*js,i,j,k;
 double d,p;
 is = new int[nLen];
 js = new int[nLen];
 for(k=0;k<nLen;k++)
 {
  d=0.0;
  for(i=k;i<nLen;i++)
  for(j=k;j<nLen;j++)
  { 
   p=fabs(pDbSrc[i*nLen + j]);
   if(p>d)
   {
    d     = p;
    is[k] = i;
    js[k] = j;
   }
  }
  if(d+1.0==1.0)
  {
   delete is;
   delete js;
   assert(0);
   //return FALSE;
  }
  if(is[k] != k)
  for(j=0;j<nLen;j++)
  {
   p = pDbSrc[k*nLen + j];
   pDbSrc[k*nLen + j] = pDbSrc[(is[k]*nLen) + j];
   pDbSrc[(is[k])*nLen + j] = p;
  }
  if(js[k] != k)
  for(i=0; i<nLen; i++)
  {
   p = pDbSrc[i*nLen + k];
   pDbSrc[i*nLen + k] = pDbSrc[i*nLen + (js[k])];
   pDbSrc[i*nLen + (js[k])] = p;
  }
  pDbSrc[k*nLen + k]=1.0/pDbSrc[k*nLen + k];
  for(j=0; j<nLen; j++)
   if(j != k)
   {
    pDbSrc[k*nLen + j]*=pDbSrc[k*nLen + k];
   }
  for(i=0; i<nLen; i++)
   if(i != k)
    for(j=0; j<nLen; j++)
     if(j!=k)
     { 
      pDbSrc[i*nLen + j] -= pDbSrc[i*nLen + k]*pDbSrc[k*nLen + j];
     }
  for(i=0; i<nLen; i++)
   if(i != k)
   {
    pDbSrc[i*nLen + k] *= -pDbSrc[k*nLen + k];
   }
  }
  for(k=nLen-1; k>=0; k--)
  {
   if(js[k] != k)
    for(j=0; j<nLen; j++)
    {
     p = pDbSrc[k*nLen + j];
     pDbSrc[k*nLen + j] = pDbSrc[(js[k])*nLen + j];
     pDbSrc[(js[k])*nLen + j] = p;
    }
  if(is[k] != k)
   for(i=0; i<nLen; i++)
   {
    p = pDbSrc[i*nLen + k];
    pDbSrc[i*nLen + k] = pDbSrc[i*nLen +(is[k])];
    pDbSrc[i*nLen + (is[k])] = p;
   }
  }
  delete is;
  delete js;
  return m;
}
Matrix  Matrix::BiCubic(long n)
{
 assert(n>1);
 delete [] p;
 width=n;
 height=n;
 p=new double[n*n];
 //首行
 *(p)=0.5; *(p+1)=2;
 for(long k=2;k<n;k++) *(p+k)=0;
 //中间行
 for(long i=1;i<n-1;i++){//行
  for(long j=0;j<n;j++){//列
   if(j==(i-1)){
    *(p+n*i+j)=0.5;
    j++;
    *(p+n*i+j)=2;
    j++;
    *(p+n*i+j)=0.5;
    j++;
   }
   *(p+n*i+j)=0;
  }
 }
 //末行
 for(long l=0;l<n-2;l++) *(p+n*(n-1)+l)=0;
 *(p+n*n-1)=2;
 *(p+n*n-2)=0.5;

 return *this;
}
Matrix Matrix::operator +(Matrix &m1)
{
 assert(m1.height==height && m1.width==width);
 long tmpHeight=m1.height;
 long tmpWidth=m1.width;
 double * t=new double[tmpWidth*tmpHeight];
 for(long i=0;i<tmpHeight;i++){
  for(long j=0;j<tmpWidth;j++){
   *(t+tmpWidth*i+j)=*(m1.p+tmpWidth*i+j)+*(p+tmpWidth*i+j);
  }
 }
 Matrix m(t,tmpHeight,tmpWidth);
 delete [] t;
 return m;
}
Matrix Matrix::operator -(Matrix &m1)
{
 assert(m1.height==height && m1.width==width);
 long tmpHeight=m1.height;
 long tmpWidth=m1.width;
 double * t=new double[tmpWidth*tmpHeight];
 for(long i=0;i<tmpHeight;i++){
  for(long j=0;j<tmpWidth;j++){
   *(t+tmpWidth*i+j)=*(p+tmpWidth*i+j)-*(m1.p+tmpWidth*i+j);
  }
 }
 Matrix m(t,tmpHeight,tmpWidth);
 delete [] t;
 return m;
}
Matrix Matrix::operator *(Matrix &m1)
{
 if(!this->isVector() && m1.isVector()){//左为数,右为矩阵
  Matrix m;
   m=p[0]*m1;
  return m;
 }else if(this->isVector() && !m1.isVector()){//左为矩阵,右为数
  Matrix m;
   m=*this*m1[0][0];
  return m;
 }else if(!this->isVector() && m1.isVector()){//左右都为数
  double * t=new double[1];
  t[0]=p[0]*m1[0][0];
  Matrix m(t,1,1);
  delete [] t;
  return m;
 }else if(this->isVector() && m1.isVector() && width==m1.height){//左为矩阵,右为矩阵
  double sum;
  double * t=new double[height*m1.width];
  for(long i=0;i<height;i++){
   for(long j=0;j<m1.width;j++){
    sum=0;
    for(long k=0;k<width;k++){
     sum+=(*(p+width*i+k))*(m1[k][j]);
    }
    *(t+m1.width*i+j)=sum;
   }
  }
  Matrix m(t,height,m1.width);
  delete [] t;
  return m;
 }else{
  assert(0);//未知运算
  return *this;
 }
}
Matrix operator*(double alpha,Matrix & m1)
{
 Matrix m=m1;
 for(long i=0;i<m.height;i++){
  for(long j=0;j<m.width;j++){
   m[i][j]=alpha*m1[i][j];
  }
 }
 return m;
}
Matrix Matrix::operator*(double alpha)
{
 return alpha*(*this);
}
Matrix Matrix::operator+=(Matrix & m)
{
 return *this+m;
}
Matrix Matrix::operator-=(Matrix & m)
{
 return *this-m;
}
Matrix Matrix::operator *=(double alpha)
{
 return *this*alpha;
}
Matrix Matrix::operator *=(Matrix & m1)
{
 return *this*m1;
}
Matrix sqrt(Matrix m)
{
 assert(m[0][0]>=0 && m.getHeight()==1 && m.getWidth()==1);
 m[0][0]=sqrt(m[0][0]);
 return m;
}
double abs(Matrix & m)
{
 double sum=0;
 for(long i=0;i<m.height;i++){
  for(long j=0;j<m.width;j++){
   sum+=m[i][j]*m[i][j];
  }
 }
 return sqrt(sum);
}
double sum(Matrix & m)
{
 double total=0;
 for(long i=0;i<m.height;i++){
  for(long j=0;j<m.width;j++){
   total+=m[i][j];
  }
 }
 return total;
}
Matrix multiply(Matrix & m1,Matrix & m2){
 assert(m1.getHeight()==m2.getHeight() && m1.getWidth()==m2.getWidth());
 Matrix m(m1.getHeight(),m1.getWidth());
 for(long i=0;i<m1.getHeight();i++){
  for(long j=0;j<m1.getWidth();j++){
   m[i][j]=m1[i][j]*m2[i][j];
  }
 }
 //cout<<m;
 return m;
}
Matrix Matrix::operator /(Matrix &m1)
{
 assert(m1.width==1 && m1.height==1);
 assert(m1[0][0]!=0);
 return *this/m1[0][0];
}
Matrix Matrix::operator /(double sub)
{
 assert(sub!=0);
 Matrix m=*this;
 for(long i=0;i<height;i++){
  for(long j=0;j<width;j++){
   m[i][j]=*(p+width*i+j)/sub;
  }
 }
 return m;
}
Matrix & Matrix::operator =(Matrix & m)
{
 if(&m==this) return *this;
 height=m.height;
 width=m.width;
 delete [] p;
 p=new double[height*width];
 for(long i=0;i<height;i++){
    for(long j=0;j<width;j++){
     *(p+width*i+j)=*(m.p+width*i+j);
    }
  }
 return *this;
}
double * Matrix::operator [](long heightPos)
{
 assert(heightPos>=0 && heightPos<=height);//报错
 return p+width*heightPos;
}
double operator+(double dbl,Matrix & m)
{
 assert(m.getHeight()==1 && m.getWidth()==1);
 return dbl+m[0][0];
}
double operator+(Matrix & m,double dbl)
{
 return dbl+m;
}
double operator-(double dbl,Matrix & m)
{
 assert(m.getHeight()==1 && m.getWidth()==1);
 return dbl-m[0][0];
}
double operator-(Matrix & m,double dbl)
{
 return -(dbl-m);
}
ostream & operator<<(ostream & os,Matrix & m)
{
 os<<"Matrix:height="<<m.height<<endl;
 os<<"Matrix:width="<<m.width<<endl;
 for(long i=0;i<m.height;i++){
  for(long j=0;j<m.width;j++){
   os.width(15);
   os<<m[i][j]; 
  }
  os<<endl;
 }
 return os;
}
Matrix::~Matrix(void)
{
 delete [] p; 
 //cout<<"Default Destructor Invoked."<<endl;
}
void Matrix::Test(){
 double arr[4][4]={{1,1,1,1},{1,2,3,4},{1,3,6,10},{1,4,10,20}};
 Matrix m1,m(&arr[0][0],4,4),m2(&arr[0][0],4,4);
 m.isPtv();
 /*           
 cout<<"arranger="<<m.Arg()<<endl;
 cout<<m;*/
}
void Matrix::info(){
 cout<<"height="<<height<<endl;
 cout<<"width="<<width<<endl;
}