matrix.cpp

基于线性规划的回归支持向量机源程序

#include "matdef.hpp"
#include "matrix.hpp"

char *Merrorstring[ElastError]=
{
	"[[error in error reporting ????]]",
	"the dimension is negative",
	"the dimensions of two matrix conflict",
	"out of data range",
	"the matrix is singular",
	"iteration failure",
	"data type error",
	"the index out of define range",
	"the divisor is close to zero",
	"the method used has failed",
	"the data has no definition" ,
	"the condition is infinite",
	"the method has not been defined",
	"capacity exceeded in construction",
	"function need a square matrix"
};

//WARNING !! the sequence of next strings is in reverse order of
//			 enum Marvec_warnings !!
char *Mwarningstring[-WlastWarning]=
{
	"[[Warning in warning reporting ????]]",
	"condition is too large,the result may be invaliable",
	"norm is clear to zero, result is inaccurate",
	"matrix is near singular, inv() result is inaccurate"
};
matrix::matrix()
{
	  p=new mrep;
      if(!p)
      throw("not enough memory");
	  p->length=0;
	  p->f=NULL;
	  p->refcnt=1;
	  p->tmppointer=NULL;
}



matrix::matrix(int xsize, int ysize, double init)
{
      if((xsize<=0)||(ysize<=0))
      {
        cout<< Merrorstring[EMAT_INVALIDSIZE];
        throw Merrorstring[EMAT_INVALIDSIZE];
    }

      int tmp=UINT_MAX /sizeof(double);
      if((xsize>=floor(UINT_MAX/4))||(ysize>=tmp))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  p=new mrep;
      if(!p)
      throw("not enough memory");
	  p->length=xsize;
      p->f=new vector *[xsize];
	  for(int i=0; i<xsize; i++)
          p->f[i]=new vector(ysize,init);
	   p->refcnt=1;
	  p->tmppointer=NULL;
}

int matrix::initial(int xsize, int ysize)
{
      if((xsize<=0)||(ysize<=0))
      {
        cout<< Merrorstring[EMAT_INVALIDSIZE];
        throw Merrorstring[EMAT_INVALIDSIZE];
    }

      int tmp=UINT_MAX /sizeof(double);
      if((xsize>=floor(UINT_MAX/4))||(ysize>=tmp))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  int retn=1;
    if(p->refcnt==1)
	{
        if( (mrow()==xsize)&&(mcol()==ysize) )
        return 1;
		if(p->f !=NULL)
		{
			int len=p->length;
			for(int i=0; i<len; i++)
                 delete p->f[i];
			if(p->tmppointer !=NULL)
				delete p->tmppointer;
            delete []p->f;
		}
	}
    else
    {
        p->refcnt--;
        p=new mrep;
    }
    if(xsize && ysize)
    {
        p->length=xsize;
        p->f=new vector *[xsize];
          if(!(p->f))
          throw("not enough memory");
        for(int i=0; i<xsize; i++)
        p->f[i]=new vector(ysize);
    }

  if(!p)
  throw("not enough memory");
    p->refcnt=1;
    p->tmppointer=NULL;
	return retn;
}

int matrix::initrow(int xsize)
{
      if((xsize<=0))
      {
        cout<< Merrorstring[EMAT_INVALIDSIZE];
        throw Merrorstring[EMAT_INVALIDSIZE];
    }

      int tmp=UINT_MAX /sizeof(double);
      if(xsize>=floor(UINT_MAX/4))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  int retn=1;
    if(p->refcnt==1)
	{
        if( mrow()==xsize)
        return 1;
		if(p->f !=NULL)
		{
			int len=p->length;
			for(int i=0; i<len; i++)
				 delete p->f[i];
			if(p->tmppointer !=NULL)
				delete p->tmppointer;
			delete p->f;
		}
	}
    else
    {
        p->refcnt--;
        p=new mrep;
    }
      if(!p)
      throw("not enough memory");
    p->length=xsize;
    if( xsize )
    {
        p->f=new vector *[xsize];
        for(int i=0; i<xsize; i++)
        p->f[i]=new vector;
    }

    p->refcnt=1;
    p->tmppointer=NULL;
	return retn;
}


matrix::matrix(double  **a,int xlength,int ylength)
{
      if((xlength<=0)||(ylength<=0))
      {
        cout<< Merrorstring[EMAT_INVALIDSIZE];
        throw Merrorstring[EMAT_INVALIDSIZE];
    }

      int tmp=UINT_MAX /sizeof(double);
      if((xlength>=floor(UINT_MAX/4))||(ylength>=tmp))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  p=new mrep;
	  p->length=xlength;
      p->f=new vector *[xlength];
	  for(int i=0; i<xlength; i++)
          p->f[i]=new vector(a[i],ylength);
	   p->refcnt=1;
	  p->tmppointer=NULL;
}

//constructs a new matrix with matrix a1 and a2. the construction mode is
//controlled by variable select:
//select=0: 		  the matrix A1(n1*m1) and B1(n2*m2) will be placed at
//					  main diagnal and the dimension of new matrix is
//					  (n1+n2)*(m1+m2);
//select=1:	   	      combines matrix A1(n1*m1) and B1(n2*m2).
//					  the dimension of new matrix is n*(m1+m2)
//					  (n=MAX(n1,n2));
//select=2(default):  combines matrix A1(n1*m1) and B1(n2*m2).
//					  the dimension of new matrix is (n1+n2)*m
//					  (m=MAX(m1,m2));
//other select:		  reports an error !

matrix::matrix(matrix  & a1,
        matrix  & a2,int select)
{
	  int xlength,ylength,minc,ninc,i,j;
	  int mrow1=a1.mrow();
	  int mrow2=a2.mrow();
	  int mcol1=a1.mcol();
	  int mcol2=a2.mcol();
	  switch(select)
	  {
		case 0:
			xlength=mrow1+mrow2;
			ylength=mcol1+mcol2;
			minc=mrow1;
			ninc=mcol1;
			break;
		case 1:
			xlength=MAX(mrow1,mrow2);
			ylength=mcol1+mcol2;
			minc=0;
			ninc=mcol1;
			break;
		case 2:
			xlength=mrow1+mrow2;
			ylength=MAX(mcol1,mcol2);
			minc=mrow1;
			ninc=0;
			break;
		default:
        cout<<Merrorstring[EMAT_ASSIGNDATAERR];
        throw Merrorstring[EMAT_ASSIGNDATAERR];
        ;
	  }
      int tmp=UINT_MAX /sizeof(double);
      if((xlength>=floor(UINT_MAX/4))||(ylength>=tmp))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  p=new mrep;
	  p->length=xlength;
      p->f=new vector *[p->length];
	  for(i=0; i<xlength; i++)
          p->f[i]=new vector(ylength);
	  for(i=0; i<mrow1; i++)
		  for(j=0; j<mcol1; j++)
			 p->f[i]->p->f[j]=a1[i][j];
	  for(i=0; i<mrow2; i++)
		  for(j=0; j<mcol2; j++)
			p->f[i+minc]->p->f[j+ninc]=a2[i][j];
	  p->refcnt=1;
	  p->tmppointer=NULL;
}

//constructs an submatrix using matrix a. the dimension of new matrix is
//(xlast-xfirst)*(ylast-yfist). the first element is a[xfirst][yfirst].

matrix:: matrix(matrix  & a,int xfirst,int xlast,
						int yfirst,int ylast)
{
	  int xsize=xlast-xfirst;
	  int ysize=ylast-yfirst;
      if((xsize<=0)||(ysize<=0))
      {
        cout<<Merrorstring[EMAT_ASSIGNDATAERR];
        throw Merrorstring[EMAT_ASSIGNDATAERR];
    }
      if((xfirst<0)||(yfirst<0)||(xlast>a.mrow())||(ylast>a.mcol()))
      {
        cout<< Merrorstring[EMAT_OUTOFRANGE];
        throw Merrorstring[EMAT_OUTOFRANGE];
    }

	  p=new mrep;
	  p->length=xsize;
      p->f=new vector *[xsize];
	  for(int i=0; i<xsize; i++)
          p->f[i]=new vector(a[xfirst+i],yfirst,ylast);
	   p->refcnt=1;
	  p->tmppointer=NULL;
}



matrix::matrix(matrix  & x)
{
	x.p->refcnt++;
	p=x.p;
}



matrix::~matrix()
{
	if(--p->refcnt==0)
	{
		if(p->f !=NULL)
		{
			int len=p->length;
			for(int i=0; i<len; i++)
				 delete p->f[i];
			if(p->tmppointer !=NULL)
				delete p->tmppointer;
			delete p->f;
		}
		delete p;
	}
}


matrix copy(matrix  & m)
{
	int mr=m.mrow();
	int mc=m.mcol();
    matrix x(mr,mc);
	for(int i=0; i<mr; i++)
		 *(x.p->f[i])=copy(*(m.p->f[i]));
	return x;
}


matrix matrix::copy()
{
	int mr=mrow();
	int mc=mcol();
    matrix x(mr,mc);
	for(int i=0; i<mr; i++)
		 *(x.p->f[i])=::copy(*(p->f[i]));
	return x;
}



matrix::operator double  ** ()
{
	if(p->tmppointer==NULL)
	{
		int n=mrow();

        p->tmppointer=new double *[n];
		for(int i=0; i<n; i++)
			  p->tmppointer[i]=p->f[i]->p->f;
	}
	return p->tmppointer;
}



matrix  &  matrix::operator=( matrix  & m)
{
	m.p->refcnt++;
	if(--p->refcnt==0)
	{
		int len=p->length;
		for(int i=0; i<len; i++)
			delete p->f[i];
		delete p->f;
		if(p->tmppointer !=NULL)
			delete p->tmppointer;
		delete p;
	}
	p=m.p;
	return *this;
}


matrix  &  matrix::operator=(vector  & x)
{
	int len=p->length;
	if(p->refcnt>1)
	{                    //x has been referenced by other matrix.
		int ven=p->f[0]->p->length;
		p->refcnt--;
		p=new mrep;      //produce a new matrix data structure.
		p->length=len;
		p->refcnt=1;
        p->f=new vector*[len];
		for(int i=0; i<len; i++)
            p->f[i]=new vector(ven);
		p->tmppointer=NULL;
	}
    vector **f=p->f;
	for(int i=0; i<len; i++)
		*f[i]=x;
	return *this;
}


matrix  &  matrix::operator=(double x)
{
	int len=p->length;
	if(p->refcnt>1)
	{
		int ven=p->f[0]->p->length;
		p->refcnt--;
		p=new mrep;
		p->length=len;
		p->refcnt=1;
        p->f=new vector*[len];
		for(int i=0; i<len; i++)
            p->f[i]=new vector(ven);
		p->tmppointer=NULL;
	}
    vector **f=p->f;
	for(int i=0; i<len; i++)
		*f[i]=x;
	return *this;
}


//m=m1+m2

matrix operator+(matrix  & m1,matrix  & m2)
{
    //allocate the new matrix
    if(m1!=m2)
    {
        cout<<"matrix's dim are not equal";
    }
	int len=MIN(m1.p->length,m2.p->length);
	int ven=MIN(vlen(m1.p->f[0]),vlen(m2.p->f[0]));

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
    vector **f2=m2.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] + *f2[i];
	return m;
}

//m=m1-m2

matrix operator-(matrix  & m1,matrix  & m2)
{
    //allocate the new matrix
    if(m1!=m2)
    {
        cout<<"matrix's dim are not equal";
    }
	int len=MIN(m1.p->length,m2.p->length);
	int ven=MIN(vlen(m1.p->f[0]),vlen(m2.p->f[0]));

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
    vector **f2=m2.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] - *f2[i];
	return m;
}


//m=m1*m2

matrix operator*(matrix  & m1,matrix  & m2)
{
	//check the dimensions
	int i,j;
	int m1rows=m1.mrow();
	int m1cols=m1.mcol();
	int m2rows=m2.mrow();
	int m2cols=m2.mcol();
    if(m1cols != m2rows)
    {
        cout<<Merrorstring[EMAT_INVALIDORDER];
        throw Merrorstring[EMAT_INVALIDORDER];
    }
    matrix m(m1rows,m2cols);
    // create the matrix with its row and columns flips
    matrix flip(m2cols, m2rows);
	for(i=0; i<m2rows; i++)
		for(j=0; j<m2cols; j++)
			flip[j][i]=m2[i][j];

	//create the cross product
	for( i=0; i<m1rows; i++)
		for( j=0; j<m2cols; j++){
			m[i][j]=m1[i].dot(flip[j]);
		}
	return m;
}


//m=m1/m2

matrix operator/(matrix  & m1,matrix  & m2)
{
    //allocate the new matrix
    if(m1!=m2)
    {
        cout<<"matrix's dim are not equal";
    }
	int len=MIN(m1.p->length,m2.p->length);
	int ven=MIN(vlen(m1.p->f[0]),vlen(m2.p->f[0]));

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
    vector **f2=m2.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] / *f2[i];
	return m;
}

//m=m1+x

matrix operator+(matrix  & m1,double x)
{
    //allocate the new matrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] + x;
	return m;
}

//m=m1-x

matrix operator-(matrix  & m1,double x)
{
    //allocate the new matrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] - x;
	return m;
}

//m=m1*x

matrix operator*(matrix  & m1,double x)
{
    //allocate the new matrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] *x;
	return m;
}


//v=m1*v2

vector operator*(matrix  & m1,vector  & v2)
{
	//check the dimensions
	int m1rows=m1.mrow();
	int m1cols=m1.mcol();
	int m2rows=v2.vlen();
    if(m1cols != m2rows)
    {
        cout<<Merrorstring[EMAT_INVALIDORDER];
        throw Merrorstring[EMAT_INVALIDORDER];
    }
    vector v(m1rows);

	//create the cross product
	for( int i=0; i<m1rows; i++){
			 v[i]=m1[i].dot(v2);
	}
	return v;
}

//v=v1*m2;

vector operator*(vector  & v1,matrix  & m2)
{
	//check the dimensions
	int i,j;
	int m1cols=v1.vlen();
	int m2rows=m2.mrow();
	int m2cols=m2.mcol();
    if(m1cols != m2rows)
    {
        cout<<Merrorstring[EMAT_INVALIDORDER];
        throw Merrorstring[EMAT_INVALIDORDER];
    }
    vector v(m2cols);
    matrix flip(m2cols, m2rows);
	for(i=0; i<m2rows; i++)
		for(j=0; j<m2cols; j++)
			flip[j][i]=m2[i][j];

	//create the cross product
	for(i=0; i<m2cols; i++){
			 v[i]=v1.dot(flip[i]);
	}
	return v;
}

//m=m1/x

matrix operator/(matrix  & m1,double x)
{
    //allocate the new matrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] / x;
	return m;
}

//m=x+m1

matrix operator+(double x, matrix  & m1)
{
    //allocate the new matrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=x+ *f1[i];
	return m;
}

//m=x-m1

matrix operator-(double x, matrix  & m1)
{
    //allocate the new matrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=x- *f1[i];
	return m;
}


//m=x*m1

matrix operator*(double x, matrix  & m1)
{
    //allocate the new matrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=x* *f1[i];
	return m;
}

//m=x/m1

matrix operator/(double x, matrix  & m1)
{
    //allocate the new matrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    matrix m(len,ven);
	//do the operations
    vector **f=m.p->f;
    vector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=x/ *f1[i];
	return m;
}

//m=+m1

matrix operator+(matrix  & m1)
{
	return m1;
}

//m=-m1

matrix operator-(matrix  & m1)
{
    //allocate the new matrix