matrix__.h

很多二维 三维几何计算算法 C++ 类库

\\|M.row_dimension() == M1.row_dimension()| and
\\|M.column_dimension() == M1.column_dimension()|
}
}*/

Matrix_<NT_,AL_> operator+ (const Matrix_<NT_,AL_>& M1); 
/*{\Mbinop Addition. \precond \dimeq.}*/

Matrix_<NT_,AL_> operator- (const Matrix_<NT_,AL_>& M1); 
/*{\Mbinop Subtraction. \precond \dimeq.}*/

Matrix_<NT_,AL_> operator-(); // unary
/*{\Munop Negation.}*/

Matrix_<NT_,AL_>& operator-=(const Matrix_<NT_,AL_>&); 

Matrix_<NT_,AL_>& operator+=(const Matrix_<NT_,AL_>&); 

Matrix_<NT_,AL_> operator*(const Matrix_<NT_,AL_>& M1) const; 
/*{\Mbinop Multiplication. \precond \\ |\Mvar.column_dimension() = M1.row_dimension()|. }*/

Vector_<NT_,AL_> 
operator*(const Vector_<NT_,AL_>& vec) const
{  return ((*this) * Matrix_<NT_,AL_>(vec)).to_vector(); }
/*{\Mbinop  Multiplication with vector. \precond \\
|\Mvar.column_dimension() = vec.dimension()|.}*/

Matrix_<NT_,AL_> compmul(const NT& x) const; 

static int compare(const Matrix_<NT_,AL_>& M1, 
                   const Matrix_<NT_,AL_>& M2);

}; // end of class

/*{\Xtext \headerline{Input and Output}}*/

template <class NT_, class AL_> 
std::ostream&  operator<<(std::ostream& os, const Matrix_<NT_,AL_>& M);
/*{\Xbinopfunc writes matrix |\Mvar| row by row to the output stream |os|.}*/

template <class NT_, class AL_> 
std::istream&  operator>>(std::istream& is, Matrix_<NT_,AL_>& M);
/*{\Xbinopfunc reads matrix |\Mvar| row by row from the input stream |is|.}*/


template <class NT_, class AL_>
Matrix_<NT_,AL_>::
Matrix_(int dim) : dm_(dim),dn_(dim)
{ 
  CGAL_assertion_msg((dim >= 0), 
    "Matrix_::constructor: negative dimension.");
  if (dm_ > 0) { 
    allocate_mat_space(v_,dm_);
    for (int i=0; i<dm_; i++) 
      v_[i] = new Vector(dn_);
  } else 
    v_ = (Vector**)0; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_>::
Matrix_(int dim1, int dim2) : dm_(dim1),dn_(dim2)
{ 
  CGAL_assertion_msg((dim1>=0 && dim2>=0), 
    "Matrix_::constructor: negative dimension.");

  if (dm_ > 0) { 
    allocate_mat_space(v_,dm_);
    for (int i=0; i<dm_; i++) 
      v_[i] = new Vector(dn_); 
  } else 
    v_ = (Vector**)0; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_>::
Matrix_(std::pair<int,int> p) : dm_(p.first),dn_(p.second)
{ 
  CGAL_assertion_msg((dm_>=0 && dn_>=0), 
    "Matrix_::constructor: negative dimension.");
  if (dm_ > 0) { 
    allocate_mat_space(v_,dm_);
    for (int i=0; i<dm_; i++) 
      v_[i] = new Vector(dn_); 
  } else 
    v_ = (Vector**)0; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_>::
Matrix_(int dim, const Identity&, const NT& x) : dm_(dim),dn_(dim)
{ CGAL_assertion_msg((dim >= 0),
    "matrix::constructor: negative dimension.");
  if (dm_ > 0) { 
    allocate_mat_space(v_,dm_);
    for (int i=0; i<dm_; i++) 
      v_[i] = new Vector(dn_); 
    if (x!=NT(0)) for (int i=0; i<dm_; ++i) elem(i,i)=x;
  } else 
    v_ = (Vector**)0; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_>::
Matrix_(int dim1, int dim2, const NT& x) : dm_(dim1),dn_(dim2)
{ CGAL_assertion_msg((dim1>=0 && dim2>=0), 
    "Matrix_::constructor: negative dimension.");
  if (dm_ > 0) { 
    allocate_mat_space(v_,dm_);
    for (int i=0; i<dm_; ++i) v_[i] = new Vector(dn_,x);
  } else 
    v_ = (Vector**)0; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_>::
Matrix_(const Matrix_<NT_,AL_>& p) : dm_(p.dm_),dn_(p.dn_)
{ if (dm_ > 0) {  
    allocate_mat_space(v_,dm_);
    for (int i=0; i<dm_; i++) 
      v_[i] = new Vector(*p.v_[i]); 
  }
  else 
    v_ = (Vector_<NT_,AL_>**)0; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_>::
Matrix_(const Vector& v) : dm_(v.d_),dn_(1)
{ if (dm_>0) allocate_mat_space(v_,dm_);
  else v_ = (Vector_<NT_,AL_>**)0; 
  for(int i = 0; i < dm_; i++) { 
    v_[i] = new Vector(1); 
    elem(i,0) = v[i]; 
  }
}


template <class NT_, class AL_>
Matrix_<NT_,AL_>::
Matrix_(int dim1, int dim2, NT** p) : dm_(dim1),dn_(dim2)
{ 
  CGAL_assertion_msg((dim1 >= 0 && dim2 >= 0), 
    "Matrix_::constructor: negative dimension.");
  if (dm_ > 0) {
    allocate_mat_space(v_,dm_);
    for(int i=0; i<dm_; i++) { 
      v_[i] = new Vector_<NT_,AL_>(dn_); 
      for(int j=0; j<dn_; j++) 
        elem(i,j) = p[i][j]; 
    }
  } else 
    v_ = (Vector_<NT_,AL_>**)0; 
}


template <class NT_, class AL_>
Matrix_<NT_,AL_>& Matrix_<NT_,AL_>::
operator=(const Matrix_<NT_,AL_>& mat)
{ 
  if (&mat == this)
    return *this;

  int i,j; 
  if (dm_ != mat.dm_ || dn_ != mat.dn_) { 
    for(i=0; i<dm_; i++) delete v_[i]; 
    if (v_) deallocate_mat_space(v_,dm_);

    dm_ = mat.dm_; dn_ = mat.dn_; 
    if (dm_>0)
      allocate_mat_space(v_,dm_);
    for(i = 0; i < dm_; i++) 
      v_[i] = new Vector(dn_); 
  }

  for(i = 0; i < dm_; i++)
    for(j = 0; j < dn_; j++) 
      elem(i,j) = mat.elem(i,j); 
  return *this; 
}


template <class NT_, class AL_>
Matrix_<NT_,AL_>::
~Matrix_()  
{ 
  if (v_) {
    for (int i=0; i<dm_; i++) 
      delete v_[i];  
    deallocate_mat_space(v_,dm_);
  }
}


template <class NT_, class AL_>
inline bool Matrix_<NT_,AL_>::
operator==(const Matrix_<NT_,AL_>& x) const
{ 
  int i,j; 
  if (dm_ != x.dm_ || dn_ != x.dn_) 
    return false; 

  for(i = 0; i < dm_; i++)
    for(j = 0; j < dn_; j++)
      if (elem(i,j) != x.elem(i,j)) 
        return false; 
  return true; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_> Matrix_<NT_,AL_>::
operator+ (const Matrix_<NT_,AL_>& mat)
{ 
  int i,j; 
  check_dimensions(mat); 
  Matrix_<NT_,AL_> result(dm_,dn_); 
  for(i=0; i<dm_; i++)
    for(j=0; j<dn_; j++)
      result.elem(i,j) = elem(i,j) + mat.elem(i,j); 
  return result; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_> Matrix_<NT_,AL_>::
operator- (const Matrix_<NT_,AL_>& mat)
{ 
  int i,j; 
  check_dimensions(mat); 
  Matrix_<NT_,AL_> result(dm_,dn_); 
  for(i=0; i<dm_; i++)
    for(j=0; j<dn_; j++)
      result.elem(i,j) = elem(i,j) - mat.elem(i,j); 
  return result; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_> Matrix_<NT_,AL_>::
operator- ()  // unary
{ 
  int i,j; 
  Matrix_<NT_,AL_> result(dm_,dn_); 
  for(i=0; i<dm_; i++)
    for(j=0; j<dn_; j++)
      result.elem(i,j) = -elem(i,j); 
  return result; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_>& Matrix_<NT_,AL_>::
operator-= (const Matrix_<NT_,AL_>& mat) 
{ 
  int i,j; 
  check_dimensions(mat); 
  for(i=0; i<dm_; i++)
    for(j=0; j<dn_; j++)
      elem(i,j) -= mat.elem(i,j); 
  return *this; 
}

template <class NT_, class AL_>
Matrix_<NT_,AL_>& Matrix_<NT_,AL_>::
operator+= (const Matrix_<NT_,AL_>& mat) 
{ 
  int i,j; 
  check_dimensions(mat); 
  for(i=0; i<dm_; i++)
    for(j=0; j<dn_; j++)
      elem(i,j) += mat.elem(i,j); 
  return *this; 
}

template <class NT_, class AL_>
inline Matrix_<NT_,AL_> Matrix_<NT_,AL_>::
operator*(const Matrix_<NT_,AL_>& M1) const
{ CGAL_assertion_msg((dn_==M1.dm_), 
    "Matrix_::operator*: incompatible matrix types."); 
  Matrix_<NT_,AL_> result(dm_,M1.dn_); 
  int i,j,l;
  for (i=0; i<dm_; ++i)
    for (j=0; j<M1.dn_; ++j)
      for (l=0; l<dn_; ++l)
        result.elem(i,j) += elem(i,l)*M1.elem(l,j);
  return result;
}

template <class NT_, class AL_>
inline Matrix_<NT_,AL_> Matrix_<NT_,AL_>::
compmul(const NT& f) const
{ 
  int i,j; 
  Matrix_<NT_,AL_> result(dm_,dn_); 
  for(i=0; i<dm_; i++)
    for(j=0; j<dn_; j++)
      result.elem(i,j) = elem(i,j) *f; 
  return result; 
}


template <class NT, class AL>

Matrix_<NT,AL>  operator*(const NT& x, const Matrix_<NT,AL>& M)
/*{\Mbinopfunc Multiplication of every entry with |x|. }*/
{ return M.compmul(x); }

template <class NT, class AL>

Matrix_<NT,AL>  operator*(const Matrix_<NT,AL>& M, const NT& x)
/*{\Mbinopfunc Multiplication of every entry with |x|. }*/
{ return M.compmul(x); }



template <class NT_, class AL_> 
int Matrix_<NT_,AL_>::
compare(const Matrix_<NT_,AL_>& M1, const Matrix_<NT_,AL_>& M2) 
{ int i; int res;
  M1.check_dimensions(M2);
  for(i=0; i < M1.row_dimension() && 
      (res = compare(M1.row(i),M2.row(i))) != 0; i++);
  return res;
}


template <class NT_, class AL_> 
std::ostream&  operator<<(std::ostream& os, const Matrix_<NT_,AL_>& M)
{ 
  /* syntax: d1 d2 
             x_0,0    ... x_0,d1-1
                  d2-times
             x_d2-1,0 ... x_d2-1,d1-1 */

    int d = M.row_dimension();
    int k = M.column_dimension();
    switch (os.iword(CGAL::IO::mode)) {
    case CGAL::IO::BINARY:
        CGAL::write( os, d);
        CGAL::write( os, k);
        for ( int i = 0; i < d; ++i) {
            for ( int j = 0; j < k; ++j) {
                CGAL::write( os, M[i][j]);
            }
        }
        break;
    case CGAL::IO::ASCII:
        os << d << ' ' << k;
        for ( int i = 0; i < d; ++i) {
            for ( int j = 0; j < k; ++j) {
                os << ' ' << M[i][j];
            }
        }
        break;
    case CGAL::IO::PRETTY:
        os << "LA::Matrix((" << d << ", " << k << " [";
        for ( int i = 0; i < d; ++i) {
            for ( int j = 0; j < k; ++j) {
                if ( j != 0)
                    os << ',' << ' ';
                os << M[i][j];
            }
            if ( i != d)
                os << ",\n";
        }
        os << "])";
        break;
    }
    return os;
}

template <class NT_, class AL_> 
std::istream&  operator>>(std::istream& is, Matrix_<NT_,AL_>& M) 
{ 
  /* syntax: d1 d2 
             x_0,0  ... x_0,d1-1
                  d2-times
             x_d2,0 ... x_d2,d1-1 */

  int cdim, rdim, i;
  switch(is.iword(CGAL::IO::mode)) {
    case CGAL::IO::BINARY : 
      CGAL::read(is,rdim);
      CGAL::read(is,cdim);
      for (i=0; i<rdim*cdim; ++i)
        CGAL::read(is,M(i/rdim,i%cdim));
      break;
    case CGAL::IO::ASCII :
      is >> rdim >> cdim;
      M = Matrix_<NT_,AL_>(rdim,cdim);
      for (i=0; i<rdim*cdim; ++i)
        is >> M(i/rdim,i%cdim);
      break; 
    default:
      std::cerr<<"\nStream must be in ascii or binary mode"<<std::endl;
      break;
  }
  return is;
}

template <class NT_, class AL_>
typename Matrix_<NT_,AL_>::allocator_type Matrix_<NT_,AL_>::MM;


/*{\Ximplementation 
The data type |\Mname| is implemented by two-dimensional arrays of
variables of type |NT|. The memory layout is row oriented. Operation
|column| takes time $O(n)$, |row|, |dim1|, |dim2| take constant time,
and all other operations take time $O(nm)$.  The space requirement is
$O(nm)$.}*/


} // CGALLA
#endif // CGAL_MATRIX___H