matrix.cpp

qrdrls算法

void MATRIX::subtract_from_col(int _col, complex2 _val)
{
  if (_col < 0 || _col >= col)
    {
      matrix_err = ILLEGAL_INDEX;
#ifdef NO_RCS
      printf("MATRIX:subtract_from_col: ERROR --- %d!!!\n", matrix_err);
      printf("Invalid column index. Cant subtract from column ...\n");
#else
      rcs_print("MATRIX:subtract_from_col: ERROR --- %d!!!\n", matrix_err);
      rcs_print("Invalid column index. Cant subtruct from column  ...\n");
#endif
    }
  else
    {
      for (int ii = 0; ii < row; ii++)
	{
	  elements[ii][_col] -= _val;
	}
    }
}

void MATRIX::multiply_col(int _col, complex2 _val)
{
  if (_col < 0 || _col >= col)
    {
      matrix_err = ILLEGAL_INDEX;
#ifdef NO_RCS
      printf("MATRIX:multiply_col: ERROR --- %d!!!\n", matrix_err);
      printf("Invalid column index. Cant multiply column ...\n");
#else
      rcs_print("MATRIX:multiply_col: ERROR --- %d!!!\n", matrix_err);
      rcs_print("Invalid column index. Cant multiply column  ...\n");
#endif
    }
  else
    {
      for (int ii = 0; ii < row; ii++)
	{
	  elements[ii][_col] *= _val;
	}
    }
}

void MATRIX::devide_col(int _col, complex2 _val)
{
  if (_col < 0 || _col >= col)
    {
      matrix_err = ILLEGAL_INDEX;
#ifdef NO_RCS
      printf("MATRIX:divide_col: ERROR --- %d!!!\n", matrix_err);
      printf("Invalid column index. Cant divide column ...\n");
#else
      rcs_print("MATRIX:divide_col: ERROR --- %d!!!\n", matrix_err);
      rcs_print("Invalid column index. Cant divide column  ...\n");
#endif
    }
  else
    {
      if (!_val.IsZero())
	{
	  for (int ii = 0; ii < row; ii++)
	    {
	      elements[ii][_col] *= _val;
		  
	    }
	}
      else
	{
	  matrix_err = DIVIDE_BY_ZERO;
#ifdef NO_RCS
	  printf("MATRIX:divide_col: ERROR --- %d!!!\n", matrix_err);
	  printf("Divide by 0. Cant divide column ...\n");
#else
	  rcs_print("MATRIX:divide_col: ERROR --- %d!!!\n", matrix_err);
	  rcs_print("Divide by 0. Cant divide column ...\n");
#endif
	}
    }
}

void MATRIX::add_rows(int _source_row, int _destination_row, complex2 gain)
{
  if ((_source_row < 0 || _source_row >= row) || 
      (_destination_row < 0 || _destination_row >= row))
    {
      matrix_err = ILLEGAL_INDEX;
#ifdef NO_RCS
      printf("MATRIX:add_rows: ERROR --- %d!!!\n", matrix_err);
      printf("Invalid index. Cant add rows ...\n");
#else
      rcs_print("MATRIX:add_rows: ERROR --- %d!!!\n", matrix_err);
      rcs_print("Invalid index. Cant add rows  ...\n");
#endif
    }
  else
    {
      for (int ii = 0; ii < col; ii++)
	{
	  elements[_destination_row][ii] += gain*elements[_source_row][ii];
	}
    }
}

void MATRIX::add_cols(int _source_col, int _destination_col, complex2 gain)
{
  if ((_source_col < 0 || _source_col >= col) || 
      (_destination_col < 0 || _destination_col >= col))
    {
      matrix_err = ILLEGAL_INDEX;
#ifdef NO_RCS
      printf("MATRIX:add_cols: ERROR --- %d!!!\n", matrix_err);
      printf("Invalid index. Cant add columns ...\n");
#else
      rcs_print("MATRIX:add_cols: ERROR --- %d!!!\n", matrix_err);
      rcs_print("Invalid index. Cant add columns  ...\n");
#endif
    }
  else
    {
      for (int ii = 0; ii < row; ii++)
	{
	  elements[ii][_destination_col] += gain*elements[ii][_source_col];
	}
    }
}

MATRIX inv(const MATRIX &_matrix)
{
  int ii, jj; // Index variales for the loops

  // Get the # of rows and columns of the input matrix
  int r = _matrix.get_num_rows();
  int c = _matrix.get_num_cols();

  // Temporary variables
  complex2 temp;
  MATRIX ans(r, c);
  MATRIX copy = _matrix;
  
  // Check whether square or not
  if (r != c)  
    { // Nonsquare matrix
      ans.matrix_err = NONSQUARE_MATRIX;
#ifdef NO_RCS
      printf("inv: ERROR --- Nonsqure matrix. Cant invert!!!\n");
      printf("Use pinv for the pseudoinverse ...\n");
#else
      rcs_print("inv: ERROR --- Nonsqure matrix. Cant invert!!!\n");
      rcs_print("Use pinv for the pseudoinverse ...\n");
#endif
      // Return matrix of zeros
      ans.set_matrix(complex2(0));
      return (ans);
    }
  else // Square
    { // Inverts using row operations. Lame, but works most times.
      
      // The answer is identity initially
      ans.set_diagonal(1);
      
      for (ii = 0; ii < r; ii++) 
	{
	  // Get the next diagonal element
	  temp = copy.get_element(ii, ii);

	  if (temp == 0)
	    {
	      ans.matrix_err = NONINVERTIBLE_MATRIX;
#ifdef NO_RCS
	      printf("inv: ERROR --- Noninvertible matrix. Cant invert!!!\n");
	      printf("Check whether your matrix is singular ...\n");
#else
	      rcs_print("inv: ERROR --- Noninvertible matrix. Cant invert!!!\n");
	      rcs_print("Check whether your matrix is singular ...\n");
#endif
	      // Return matrix of zeros
	      ans.set_matrix(complex2(0));
	      return (ans);
	    }
	  else 
	    {
	      ans.multiply_row(ii, complex2(1)/temp);
	      copy.multiply_row(ii, complex2(1)/temp);
	    }
	  for (jj = 0; jj < r; jj++)
	    if (ii != jj) 
	      {
		ans.add_rows(ii, jj, complex2(0)-copy.get_element(jj, ii));
		copy.add_rows(ii, jj, complex2(0)-copy.get_element(jj, ii));
	      }
	}
      return (ans);
    }
}

MATRIX pinv(const MATRIX &_matrix)
{
  // Get the # of rows and columns of the input matrix
  int r = _matrix.get_num_rows();
  int c = _matrix.get_num_cols();
  
  // Inverts using row operations. Lame, but works most times.
  if (r == c)
    {
      return(inv(_matrix));
    }
  else
    {
      if (r < c)
	{ // Return the right pseudoinverse
	  return ( transpose(_matrix) * inv(_matrix * transpose(_matrix)) );
	}        
      else
	{ // Return the left pseudoinverse
	  return ( inv(transpose(_matrix) * _matrix) * transpose(_matrix) );
	}        
    }
}

MATRIX transpose(const MATRIX &_matrix)
{
  // Get the # of rows and columns of the input matrix
  int r = _matrix.get_num_rows();
  int c = _matrix.get_num_cols();
  
  int ii, jj; // Index variales for the loops

  MATRIX ans(c,r);

  for (ii = 0; ii < r; ii++)
    {
      for (jj = 0; jj < c; jj++)
	{
	  ans.set_element(jj, ii, _matrix.get_element(ii, jj));
	}
    }
  return (ans);
}

MATRIX transpose_col(const MATRIX &_matrix)

{
  // Get the # of rows and columns of the input matrix
  int r = _matrix.get_num_rows();
  int c = _matrix.get_num_cols();
  
  int ii, jj; // Index variales for the loops

  MATRIX ans(r,c);

  for (ii = 0; ii < r; ii++)
    {
      for (jj = 0; jj < c; jj++)
	{
	  ans.set_element(ii,jj, _matrix.get_element(ii, c-1-jj));
	}
    }
  return (ans);
}



MATRIX hermite(const MATRIX &_matrix)
{
  // Get the # of rows and columns of the input matrix
  int r = _matrix.get_num_rows();
  int c = _matrix.get_num_cols();

  int ii, jj; // Index variales for the loops

  MATRIX temp(c,r);  
  MATRIX ans(c,r);

  for (ii = 0; ii < r; ii++)
    {
      for (jj = 0; jj < c; jj++)
	  { 
        temp.set_element(ii,jj,_matrix.get_element(ii, jj).conjugate());
		ans=transpose(temp);
	}
    }
  return (ans);
}

complex2 trace(const MATRIX &_matrix)
{
  // Get the # of rows and columns of the input matrix
  int r = _matrix.get_num_rows();
  int c = _matrix.get_num_cols();

  int ii; // Index variales for the loops
  complex2 ans = 0;  // The variabe to hold the answer

  // Find the smaller of the rows and columns
  int s = (r < c ? r : c);

  for (ii = 0; ii < s; ii++)
    {
      ans += _matrix.get_element(ii, ii);
    }
  return (ans);
}

complex2 norm(const MATRIX &_matrix)
{
  return ( trace(_matrix * transpose(_matrix)) );
}

MATRIX &MATRIX::operator=(const MATRIX &_matrix2)
{
  int ii, jj; // Index variales for the loops

  if (&_matrix2 != this)
  {
    free_memory();
    row = _matrix2.row;
    col = _matrix2.col;
    allocate_memory();
    for (ii = 0; ii < row; ii++) 
    {
      for (jj = 0; jj < col; jj++)
      {
	elements[ii][jj] = _matrix2.elements[ii][jj];
      }
    }
  }
  return (*this);
}

const MATRIX &MATRIX::operator+=(const MATRIX &_matrix2)
{
  int ii, jj; // Index variales for the loops

  if ((row != _matrix2.row) || (col != _matrix2.col))
    {
      matrix_err = SIZE_MISMATCH;
    }
  else
    {
      for (ii = 0; ii < row; ii++)
	{
	  for (jj = 0; jj < col; jj++)
	    {
	      elements[ii][jj] += _matrix2.elements[ii][jj];
	    }
	}
    }
  return (*this);
}


const MATRIX &MATRIX::operator-=(const MATRIX &_matrix2)
{
  int ii, jj; // Index variales for the loops

  if ((row != _matrix2.row) || (col != _matrix2.col))
    {
      matrix_err = SIZE_MISMATCH;
    }
  else
    {
      for (ii = 0; ii < row; ii++)
	{
	  for (jj = 0; jj < col; jj++)
	    {
	      elements[ii][jj] -= _matrix2.elements[ii][jj];
	    }
	}
    }
  return (*this);
}


const MATRIX &MATRIX::operator*=(const complex2 _val)
{
  int ii, jj; // Index variales for the loops

  for (ii = 0; ii < row; ii++)
    {
      for (jj = 0; jj < col; jj++)
	{
	  elements[ii][jj] = elements[ii][jj]*_val;
	}
    }
  return (*this);
}

MATRIX MATRIX::operator+(const MATRIX &_matrix2) const
{
  int ii, jj; // Index variales for the loops

  MATRIX ans = *this;

  if ((row != _matrix2.row) || (col != _matrix2.col))
    {
      ans.matrix_err = SIZE_MISMATCH;
    }
  else
    {
      for (ii = 0; ii < row; ii++)
	{
	  for (jj = 0; jj < col; jj++)
	    {
	      ans.elements[ii][jj] += _matrix2.elements[ii][jj];
	    }
	}
    }
  return (ans);
}

MATRIX MATRIX::operator-(const MATRIX &_matrix2) const
{
  int ii, jj; // Index variales for the loops

  MATRIX ans = *this;

  if ((row != _matrix2.row) || (col != _matrix2.col))
    {
      ans.matrix_err = SIZE_MISMATCH;
    }
  else
    {
      for (ii = 0; ii < row; ii++)
	{
	  for (jj = 0; jj < col; jj++)
	    {
	      ans.elements[ii][jj] -= _matrix2.elements[ii][jj];
	    }
	}
    }
  return (ans);
}


MATRIX MATRIX::operator*(const MATRIX &_matrix2) const
{
  int ii, jj, kk; // Index variales for the loops

  MATRIX ans(row, _matrix2.col);

  if (col != _matrix2.row)
    {
      ans.matrix_err = SIZE_MISMATCH;
    }
  else
    {
      for (ii = 0; ii < row; ii++)
	{
	  for (jj = 0; jj < _matrix2.col; jj++)
	    {
	      for (kk = 0; kk < col; kk++)
		{
		  ans.elements[ii][jj] += 
		    elements[ii][kk] * _matrix2.elements[kk][jj];
		}
	    }
	}
    }
  return (ans);
}



MATRIX  MATRIX::operator*(const complex2 _val)
{
  int ii, jj; // Index variales for the loops

  MATRIX ans(row,col);

  for (ii = 0; ii < row; ii++)
    {
      for (jj = 0; jj < col; jj++)
	{
	  ans.elements[ii][jj] = elements[ii][jj]*_val;
	}
    }
  return (ans);
}

/************************************************************/


MATRIX operator*(const complex2 _val, const MATRIX &_matrix)
{
  int ii, jj; // Index variales for the loops
 
  MATRIX ans(_matrix.row,_matrix.col);

  for (ii = 0; ii < _matrix.row; ii++)
    {
      for (jj = 0; jj < _matrix.col; jj++)
	{
	  ans.elements[ii][jj] = _matrix.elements[ii][jj]*_val;
	}
    }
  return (ans);
}



void MATRIX::print_matrix()
{   
	int i,j;
	//MATRIX ans(_matrix.row,_matrix.col);
	for( i=0; i<row;i++)
	{
		for( j=0;j<col;j++)
		{
			print(elements[i][j]);
			cout<<"\t";
		}

		//	{
			 // ans.elements[i][j] = _matrix.elements[i][j];
			//  cout<<" ans.elements[i][j]";
		//	}
		  cout<<"\n";
		}
	cout<<endl;
}


 
/**********************************************************************/

/*void main()
{ 
  complex a00=complex(1,1);
  complex a01=complex(1,2);
  complex a10=complex(2,2);
  complex a11=complex(2,3);
  complex b00=complex(4,5);
  complex b01=complex(2,4);
  complex b10=complex(3,4);
  complex b11=complex(4,2);
  complex val=complex(1,1);

  complex row0[2]={a00,a01};
  complex row1[2]={a10,a11};
  complex row2[2]={b00,b01};
  complex row3[2]={b10,b11};

  MATRIX R(2, 2);
  MATRIX T(2, 2);
  R.set_row(0,row0);
  R.set_row(1,row1);
  T.set_row(0,row2);
  T.set_row(1,row3);
  MATRIX A=R+T;
  MATRIX S=R-T;
  MATRIX M=R*T;
  MATRIX V=T*val;
  MATRIX I,B,C,Q;
  complex2 D,E;
  
  I= inv( T );
  B=transpose(T);
  C=pinv(T);
  D=trace(T);
  E=norm(T);
  Q=hermite(T);

  cout<<"\nR=\n";
  R.print_matrix();
  cout<<"\nT=\n";
  T.print_matrix();
  cout<<"\nA=R+T=\n";
  A.print_matrix();
  cout<<"\nS=R-T=\n";
  S.print_matrix();
  cout<<"\nM=R*T=\n";
  M.print_matrix();
  cout<<"\nV=T*val=\n";
  V.print_matrix();
  cout<<"\nI=InvT=\n";
  I.print_matrix();
  cout<<"\nB=transposeT=\n";
  B.print_matrix();
  cout<<"\nC=pinvT=\n";
  C.print_matrix();
  cout<<"\nD=traceT=\n";
  print(D);
  cout<<"\nE=normT=\n";
  print(E);
  cout<<"\nQ=hermiteT=\n";
  Q.print_matrix();
}*/