matrix__.h

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

// Copyright (c) 1997-2000  Utrecht University (The Netherlands),
// ETH Zurich (Switzerland), Freie Universitaet Berlin (Germany),
// INRIA Sophia-Antipolis (France), Martin-Luther-University Halle-Wittenberg
// (Germany), Max-Planck-Institute Saarbruecken (Germany), RISC Linz (Austria),
// and Tel-Aviv University (Israel).  All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 of the License.
// See the file LICENSE.LGPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Kernel_d/include/CGAL/Kernel_d/Matrix__.h $
// $Id: Matrix__.h 36979 2007-03-10 10:42:31Z spion $
// 
//
// Author(s)     : Michael Seel <[e-mail]>

#ifndef CGAL_MATRIX___H
#define CGAL_MATRIX___H

#include <CGAL/Kernel_d/Vector__.h>
#include <new>
#include <cstddef>                 // for std::size_t, std::ptrdiff_t

namespace CGALLA {

template <typename ROW_, typename V_, typename R_, typename P_> 
class column_iterator_ {
  ROW_ row_; unsigned i_;
public:
  typedef column_iterator_ Self;
  typedef V_  value_type;
  typedef R_  reference;
  typedef P_  pointer;
  typedef std::size_t     size_type;
  typedef std::ptrdiff_t  difference_type;
  typedef std::random_access_iterator_tag iterator_category;

  column_iterator_() : row_(),i_() {}
  column_iterator_(ROW_ row, unsigned i) : row_(row),i_(i) {}

  bool  operator==( const Self& x) const 
  { return row_ == x.row_ && i_ == x.i_; }
  bool  operator!=( const Self& x) const 
  { return !(*this == x); }

  R_ operator*()  const { return (**row_)[i_]; }
  P_ operator->() const { return (**row_)+i_; }

  Self& operator++() { ++row_; return *this; }
  Self  operator++(int) 
  { Self tmp = *this; ++*this; return tmp; }

  Self& operator--() { --row_; return *this; }
  Self  operator--(int) 
  { Self tmp = *this; --*this; return tmp; }

  Self operator+(difference_type i) const
  { return Self(row_+i,i_); }
  Self operator-(difference_type i) const
  { return Self(row_-i,i_); }
  difference_type operator-(const Self& x) const
  { return (row_ - x.row_); }
      
};

template <typename ROW_, typename V_, typename R_, typename P_> 
class component_iterator_ {
  ROW_ row_;  // pointer to row
  int i_, n_; // offset and limit
public:
  typedef component_iterator_ Self;
  typedef V_  value_type;
  typedef R_  reference;
  typedef P_  pointer;
  typedef std::size_t     size_type;
  typedef std::ptrdiff_t  difference_type;
  typedef std::bidirectional_iterator_tag iterator_category;

  component_iterator_() : row_(),i_(),n_() {}
  component_iterator_(ROW_ row, int i, int n) 
    : row_(row),i_(i),n_(n) {}

  bool  operator==( const Self& x) const 
  { return row_==x.row_ && i_==x.i_; }
  bool  operator!=( const Self& x) const   
  { return !(*this == x); }

  R_    operator*()  const { return (**row_)[i_]; }
  P_    operator->() const { return (**row_)+i_; }

  Self& operator++() { ++i_; if (i_==n_) { ++row_; i_=0; } return *this; }
  Self  operator++(int) { Self tmp = *this; ++*this; return tmp; }
  Self& operator--() { --i_; if (i_<0) { --row_; i_=n_-1; } return *this; }
  Self  operator--(int) { Self tmp = *this; --*this; return tmp; }
     
};


/*{\Msubst
<NT_,AL_>#
<NT,AL>#
Vector_#Vector
Matrix_#Matrix
}*/
/*{\Moptions print_title=yes}*/
/*{\Moptions outfile=Matrix.man}*/
/*{\Manpage {Matrix}{}{Matrices with NT Entries}{M}}*/

template <class NT_, class AL_>
class Matrix_ 
{ 
/*{\Mdefinition An instance of data type |\Mname| is a matrix of
variables of number type |NT|. The types |\Mname| and |Vector_|
together realize many functions of basic linear algebra.}*/

public:

/*{\Mtypes 6}*/

typedef Vector_<NT_,AL_>* vector_pointer;
typedef const Vector_<NT_,AL_>* const_vector_pointer;

typedef NT_ NT;
/*{\Mtypemember the ring type of the components.}*/ 

typedef component_iterator_<vector_pointer*,NT,NT&,NT*> iterator;
/*{\Mtypemember bidirectional iterator for accessing all components
row-wise.}*/
typedef component_iterator_<vector_pointer*,NT,const NT&,const NT*> 
  const_iterator;

typedef NT* row_iterator;
/*{\Mtypemember random access iterator for accessing row
  entries.}*/
typedef const NT* row_const_iterator;

typedef column_iterator_<vector_pointer*,NT,NT&,NT*> column_iterator;
/*{\Mtypemember random access iterator for accessing column entries.}*/ 
typedef column_iterator_<vector_pointer*,NT,const NT&, const NT*> 
  column_const_iterator;

/*{\Mtext There also constant versions of the above iterators:
|const_iterator|, |row_const_iterator|, and |column_const_iterator|.}*/

class Identity {};
/*{\Mtypemember a tag class for identity initialization}*/

typedef Vector_<NT_,AL_> Vector;
/*{\Mtypemember the vector type used.}*/ 

protected:
vector_pointer* v_; int dm_,dn_; 

NT& elem(int i, int j) const { return v_[i]->v_[j]; }

typedef typename AL_::template rebind<vector_pointer>::other 
        allocator_type;
static allocator_type MM;

inline void allocate_mat_space(vector_pointer*& vi, int d)
{
  /* We use this procedure to allocate memory. We use our allocator
     memory allocation scheme. There we first get an appropriate piece
     of memory and then initialize each cell by an inplace new. */

  vi = MM.allocate(d); 
  vector_pointer* p = vi + d - 1; 
  while (p >= vi) { 
    new (p) vector_pointer*(0); p--;
  }
}

inline void deallocate_mat_space(vector_pointer*& vi, int d)
{
  /* deallocate memory via our AL_ object. */

  MM.deallocate(vi,d);
  vi = (vector_pointer*)0;
}

inline void check_dimensions(const Matrix_<NT_,AL_>& mat) const
{ 
  CGAL_assertion_msg((dm_ == mat.dm_ && dn_ == mat.dn_), 
    "Matrix::check_dimensions: incompatible matrix dimensions.") ;
}

public:

/*{\Mcreation 5}*/

Matrix_() : dm_(0),dn_(0) { v_ = (Vector**)0; }
/*{\Mcreate creates an instance |\Mvar| of type |\Mname|.}*/
Matrix_(int n); 
/*{\Mcreate creates an instance |\Mvar| of type |\Mname| of
dimension $n \times n$ initialized to the zero matrix.}*/
Matrix_(int m, int n); 
/*{\Mcreate creates an instance |\Mvar| of type |\Mname| of 
dimension $m \times n$ initialized to the zero matrix.}*/
Matrix_(std::pair<int,int> p);
/*{\Mcreate creates an instance |\Mvar| of type |\Mname| of dimension
|p.first|$\times$|p.second| initialized to the zero matrix.}*/
Matrix_(int n , const Identity&, const NT& x = NT(1) ); 
/*{\Mcreate creates an instance |\Mvar| of type |\Mname| of
dimension $n \times n$ initialized to the identity matrix
(times |x|).}*/
Matrix_(int m, int n, const NT& x);
/*{\Mcreate creates an instance |\Mvar| of type |\Mname| of 
dimension $m \times n$ initialized to the matrix with |x|
entries.}*/

template <class RAIterator>
void range_initialize(RAIterator first, RAIterator last,  
                      std::random_access_iterator_tag) 
{ typedef typename std::iterator_traits<RAIterator>::value_type value_type;
  typedef typename value_type::const_iterator const_iterator;
  dn_ = last-first;
  if (dn_ == 0) { dm_=0; v_=0; return; }
  dm_ = first->dimension(); 
  if (dm_ > 0) { 
    int i,j;
    allocate_mat_space(v_,dm_);
    for (i=0; i<dm_; i++) {
      v_[i] = new Vector(dn_);
      // for (int j = 0; j < dn_; j++) elem(i,j) = (*(first+j))[i];
    }
    const_iterator it;
    for (j=0; first != last; ++j, ++first) // column wise
      for (i=0, it=first->begin(); it != first->end(); ++i, ++it) // row wise
        elem(i,j) = *it;
  } else 
    v_ = (Vector**)0; 
}

template <class InputIterator>
void range_initialize(InputIterator first, InputIterator last, 
                 std::forward_iterator_tag) 
{ typedef typename std::iterator_traits<InputIterator>::value_type 
    value_type;
  std::vector<value_type> V(first,last);
  range_initialize(V.begin(),V.end(),std::random_access_iterator_tag());
}

template <class Forward_iterator>
Matrix_(Forward_iterator first, Forward_iterator last)
/*{\Mcreate creates an instance |\Mvar| of type |\Mname|. Let $S$ be
the ordered set of $n$ column-vectors of common dimension $m$ as given
by the iterator range |[first,last)|.  |\Mvar| is initialized to an $m
\times n$ matrix with the columns as specified by $S$.  \precond
|Forward_iterator| has a value type |V| from which we require to
provide a iterator type |V::const_iterator|, to have |V::value_type ==
NT|.\\ Note that |Vector_| or |std::vector<NT>| fulfill these
requirements.}*/
{ typedef typename std::iterator_traits<Forward_iterator>::iterator_category 
    iterator_category;
  range_initialize(first,last,iterator_category()); }

Matrix_(const std::vector< Vector >& A) 
/*{\Mcreate creates an instance |\Mvar| of type |\Mname|. Let $A$ be
an array of $n$ column-vectors of common dimension $m$.  |\Mvar| is
initialized to an $m \times n$ matrix with the columns as specified by
$A$. }*/
{ range_initialize(A.begin(),A.end(),
    std::random_access_iterator_tag()); }

Matrix_(const Matrix_<NT_,AL_>&); 

Matrix_(const Vector&); 
/* creates a $d \times 1$ matrix */

Matrix_(int, int, NT**); 

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

~Matrix_(); 

/*{\Moperations 3 4}*/

int row_dimension()  const {  return dm_; }
/*{\Mop returns $n$, the number of rows of |\Mvar|.}*/

int column_dimension()  const { return dn_; }
/*{\Mop returns $m$, the number of columns of |\Mvar|.}*/

std::pair<int,int> dimension() const 
/*{\Mop returns $(m,n)$, the dimension pair of |\Mvar|.}*/
{ return std::pair<int,int>(dm_,dn_); }

Vector& row(int i) const
/*{\Mop returns the $i$-th row of |\Mvar| (an $m$ - vector).\\
\precond  $0 \le i \le m - 1$. }*/
{ CGAL_assertion_msg((0<=i && i<dm_),"Matrix_: row index out of range.");
  return *v_[i]; 
}

Vector column(int i) const 
/*{\Mop returns the $i$-th column of |\Mvar| (an $n$ - vector).\\
\precond  $0 \le i \le n - 1$. }*/
{ return Vector(column_begin(i),column_end(i)); }

Vector to_vector() const 
{ 
  CGAL_assertion_msg((dn_==1), 
    "Matrix_::to_vector: cannot make vector from matrix.");
  return column(0); 
}

Vector_<NT_,AL_>& operator[](int i) const  
{ 
  CGAL_assertion_msg((0<=i && i<dm_), 
    "Matrix_::operator[]: index out of range.");
  return row(i); 
}

NT& operator()(int i, int j)
/*{\Mfunop returns $M_{ i,j }$. \\
\precond $0\le i\le m-1$ and $0\le j\le n-1$. }*/
{ CGAL_assertion_msg((0<=i && i<dm_), 
    "Matrix_::operator(): row index out of range.");
  CGAL_assertion_msg((0<=j && j<dn_), 
    "Matrix_::operator(): column index out of range.");
  return elem(i,j); 
}

NT operator()(int i, int j) const
{ 
  CGAL_assertion_msg((0<=i && i<dm_), 
    "Matrix_::operator(): row index out of range.");
  CGAL_assertion_msg((0<=j && j<dn_), 
    "Matrix_::operator(): column index out of range.");
  return elem(i,j); 
}

void swap_rows(int i, int j)
/*{\Mop swaps rows $i$ and $j$.
\precond $0\le i\le m-1$ and $0\le j\le m-1$.}*/
{ CGAL_assertion(0<=i && i<dm_ && 0<=j && j<dm_);
  std::swap(v_[i],v_[j]); 
}

void swap_columns(int i, int j) 
/*{\Mop swaps columns $i$ and $j$.
\precond $0\le i\le n-1$ and $0\le j\le n-1$.}*/
{ CGAL_assertion(0<=i && i<dn_ && 0<=j && j<dn_);
  for(int l = 0; l < dm_; l++) std::swap(elem(l,i),elem(l,j)); 
}

row_iterator row_begin(int i) 
/*{\Mop an iterator pointing to the first entry of the $i$th row.
\precond $0\le i\le m-1$.}*/
{ CGAL_assertion_msg((0<=i&&i<dm_),"Matrix: row index out of range.");
  return v_[i]->begin(); }

row_iterator row_end(int i)   
/*{\Mop an iterator pointing beyond the last entry of the $i$th row.
\precond $0\le i\le m-1$.}*/
{ CGAL_assertion_msg((0<=i&&i<dm_),"Matrix: row index out of range.");
  return v_[i]->end(); }

row_const_iterator row_begin(int i) const 
{ CGAL_assertion_msg(0<=i&&i<dm_,"Matrix: row index out of range.");
  return v_[i]->begin(); }
row_const_iterator row_end(int i) const  
{ CGAL_assertion_msg(0<=i&&i<dm_,"Matrix: row index out of range.");
  return v_[i]->end(); }

column_iterator column_begin(int i) 
/*{\Mop an iterator pointing to the first entry of the $i$th column.
\precond $0\le i\le n-1$.}*/
{ CGAL_assertion_msg(0<=i&&i<dn_,"Matrix: column index out of range.");
  return column_iterator(v_,i); }

column_iterator column_end(int i)   
/*{\Mop an iterator pointing beyond the last entry of the $i$th column.
\precond $0\le i\le n-1$.}*/
{ return column_begin(i)+dm_; }

column_const_iterator column_begin(int i) const 
{ CGAL_assertion_msg(0<=i&&i<dn_,"Matrix: column index out of range.");
  return column_const_iterator(v_,i); }
column_const_iterator column_end(int i)   const 
{ return column_begin(i)+dm_; }

iterator begin() { return iterator(v_,0,dn_); }
/*{\Mop an iterator pointing to the first entry of |\Mvar|.}*/
iterator end() { return iterator(v_+dm_,0,dn_); }
/*{\Mop an iterator pointing beyond the last entry of |\Mvar|.}*/

const_iterator begin() const { return const_iterator(v_,0,dn_); }
const_iterator end() const { return const_iterator(v_+dm_,0,dn_); }

/*{\Mtext The same operations exist for |row_const_iterator| and
|column_const_iterator|.}*/

bool  operator==(const Matrix_<NT_,AL_>& M1)  const; 
/*{\Mbinop Test for equality. }*/

bool  operator!=(const Matrix_<NT_,AL_>& M1)  const 
/*{\Mbinop Test for inequality. }*/
{ return !(*this == M1); }

/*{\Mtext \headerline{Arithmetic Operators}}*/
/*{\Mtext
\settowidth{\typewidth}{|Matrix_<NT,LA>m|}
\addtolength{\typewidth}{\colsep}
\callwidth2cm
\computewidths
\newcommand{\dimeq}[2]{