From 2e2fe30fbd6b60c2dbbec4f3229d6cee698022cc Mon Sep 17 00:00:00 2001 From: Michael Seel Date: Mon, 19 Feb 2001 19:34:10 +0000 Subject: [PATCH] added Herves Cartesian linear algebra and debugged it --- .../Kernel_d/include/CGAL/Homogeneous_d.h | 4 +- .../CGAL/Kernel_d/Aff_transformation_d.h | 48 + .../include/CGAL/Kernel_d/Direction_d.h | 40 + .../include/CGAL/Kernel_d/Hyperplane_d.h | 56 + .../Kernel_d/include/CGAL/Kernel_d/Imatrix.h | 422 +++--- .../Kernel_d/include/CGAL/Kernel_d/Ivector.h | 270 ++-- .../include/CGAL/Kernel_d/Linear_algebraCd.C | 468 ++++++ .../include/CGAL/Kernel_d/Linear_algebraHd.C | 1036 +++++++++++++ .../Kernel_d/include/CGAL/Kernel_d/Vector_d.h | 77 + .../Kernel_d/include/CGAL/Kernel_d/d_utils.h | 67 + .../Kernel_d/include/CGAL/Linear_algebra.h | 1323 ----------------- .../Kernel_d/include/CGAL/Linear_algebraCd.h | 122 ++ .../Kernel_d/include/CGAL/Linear_algebraHd.h | 213 +++ .../test/Kernel_d/Linear_algebra-test.C | 479 ++++-- .../test/Kernel_d/Linear_algebra-test.cmd | 1 + 15 files changed, 2840 insertions(+), 1786 deletions(-) create mode 100644 Packages/Kernel_d/include/CGAL/Kernel_d/Aff_transformation_d.h create mode 100644 Packages/Kernel_d/include/CGAL/Kernel_d/Direction_d.h create mode 100644 Packages/Kernel_d/include/CGAL/Kernel_d/Hyperplane_d.h create mode 100644 Packages/Kernel_d/include/CGAL/Kernel_d/Linear_algebraCd.C create mode 100644 Packages/Kernel_d/include/CGAL/Kernel_d/Linear_algebraHd.C create mode 100644 Packages/Kernel_d/include/CGAL/Kernel_d/Vector_d.h create mode 100644 Packages/Kernel_d/include/CGAL/Kernel_d/d_utils.h delete mode 100644 Packages/Kernel_d/include/CGAL/Linear_algebra.h create mode 100644 Packages/Kernel_d/include/CGAL/Linear_algebraCd.h create mode 100644 Packages/Kernel_d/include/CGAL/Linear_algebraHd.h create mode 100644 Packages/Kernel_d/test/Kernel_d/Linear_algebra-test.cmd diff --git a/Packages/Kernel_d/include/CGAL/Homogeneous_d.h b/Packages/Kernel_d/include/CGAL/Homogeneous_d.h index 782bbee37ee..61fdca269c5 100644 --- a/Packages/Kernel_d/include/CGAL/Homogeneous_d.h +++ b/Packages/Kernel_d/include/CGAL/Homogeneous_d.h @@ -8,7 +8,7 @@ #include #include #endif -#include +#include #include #include #include @@ -34,7 +34,7 @@ template class Ray_d; template class Line_d; template class Aff_transformation_d; -template > +template > class Homogeneous_d { public: diff --git a/Packages/Kernel_d/include/CGAL/Kernel_d/Aff_transformation_d.h b/Packages/Kernel_d/include/CGAL/Kernel_d/Aff_transformation_d.h new file mode 100644 index 00000000000..49bfea608c6 --- /dev/null +++ b/Packages/Kernel_d/include/CGAL/Kernel_d/Aff_transformation_d.h @@ -0,0 +1,48 @@ +#ifndef CGAL_AFF_TRANSFORMATION_D_H +#define CGAL_AFF_TRANSFORMATION_D_H + +CGAL_BEGIN_NAMESPACE + +template +class Aff_transformation_d : public pR::Aff_transformation_d_base +{ public: + typedef typename pR::Aff_transformation_d_base Base; + typedef Aff_transformation_d Self; + typedef pR R; + typedef typename R::RT RT; + typedef typename R::FT FT; + typedef typename R::LA LA; + + Aff_transformation_d(int d = 0, bool identity=false) + : Base(d,identity) {} + Aff_transformation_d(const typename LA::Matrix& M) + : Base(M) {} + template + Aff_transformation_d( + Forward_iterator start, Forward_iterator end) : Base(start,end) {} + Aff_transformation_d(const Vector_d& v) : Base(v) {} + Aff_transformation_d(int d, const RT& num, const RT& den) + : Base(d,num,den) {} + Aff_transformation_d(int d, const RT& sin_num, const RT& cos_num, + const RT& den, int e1 = 0, int e2 = 1) + : Base(d,sin_num,cos_num,den,e1,e2) {} + Aff_transformation_d(int d, const Direction_d& dir, + const RT& num, const RT& den, + int e1 = 0, int e2 = 1) + : Base(d,dir,num,den,e1,e2) {} + Aff_transformation_d(const Base& a) : Base(a) {} + Aff_transformation_d(const Self& a) : Base(a) {} + + Self operator*(const Self& a) + { return Base::operator*(a); } + Self inverse() const { return Base::inverse(); } + + bool operator==(const Self& a) const + { return Base::operator==(a); } + bool operator!=(const Self& a) const + { return Base::operator!=(a); } + +}; + +CGAL_END_NAMESPACE +#endif //CGAL_AFF_TRANSFORMATION_D_H diff --git a/Packages/Kernel_d/include/CGAL/Kernel_d/Direction_d.h b/Packages/Kernel_d/include/CGAL/Kernel_d/Direction_d.h new file mode 100644 index 00000000000..4802cdf3707 --- /dev/null +++ b/Packages/Kernel_d/include/CGAL/Kernel_d/Direction_d.h @@ -0,0 +1,40 @@ +#ifndef CGAL_DIRECTION_D_H +#define CGAL_DIRECTION_D_H + +CGAL_BEGIN_NAMESPACE + +template +class Direction_d : public pR::Direction_d_base +{ public: + typedef typename pR::Direction_d_base Base; + typedef Direction_d Self; + typedef pR R; + typedef typename R::RT RT; + typedef typename R::FT FT; + typedef typename R::LA LA; + + Direction_d(int d=0) : Base(d) {} + Direction_d(int a, int b) : Base(a,b) {} + Direction_d(const RT& a, const RT& b) : Base(a,b) {} + Direction_d(int a, int b, int c) : Base(a,b,c) {} + Direction_d(const RT& a, const RT& b, const RT& c) : Base(a,b,c) {} + template + Direction_d (int d, InputIterator first, InputIterator last) + : Base(d, first, last) {} + Direction_d(const Direction_d &d) : Base(d) {} + Direction_d(const Vector_d &v) : Base(v) {} + Direction_d(typename Base::Base_direction, int d, int i) : + Base(typename Base::Base_direction(),d,i) {} + Direction_d(const Base& p) : Base(p) {} + + Self operator-() const { return Base::operator-(); } + Vector_d vector() const { return Base::vector(); } + + bool operator==(const Self& w) const + { return Base::operator==(w); } + bool operator!=(const Self& w) const + { return Base::operator!=(w); } +}; + +CGAL_END_NAMESPACE +#endif //CGAL_DIRECTION_D_H diff --git a/Packages/Kernel_d/include/CGAL/Kernel_d/Hyperplane_d.h b/Packages/Kernel_d/include/CGAL/Kernel_d/Hyperplane_d.h new file mode 100644 index 00000000000..e52b65f04fc --- /dev/null +++ b/Packages/Kernel_d/include/CGAL/Kernel_d/Hyperplane_d.h @@ -0,0 +1,56 @@ +#ifndef CGAL_HYPERPLANE_D_H +#define CGAL_HYPERPLANE_D_H + +CGAL_BEGIN_NAMESPACE + +template +class Hyperplane_d : public pR::Hyperplane_d_base +{ public: + typedef typename pR::Hyperplane_d_base Base; + typedef Hyperplane_d Self; + typedef pR R; + typedef typename R::RT RT; + typedef typename R::FT FT; + typedef typename R::LA LA; + + Hyperplane_d(int d=0) : Base(d) {} + Hyperplane_d(int a, int b, int c) : Base(a,b,c) {} + Hyperplane_d(const RT& a, const RT& b, const RT& c) : + Base(a,b,c) {} + Hyperplane_d(int a, int b, int c, int d) : Base(a,b,c,d) {} + Hyperplane_d(const RT& a, const RT& b, const RT& c, const RT& d) : + Base(a,b,c,d) {} + + Hyperplane_d(const Point_d& p, const Direction_d& dir) : + Base(p,dir) {} + + Hyperplane_d(const Hyperplane_d &h) : Base(h) {} + Hyperplane_d(const Base& p) : Base(p) {} + + template + Hyperplane_d(int d, InputIterator first, InputIterator last) + : Base (d, first, last) {} + + template + Hyperplane_d(int d, InputIterator first, InputIterator last, + const RT& D) + : Base (d, first, last, D) {} + + template + Hyperplane_d(ForwardIterator first, ForwardIterator last, + const Point_d& o, Oriented_side side = Oriented_side(0)) : + Base(first,last,o,side) {} + + Vector_d orthogonal_vector() const + { return Base::orthogonal_vector(); } + Direction_d orthogonal_direction() const + { return Base::orthogonal_direction(); } + + bool operator==(const Self& w) const + { return Base::operator==(w); } + bool operator!=(const Self& w) const + { return Base::operator!=(w); } +}; + +CGAL_END_NAMESPACE +#endif //CGAL_HYPERPLANE_D_H diff --git a/Packages/Kernel_d/include/CGAL/Kernel_d/Imatrix.h b/Packages/Kernel_d/include/CGAL/Kernel_d/Imatrix.h index 4f9cfcca93c..877e0f6d223 100644 --- a/Packages/Kernel_d/include/CGAL/Kernel_d/Imatrix.h +++ b/Packages/Kernel_d/include/CGAL/Kernel_d/Imatrix.h @@ -33,48 +33,46 @@ // debugging and templatization: M. Seel //--------------------------------------------------------------------- -#ifndef LINALG_IMATRIX_H -#define LINALG_IMATRIX_H +#ifndef CGAL_IMATRIX_H +#define CGAL_IMATRIX_H #include #include +CGAL_BEGIN_NAMESPACE + /*{\Msubst -<># -<_RT,_ALLOC># +<_NT,_ALLOC># +Ivector#Vector +Imatrix#Matrix }*/ +/*{\Moptions print_title=yes}*/ +/*{\Moptions outfile=Imatrix.man}*/ +/*{\Manpage {Matrix}{NT}{Matrices with NT Entries}{M}}*/ -LA_BEGIN_NAMESPACE - -/*{\Moptions -outfile=Imatrix.man -}*/ -/*{\Manpage {Imatrix} {RT} {Matrices with RT Entries} {M}}*/ - -template +template class Imatrix { -/*{\Mdefinition -An instance of data type |Imatrix| is a matrix of variables of type -|RT|, the so-called ring type. The arithmetic type |RT| is required to -behave like integers in the mathematical sense. +/*{\Mdefinition An instance of data type |\Mname| is a matrix of +variables of number type |NT|. The types |\Mname| and |Ivector| +together realize many functions of basic linear algebra. All +functions on integer matrices compute the exact result, i.e., there is +no rounding error. |\Mname| offers all simple matrix operations. The +more sophisticated ones are provided by the class |Linear_algebra|. +Preconditions are checked by default and can be switched off by the +compile flag [[CGAL_LA_PRECOND_OFF]].}*/ -The types |Imatrix| and |Ivector| together realize many functions of -basic linear algebra. All functions on integer matrices compute the -exact result, i.e., there is no rounding error. |\Mname| offers all -simple matrix operations. The more sophisticated ones are provided by -the class |Linear_algebra|. Preconditions are checked by default and -can be switched off by the compile flag [[LA_PRECOND_OFF]].}*/ public: +public: /*{\Mtypes 5.5}*/ -typedef Ivector<_RT,_ALLOC>* vector_pointer; -typedef const Ivector<_RT,_ALLOC>* const_vector_pointer; +typedef Ivector<_NT,_ALLOC>* vector_pointer; +typedef const Ivector<_NT,_ALLOC>* const_vector_pointer; -typedef _RT RT; +typedef _NT NT; /*{\Mtypemember the ring type of the components.}*/ -typedef Ivector<_RT,_ALLOC> Vector; +typedef Ivector<_NT,_ALLOC> Vector; /*{\Mtypemember the vector type used.}*/ typedef vector_pointer iterator; @@ -83,6 +81,10 @@ typedef vector_pointer iterator; typedef const_vector_pointer const_iterator; /*{\Mtypemember the const iterator type for accessing rows.}*/ +typedef NT* Row_iterator; +/*{\Xtypemember accessing row entries.}*/ +typedef const NT& Row_const_iterator; + class Identity {}; /*{\Mtypemember a tag class for identity initialization}*/ @@ -94,7 +96,7 @@ protected: int d1; int d2; - RT& elem(int i, int j) const { return v[i]->v[j]; } + NT& elem(int i, int j) const { return v[i]->v[j]; } #ifndef _MSC_VER typedef typename _ALLOC::template rebind::other allocator_type; @@ -133,9 +135,9 @@ protected: #endif - inline void check_dimensions(const Imatrix<_RT,_ALLOC>& mat) const + inline void check_dimensions(const Imatrix<_NT,_ALLOC>& mat) const { - LA_PRECOND((d1 == mat.d1 && d2 == mat.d2), + CGAL_LA_PRECOND((d1 == mat.d1 && d2 == mat.d2), "Imatrix::check_dimensions: incompatible matrix dimensions.") } @@ -153,11 +155,12 @@ dimension $n \times n$ initialized to the zero matrix.}*/ Imatrix(int n, int m); /*{\Mcreate creates an instance |\Mvar| of type |\Mname| of dimension $n \times m$ initialized to the zero matrix.}*/ -Imatrix(int n , const Identity&, const RT& x = RT(1) ); +Imatrix(std::pair p); +Imatrix(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|).}*/ -Imatrix(int n, int m, const Initialize&, const RT& x); +Imatrix(int n, int m, const Initialize&, const NT& x); /*{\Mcreate creates an instance |\Mvar| of type |\Mname| of dimension $n \times m$ initialized to the matrix with |x| entries.}*/ @@ -187,7 +190,8 @@ void range_initialize(RAIterator first, RAIterator last, template void range_initialize(InputIterator first, InputIterator last, std::forward_iterator_tag) -{ typedef typename std::iterator_traits::value_type value_type; +{ typedef typename std::iterator_traits::value_type + value_type; std::vector V(first,last); range_initialize(V.begin(),V.end(),std::random_access_iterator_tag()); } @@ -200,7 +204,7 @@ by the iterator range |[first,last)|. |\Mvar| is initialized to an $n \times m$ 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 == -RT|.\\ Note that |Ivector| or |std::vector| fulfill these +NT|.\\ Note that |Ivector| or |std::vector| fulfill these requirements.}*/ { range_initialize(first,last, std::iterator_traits::iterator_category()); } @@ -210,82 +214,110 @@ Imatrix(const std::vector< Vector >& A) be an array of $m$ column-vectors of common dimension $n$. |\Mvar| is initialized to an $n \times m$ matrix with the columns as specified by $A$. }*/ -{ range_initialize(A.begin(),A.end(),std::random_access_iterator_tag()); } +{ range_initialize(A.begin(),A.end(), + std::random_access_iterator_tag()); } -Imatrix(const Imatrix<_RT,_ALLOC>&); +Imatrix(const Imatrix<_NT,_ALLOC>&); -Imatrix(const Ivector<_RT,_ALLOC>&); +Imatrix(const Vector&); /* creates a $d \times 1$ matrix */ -Imatrix(int, int, RT**); -Imatrix<_RT,_ALLOC>& operator=(const Imatrix<_RT,_ALLOC>&); + +Imatrix(int, int, NT**); + +Imatrix<_NT,_ALLOC>& operator=(const Imatrix<_NT,_ALLOC>&); + ~Imatrix(); /*{\Moperations 1.7 3.5}*/ int row_dimension() const { return d1; } -/*{\Mop returns $n$, the number of rows of |\Mvar|. }*/ +/*{\Mop returns $n$, the number of rows of |\Mvar|.}*/ int column_dimension() const { return d2; } -/*{\Mop returns $m$, the number of columns of |\Mvar|. }*/ +/*{\Mop returns $m$, the number of columns of |\Mvar|.}*/ + +std::pair dimension() const +/*{\Mop returns $(m,n)$, the dimension pair of |\Mvar|.}*/ +{ return std::pair(d1,d2); } Vector& row(int i) const /*{\Mop returns the $i$-th row of |\Mvar| (an $m$ - vector).\\ \precond $0 \le i \le n - 1$. }*/ { - LA_PRECOND((0<=i && i to_vector() const + +Vector to_vector() const { - LA_PRECOND((d2==1), + CGAL_LA_PRECOND((d2==1), "Imatrix::to_vector: cannot make vector from matrix.") return column(0); } -friend Ivector<_RT,_ALLOC> to_vector(const Imatrix& M) +friend Ivector<_NT,_ALLOC> to_vector(const Imatrix& M) { return M.to_vector(); } -Ivector<_RT,_ALLOC>& operator[](int i) const +Ivector<_NT,_ALLOC>& operator[](int i) const { - LA_PRECOND((0<=i && i& M1) const; +void swap_columns(int i, int j) +/*{\Mop swaps columns $i$ and $j$.}*/ +{ CGAL_assertion(0<=i && ibegin(); } +Row_iterator row_end(int i) { return v[i]->end(); } + +Row_const_iterator row_begin(int i) const { return v[i]->begin(); } +Row_const_iterator row_end(int i) const { return v[i]->end(); } + + +bool operator==(const Imatrix<_NT,_ALLOC>& M1) const; /*{\Mbinop Test for equality. }*/ -bool operator!=(const Imatrix<_RT,_ALLOC>& M1) const +bool operator!=(const Imatrix<_NT,_ALLOC>& M1) const /*{\Mbinop Test for inequality. }*/ { return !(*this == M1); } /*{\Mtext \headerline{Arithmetic Operators}}*/ /*{\Mtext -\settowidth{\typewidth}{|Imatrixm|} +\settowidth{\typewidth}{|Imatrixm|} \addtolength{\typewidth}{\colsep} \callwidth2cm \computewidths @@ -295,60 +327,60 @@ bool operator!=(const Imatrix<_RT,_ALLOC>& M1) const } }*/ -Imatrix<_RT,_ALLOC> operator+ (const Imatrix<_RT,_ALLOC>& M1); +Imatrix<_NT,_ALLOC> operator+ (const Imatrix<_NT,_ALLOC>& M1); /*{\Mbinop Addition. \precond \dimeq.}*/ -Imatrix<_RT,_ALLOC> operator- (const Imatrix<_RT,_ALLOC>& M1); +Imatrix<_NT,_ALLOC> operator- (const Imatrix<_NT,_ALLOC>& M1); /*{\Mbinop Subtraction. \precond \dimeq.}*/ -Imatrix<_RT,_ALLOC> operator-(); // unary +Imatrix<_NT,_ALLOC> operator-(); // unary /*{\Munop Negation.}*/ -Imatrix<_RT,_ALLOC>& operator-=(const Imatrix<_RT,_ALLOC>&); +Imatrix<_NT,_ALLOC>& operator-=(const Imatrix<_NT,_ALLOC>&); -Imatrix<_RT,_ALLOC>& operator+=(const Imatrix<_RT,_ALLOC>&); +Imatrix<_NT,_ALLOC>& operator+=(const Imatrix<_NT,_ALLOC>&); -Imatrix<_RT,_ALLOC> -operator*(const Imatrix<_RT,_ALLOC>& M1) const; +Imatrix<_NT,_ALLOC> +operator*(const Imatrix<_NT,_ALLOC>& M1) const; /*{\Mbinop Multiplication. \precond \\ |\Mvar.column_dimension() = M1.row_dimension()|. }*/ -Ivector<_RT,_ALLOC> -operator*(const Ivector<_RT,_ALLOC>& vec) const -{ return ((*this) * Imatrix<_RT,_ALLOC>(vec)).to_vector(); } +Ivector<_NT,_ALLOC> +operator*(const Ivector<_NT,_ALLOC>& vec) const +{ return ((*this) * Imatrix<_NT,_ALLOC>(vec)).to_vector(); } /*{\Mbinop Multiplication with vector. \precond \\ - |\Mvar.column_dimension() = vec.dimension()|.}*/ +|\Mvar.column_dimension() = vec.dimension()|.}*/ -Imatrix<_RT,_ALLOC> compmul(const RT& x) const; +Imatrix<_NT,_ALLOC> compmul(const NT& x) const; -static int compare(const Imatrix<_RT,_ALLOC>& M1, - const Imatrix<_RT,_ALLOC>& M2); +static int compare(const Imatrix<_NT,_ALLOC>& M1, + const Imatrix<_NT,_ALLOC>& M2); }; // end of class /*{\Xtext \headerline{Input and Output}}*/ -template -std::ostream& operator<<(std::ostream& os, const Imatrix& M); +template +std::ostream& operator<<(std::ostream& os, const Imatrix& M); /*{\Xbinopfunc writes matrix |\Mvar| row by row to the output stream |os|.}*/ -template -std::istream& operator>>(std::istream& is, Imatrix& M); +template +std::istream& operator>>(std::istream& is, Imatrix& M); /*{\Xbinopfunc reads matrix |\Mvar| row by row from the input stream |is|.}*/ /*{\Mimplementation The data type |\Mname| is implemented by two-dimensional arrays of -variables of type |RT|. The memory layout is row oriented. Operation +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)$.}*/ -template -Imatrix:: +template +Imatrix:: Imatrix(int dim) : d1(dim),d2(dim) { - LA_PRECOND((dim >= 0), "Imatrix::constructor: negative dimension.") + CGAL_LA_PRECOND((dim >= 0), "Imatrix::constructor: negative dimension.") if (d1 > 0) { allocate_mat_space(v,d1); for (int i=0; i -Imatrix:: +template +Imatrix:: Imatrix(int dim1, int dim2) : d1(dim1),d2(dim2) { - LA_PRECOND((dim1>=0 && dim2>=0), + CGAL_LA_PRECOND((dim1>=0 && dim2>=0), "Imatrix::constructor: negative dimension.") if (d1 > 0) { @@ -372,26 +404,41 @@ Imatrix(int dim1, int dim2) : d1(dim1),d2(dim2) v = (Vector**)0; } -template -Imatrix:: -Imatrix(int dim, const Identity&, const RT& x) : d1(dim),d2(dim) +template +Imatrix:: +Imatrix(std::pair p) : d1(p.first),d2(p.second) { - LA_PRECOND((dim >= 0), "Imatrix::constructor: negative dimension.") + CGAL_LA_PRECOND((d1>=0 && d2>=0), + "Imatrix::constructor: negative dimension.") if (d1 > 0) { allocate_mat_space(v,d1); for (int i=0; i -Imatrix:: -Imatrix(int dim1, int dim2, const Initialize&, const RT& x) : +template +Imatrix:: +Imatrix(int dim, const Identity&, const NT& x) : d1(dim),d2(dim) +{ + CGAL_LA_PRECOND((dim >= 0), + "matrix::constructor: negative dimension.") + if (d1 > 0) { + allocate_mat_space(v,d1); + for (int i=0; i +Imatrix:: +Imatrix(int dim1, int dim2, const Initialize&, const NT& x) : d1(dim1),d2(dim2) { - LA_PRECOND((dim1>=0 && dim2>=0), + CGAL_LA_PRECOND((dim1>=0 && dim2>=0), "Imatrix::constructor: negative dimension.") if (d1 > 0) { @@ -402,10 +449,9 @@ Imatrix(int dim1, int dim2, const Initialize&, const RT& x) : v = (Vector**)0; } - -template -Imatrix:: -Imatrix(const Imatrix& p) : d1(p.d1),d2(p.d2) +template +Imatrix:: +Imatrix(const Imatrix& p) : d1(p.d1),d2(p.d2) { if (d1 > 0) { allocate_mat_space(v,d1); @@ -413,43 +459,45 @@ Imatrix(const Imatrix& p) : d1(p.d1),d2(p.d2) v[i] = new Vector(*p.v[i]); } else - v = (Ivector**)0; + v = (Ivector**)0; } -template -Imatrix:: -Imatrix(const Ivector& vec) : d1(vec.dim),d2(1) +template +Imatrix:: +Imatrix(const Vector& vec) : d1(vec.dim),d2(1) { if (d1>0) allocate_mat_space(v,d1); else - v = (Ivector**)0; + v = (Ivector**)0; for(int i = 0; i < d1; i++) { v[i] = new Vector(1); elem(i,0) = vec[i]; } } -template -Imatrix:: -Imatrix(int dim1, int dim2, RT** p) : d1(dim1),d2(dim2) + +template +Imatrix:: +Imatrix(int dim1, int dim2, NT** p) : d1(dim1),d2(dim2) { - LA_PRECOND((dim1 >= 0 && dim2 >= 0), + CGAL_LA_PRECOND((dim1 >= 0 && dim2 >= 0), "Imatrix::constructor: negative dimension.") if (d1 > 0) { allocate_mat_space(v,d1); for(int i=0; i(d2); + v[i] = new Ivector(d2); for(int j=0; j**)0; + v = (Ivector**)0; } -template -Imatrix& Imatrix:: -operator=(const Imatrix& mat) + +template +Imatrix& Imatrix:: +operator=(const Imatrix& mat) { register int i,j; if (d1 != mat.d1 || d2 != mat.d2) { @@ -470,8 +518,8 @@ operator=(const Imatrix& mat) } -template -Imatrix:: +template +Imatrix:: ~Imatrix() { if (v) { @@ -481,13 +529,11 @@ Imatrix:: } } - - -template -Ivector Imatrix:: +template +Ivector Imatrix:: column(int i) const { - LA_PRECOND((i>=0 && i=0 && i -bool Imatrix:: -operator==(const Imatrix& x) const +template +bool Imatrix:: +operator==(const Imatrix& x) const { register int i,j; if (d1 != x.d1 || d2 != x.d2) @@ -511,47 +557,47 @@ operator==(const Imatrix& x) const return true; } -template -Imatrix Imatrix:: -operator+ (const Imatrix& mat) +template +Imatrix Imatrix:: +operator+ (const Imatrix& mat) { register int i,j; check_dimensions(mat); - Imatrix result(d1,d2); + Imatrix result(d1,d2); for(i=0; i -Imatrix Imatrix:: -operator- (const Imatrix& mat) +template +Imatrix Imatrix:: +operator- (const Imatrix& mat) { register int i,j; check_dimensions(mat); - Imatrix result(d1,d2); + Imatrix result(d1,d2); for(i=0; i -Imatrix Imatrix:: +template +Imatrix Imatrix:: operator- () // unary { register int i,j; - Imatrix result(d1,d2); + Imatrix result(d1,d2); for(i=0; i -Imatrix& Imatrix:: -operator-= (const Imatrix& mat) +template +Imatrix& Imatrix:: +operator-= (const Imatrix& mat) { register int i,j; check_dimensions(mat); @@ -561,9 +607,9 @@ operator-= (const Imatrix& mat) return *this; } -template -Imatrix& Imatrix:: -operator+= (const Imatrix& mat) +template +Imatrix& Imatrix:: +operator+= (const Imatrix& mat) { register int i,j; check_dimensions(mat); @@ -573,14 +619,14 @@ operator+= (const Imatrix& mat) return *this; } -template -Imatrix Imatrix:: -operator*(const Imatrix& mat) const +template +Imatrix Imatrix:: +operator*(const Imatrix& mat) const { - LA_PRECOND((d2==mat.d1), + CGAL_LA_PRECOND((d2==mat.d1), "Imatrix::operator*: incompatible matrix types.") - Imatrix result(d1, mat.d2); + Imatrix result(d1, mat.d2); register int i,j; for (i=0; i& mat) const return result; } -template -Imatrix Imatrix:: -compmul(const RT& f) const +template +Imatrix Imatrix:: +compmul(const NT& f) const { register int i,j; - Imatrix result(d1,d2); + Imatrix result(d1,d2); for(i=0; i -Imatrix operator*(const Imatrix& M, const RT& x) +template +Imatrix operator*(const Imatrix& M, const NT& x) { return M.compmul(x); } /*{\Mbinopfunc Multiplication of every entry with |x|. }*/ -template -Imatrix operator*(const Imatrix& M, int x) +template +Imatrix operator*(const Imatrix& M, int x) { return M.compmul(x); } -template -Imatrix operator*(const RT& x, const Imatrix& M) +template +Imatrix operator*(const NT& x, const Imatrix& M) { return M.compmul(x); } /*{\Mbinopfunc Multiplication of every entry with |x|. }*/ -template -Imatrix operator*(int x, const Imatrix& M) +template +Imatrix operator*(int x, const Imatrix& M) { return M.compmul(x); } -template -int Imatrix:: -compare(const Imatrix& M1, const Imatrix& M2) +template +int Imatrix:: +compare(const Imatrix& M1, const Imatrix& M2) { register int i; int res; @@ -634,49 +680,63 @@ compare(const Imatrix& M1, const Imatrix& M2) } -template -std::ostream& operator<<(std::ostream& O, const Imatrix& M) +template +std::ostream& operator<<(std::ostream& os, const Imatrix& M) { /* syntax: d1 d2 x_0,0 ... x_0,d1-1 d2-times x_d2,0 ... x_d2,d1-1 */ - O << M.row_dimension() << ' ' << M.column_dimension() << std::endl; - for (register int i=0; i < M.row_dimension(); i++) { - for (register int j=0; j < M.column_dimension(); j++) - O << M(i,j) << " "; - O << std::endl; + + CGAL::print_d prt(&os); + if (os.iword(IO::mode)==IO::PRETTY) os << "LA::Matrix("; + prt(M.row_dimension()); + prt(M.column_dimension()); + if (os.iword(IO::mode)==IO::PRETTY) { os << " [\n"; prt.reset(); } + for (register int i=0; i -std::istream& operator>>(std::istream& in, Imatrix& M) +template +std::istream& operator>>(std::istream& is, Imatrix& M) { /* syntax: d1 d2 x_0,0 ... x_0,d1-1 d2-times x_d2,0 ... x_d2,d1-1 */ - int dim1 = 0, dim2 = 0; - if (!(in >> dim1 >> dim2)) - return in; - if (M.row_dimension() != dim1 || M.column_dimension() != dim2) - M = Imatrix(dim1,dim2); - - for (register int i=0; i> M(i,j))) - return in; - return in; + int cdim, rdim; + switch(is.iword(IO::mode)) { + case IO::BINARY : + CGAL::read(is,rdim); + CGAL::read(is,cdim); + for (register int i=0; i> rdim >> cdim; + M = Imatrix(rdim,cdim); + for (register int i=0; i> M(i/rdim,i%cdim); + break; + default: + std::cerr<<"\nStream must be in ascii or binary mode"< -Imatrix::allocator_type Imatrix::MM; +template +Imatrix::allocator_type Imatrix::MM; #endif -LA_END_NAMESPACE -#endif // LINALG_IMATRIX_H +CGAL_END_NAMESPACE +#endif // CGAL_IMATRIX_H diff --git a/Packages/Kernel_d/include/CGAL/Kernel_d/Ivector.h b/Packages/Kernel_d/include/CGAL/Kernel_d/Ivector.h index ed24a35ba54..e22def5787a 100644 --- a/Packages/Kernel_d/include/CGAL/Kernel_d/Ivector.h +++ b/Packages/Kernel_d/include/CGAL/Kernel_d/Ivector.h @@ -33,43 +33,36 @@ // debugging and templatization: M. Seel //--------------------------------------------------------------------- -#ifndef LINALG_IVECTOR_H -#define LINALG_IVECTOR_H +#ifndef CGAL_IVECTOR_H +#define CGAL_IVECTOR_H #include #include +#include + #include #include #include #include #include -#define LA_BEGIN_NAMESPACE namespace CGAL { -#define LA_END_NAMESPACE } // namespace CGAL +CGAL_BEGIN_NAMESPACE -LA_BEGIN_NAMESPACE - -template +template class Imatrix; -template +template class Ivector; -#if 0 -#define LA_PRECOND(cond,s)\ -if (!(cond)) {\ - std::cout << "precondition " << #cond << " failed" << std::endl;\ - std::cout << s << std::endl; exit(1); } -#define ERROR_HANDLER(n,s)(std::cerr << s << std::endl, exit(n)) -#else -#define LA_PRECOND(cond,s) CGAL_assertion_msg((cond),s); +#define CGAL_LA_PRECOND(cond,s) CGAL_assertion_msg((cond),s); #define ERROR_HANDLER(n,s) CGAL_assertion_msg(!(n),s) -#endif - + /*{\Msubst <># -<_RT,_ALLOC># +<_NT,_ALLOC># +Ivector#Vector +Imatrix#Matrix }*/ - +/*{\Moptions print_title=yes}*/ /*{\Moptions outfile=Ivector.man}*/ /*{\Mtext \headerline{Common Notation} @@ -84,25 +77,23 @@ iterators to input iterators. If we index the tuple as above then we require that $|++|^{(d)}|first == last|$ (note that |last| points beyond the last element to be accepted).}*/ -/*{\Manpage {Ivector}{RT}{Vectors with RT Entries}{v}}*/ +/*{\Manpage {Vector}{NT}{Vectors with NT Entries}{v}}*/ -template +template class Ivector { - /*{\Mdefinition An instance of data type |Ivector| is a vector of -variables of type |RT|, the so-called ring type. Together with the -type |Imatrix| it realizes the basic operations of linear -algebra. Internal correctness tests are executed if compiled with the -flag [[CGAL_LA_SELFTEST]].}*/ +variables of number type |NT|. Together with the type |Imatrix| it +realizes the basic operations of linear algebra. Internal correctness +tests are executed if compiled with the flag [[CGAL_LA_SELFTEST]].}*/ public: /*{\Mtypes 5.5}*/ -typedef _RT* pointer; -typedef const _RT* const_pointer; +typedef _NT* pointer; +typedef const _NT* const_pointer; -typedef _RT RT; +typedef _NT NT; /*{\Mtypemember the ring type of the components.}*/ typedef pointer iterator; @@ -115,46 +106,43 @@ class Initialize {}; /*{\Mtypemember a tag class for homogeneous initialization}*/ protected: - - friend class Imatrix<_RT,_ALLOC>; - RT* v; + friend class Imatrix<_NT,_ALLOC>; + NT* v; int dim; - typedef _ALLOC allocator_type; static allocator_type MM; - inline void allocate_vec_space(RT*& vi, int di) + inline void allocate_vec_space(NT*& vi, int di) { /* We use this procedure to allocate memory. We first get an appropriate piece of memory from the allocator and then initialize each cell by an inplace new. */ vi = MM.allocate(di); - RT* p = vi + di - 1; - while (p >= vi) { new (p) RT(0); p--; } + NT* p = vi + di - 1; + while (p >= vi) { new (p) NT(0); p--; } } - inline void deallocate_vec_space(RT*& vi, int di) + inline void deallocate_vec_space(NT*& vi, int di) { /* We use this procedure to deallocate memory. We have to free it by - the allocator scheme. We first call the destructor for type RT for each + the allocator scheme. We first call the destructor for type NT for each cell of the array and then return the piece of memory to the memory manager. */ - RT* p = vi + di - 1; - while (p >= vi) { p->~RT(); p--; } + NT* p = vi + di - 1; + while (p >= vi) { p->~NT(); p--; } MM.deallocate(vi, di); - vi = (RT*)0; + vi = (NT*)0; } inline void - check_dimensions(const Ivector<_RT,_ALLOC>& vec) const + check_dimensions(const Ivector<_NT,_ALLOC>& vec) const { - LA_PRECOND((dim == vec.dim), "Ivector::check_dimensions:\ + CGAL_LA_PRECOND((dim == vec.dim), "Ivector::check_dimensions:\ object dimensions disagree.") } - public: /*{\Mcreation v 3}*/ @@ -163,22 +151,23 @@ Ivector(int d = 0) /*{\Mcreate creates an instance |\Mvar| of type |\Mname|. |\Mvar| is initialized to a vector of dimension $d$.}*/ { - LA_PRECOND( d >= 0 , "Ivector::constructor: negative dimension.") + CGAL_LA_PRECOND( d >= 0 , + "Ivector::constructor: negative dimension.") dim = d; - v = (RT*)0; + v = (NT*)0; if (dim > 0){ allocate_vec_space(v,dim); - while (d--) v[d] = RT(0); + while (d--) v[d] = NT(0); } } -Ivector(int d, const Initialize&, const RT& x) +Ivector(int d, const Initialize&, const NT& x) /*{\Mcreate creates an instance |\Mvar| of type |\Mname|. |\Mvar| is initialized to a vector of dimension $d$ with entries |x|.}*/ { - LA_PRECOND( d >= 0 , "Ivector::constructor: negative dimension.") + CGAL_LA_PRECOND( d >= 0 , "Ivector::constructor: negative dimension.") dim = d; - v = (RT*)0; + v = (NT*)0; if (dim > 0){ allocate_vec_space(v,dim); while (d--) v[d] = x; @@ -189,7 +178,7 @@ template Ivector(Forward_iterator first, Forward_iterator last) /*{\Mcreate creates an instance |\Mvar| of type |\Mname|; |\Mvar| is initialized to the vector with entries -|set [first,last)|. \precond |Forward_iterator| has value type |RT|.}*/ +|set [first,last)|. \precond |Forward_iterator| has value type |NT|.}*/ { dim = std::distance(first,last); allocate_vec_space(v,dim); iterator it = begin(); @@ -197,22 +186,22 @@ Ivector(Forward_iterator first, Forward_iterator last) } private: -void init(int d, const RT& x0, const RT& x1, const RT& x2=0, const RT& x3=0) +void init(int d, const NT& x0, const NT& x1, const NT& x2=0, const NT& x3=0) { dim = d; allocate_vec_space(v,dim); v[0]=x0; v[1]=x1; ( d>2 ? (v[2]=x2) : 0); ( d>3 ? (v[3]=x3) : 0); } public: -Ivector(const RT& a, const RT& b) { init(2,a,b); } +Ivector(const NT& a, const NT& b) { init(2,a,b); } /*{\Mcreate creates an instance |\Mvar| of type |\Mname|. |\Mvar| is initialized to the two-dimensional vector $(a,b)$.}*/ -Ivector(const RT& a, const RT& b, const RT& c) { init(3,a,b,c); } +Ivector(const NT& a, const NT& b, const NT& c) { init(3,a,b,c); } /*{\Mcreate creates an instance |\Mvar| of type |\Mname|. |\Mvar| is initialized to the three-dimensional vector $(a,b,c)$.}*/ -Ivector(const RT& a, const RT& b, const RT& c, const RT& d) +Ivector(const NT& a, const NT& b, const NT& c, const NT& d) /*{\Mcreate creates an instance |\Mvar| of type |\Mname|; |\Mvar| is initialized to the four-dimensional vector $(a,b,c,d)$.}*/ @@ -223,15 +212,15 @@ Ivector(int a, int b, int c) { init(3,a,b,c); } Ivector(int a, int b, int c, int d) { init(4,a,b,c,d); } -Ivector(const Ivector<_RT,_ALLOC>& p) +Ivector(const Ivector<_NT,_ALLOC>& p) { dim = p.dim; if (dim > 0) allocate_vec_space(v,dim); - else v = (RT*)0; + else v = (NT*)0; for(int i=0; i& operator=(const Ivector<_RT,_ALLOC>& vec) +Ivector<_NT,_ALLOC>& operator=(const Ivector<_NT,_ALLOC>& vec) { register int n = vec.dim; if (n != dim) { @@ -239,7 +228,7 @@ Ivector<_RT,_ALLOC>& operator=(const Ivector<_RT,_ALLOC>& vec) dim=n; } if (n > 0) allocate_vec_space(v,n); - else v = (RT*)0; + else v = (NT*)0; while (n--) v[n] = vec.v[n]; return *this; @@ -250,47 +239,51 @@ Ivector<_RT,_ALLOC>& operator=(const Ivector<_RT,_ALLOC>& vec) /*{\Moperations 3 3}*/ - int dimension() const { return dim; } /*{\Mop returns the dimension of |\Mvar|.}*/ + +bool is_zero() const +/*{\Mop returns true iff |\Mvar| is the zero vector.}*/ +{ for(int i=0; i& operator+=(const Ivector<_RT,_ALLOC>& v1); +Ivector<_NT,_ALLOC>& operator+=(const Ivector<_NT,_ALLOC>& v1); /*{\Mbinop Addition plus assignment. \precond\\ |v.dimension() == v1.dimension()|.}*/ -Ivector<_RT,_ALLOC>& operator-=(const Ivector<_RT,_ALLOC>& v1); +Ivector<_NT,_ALLOC>& operator-=(const Ivector<_NT,_ALLOC>& v1); /*{\Mbinop Subtraction plus assignment. \precond\\ |v.dimension() == v1.dimension()|.}*/ -Ivector<_RT,_ALLOC> operator+(const Ivector<_RT,_ALLOC>& v1) const; +Ivector<_NT,_ALLOC> operator+(const Ivector<_NT,_ALLOC>& v1) const; /*{\Mbinop Addition. \precond\\ |v.dimension() == v1.dimension()|.}*/ -Ivector<_RT,_ALLOC> operator-(const Ivector<_RT,_ALLOC>& v1) const; +Ivector<_NT,_ALLOC> operator-(const Ivector<_NT,_ALLOC>& v1) const; /*{\Mbinop Subtraction. \precond\\ |v.dimension() = v1.dimension()|.}*/ -RT operator*(const Ivector<_RT,_ALLOC>& v1) const; +NT operator*(const Ivector<_NT,_ALLOC>& v1) const; /*{\Mbinop Inner Product. \precond\\ |v.dimension() = v1.dimension()|.}*/ -Ivector<_RT,_ALLOC> compmul(const RT& r) const; +Ivector<_NT,_ALLOC> compmul(const NT& r) const; -Ivector<_RT,_ALLOC> operator-() const; +Ivector<_NT,_ALLOC> operator-() const; /*{\Munopfunc Negation.}*/ /*{\Mtext We provide component access via |iterator| and |const_iterator| @@ -301,43 +294,43 @@ const_iterator begin() const { return v; } const_iterator end() const { return v+dim; } -bool operator==(const Ivector<_RT,_ALLOC>& w) const; -bool operator!=(const Ivector<_RT,_ALLOC>& w) const +bool operator==(const Ivector<_NT,_ALLOC>& w) const; +bool operator!=(const Ivector<_NT,_ALLOC>& w) const { return !(*this == w); } -static int compare(const Ivector<_RT,_ALLOC>&, - const Ivector<_RT,_ALLOC>&); +static int compare(const Ivector<_NT,_ALLOC>&, + const Ivector<_NT,_ALLOC>&); }; /*{\Mimplementation Vectors are implemented by arrays of type -|RT|. All operations on a vector |v| take time $O(|v.dimension()|)$, +|NT|. All operations on a vector |v| take time $O(|v.dimension()|)$, except for |dimension()| and $[\ ]$ which take constant time. The space requirement is $O(|v.dimension()|)$. }*/ -template -Ivector operator*(const RT& r, const Ivector& v) +template +Ivector operator*(const NT& r, const Ivector& v) /*{\Mbinopfunc Componentwise multiplication with number $r$.}*/ { return v.compmul(r); } -template -Ivector operator*(const Ivector& v, const RT& r) +template +Ivector operator*(const Ivector& v, const NT& r) /*{\Mbinopfunc Componentwise multiplication with number $r$.}*/ { return v.compmul(r); } -template -Ivector operator*(int r, const Ivector& v) +template +Ivector operator*(int r, const Ivector& v) { return v.compmul(r); } -template -Ivector operator*(const Ivector& v, int r) +template +Ivector operator*(const Ivector& v, int r) { return v.compmul(r); } -template -Ivector& Ivector:: -operator+=(const Ivector& vec) +template +Ivector& Ivector:: +operator+=(const Ivector& vec) { check_dimensions(vec); register int n = dim; @@ -345,9 +338,9 @@ operator+=(const Ivector& vec) return *this; } -template -Ivector& Ivector:: -operator-=(const Ivector& vec) +template +Ivector& Ivector:: +operator-=(const Ivector& vec) { check_dimensions(vec); register int n = dim; @@ -355,64 +348,64 @@ operator-=(const Ivector& vec) return *this; } -template -Ivector Ivector:: -operator+(const Ivector& vec) const +template +Ivector Ivector:: +operator+(const Ivector& vec) const { check_dimensions(vec); register int n = dim; - Ivector result(n); + Ivector result(n); while (n--) result.v[n] = v[n]+vec.v[n]; return result; } -template -Ivector Ivector:: -operator-(const Ivector& vec) const +template +Ivector Ivector:: +operator-(const Ivector& vec) const { check_dimensions(vec); register int n = dim; - Ivector result(n); + Ivector result(n); while (n--) result.v[n] = v[n]-vec.v[n]; return result; } -template -Ivector Ivector:: +template +Ivector Ivector:: operator-() const // unary minus { register int n = dim; - Ivector result(n); + Ivector result(n); while (n--) result.v[n] = -v[n]; return result; } -template -Ivector Ivector:: -compmul(const RT& x) const +template +Ivector Ivector:: +compmul(const NT& x) const { int n = dim; - Ivector result(n); + Ivector result(n); while (n--) result.v[n] = v[n] * x; return result; } -template -RT Ivector:: -operator*(const Ivector& vec) const +template +NT Ivector:: +operator*(const Ivector& vec) const { check_dimensions(vec); - RT result=0; + NT result=0; register int n = dim; while (n--) result = result+v[n]*vec.v[n]; return result; } -template -bool Ivector:: -operator==(const Ivector& vec) const +template +bool Ivector:: +operator==(const Ivector& vec) const { if (vec.dim != dim) return false; int i = 0; @@ -420,9 +413,9 @@ operator==(const Ivector& vec) const return (i==dim); } -template -int Ivector:: -compare(const Ivector& v1, const Ivector& v2) +template +int Ivector:: +compare(const Ivector& v1, const Ivector& v2) { register int i; v1.check_dimensions(v2); @@ -431,34 +424,45 @@ compare(const Ivector& v1, const Ivector& v2) return (v1[i] < v2[i]) ? -1 : 1; } -template -std::ostream& operator<<(std::ostream& O, const Ivector<_RT,_ALLOC>& v) +template +std::ostream& operator<<(std::ostream& os, const Ivector<_NT,_ALLOC>& v) /*{\Xbinopfunc writes |\Mvar| componentwise to the output stream $O$.}*/ { /* syntax: d x_0 x_1 ... x_d-1 */ - O << v.dimension() << ' '; - for (register int i = 0; i < v.dimension(); i++) O << v[i] << ' '; - return O; + CGAL::print_d<_NT> prt(&os); + if (os.iword(IO::mode)==IO::PRETTY) os << "LA::Vector("; + prt(v.dimension()); + if (os.iword(IO::mode)==IO::PRETTY) { os << " ["; prt.reset(); } + std::for_each(v.begin(),v.end(),prt); + if (os.iword(IO::mode)==IO::PRETTY) os << "])"; + return os; } -template -std::istream& operator>>(std::istream& in, Ivector<_RT,_ALLOC>& v) +template +std::istream& operator>>(std::istream& is, Ivector<_NT,_ALLOC>& v) /*{\Xbinopfunc reads |\Mvar| componentwise from the input stream $I$.}*/ { /* syntax: d x_0 x_1 ... x_d-1 */ - int dim = 0; - in >> dim; - if (v.dimension() != dim) v = Ivector<_RT,_ALLOC>(dim); - - int i=0; while (i < v.dimension() && in >> v[i++]); - return in; + int dim; + switch (is.iword(IO::mode)) { + case IO::ASCII : + case IO::BINARY : + is >> dim; + v = Ivector<_NT,_ALLOC>(dim); + std::copy_n(std::istream_iterator<_NT>(is),dim,v.begin()); + break; + default: + std::cerr<<"\nStream must be in ascii or binary mode"< -Ivector::allocator_type Ivector::MM; +template +Ivector::allocator_type Ivector::MM; -LA_END_NAMESPACE -#endif // LINALG_IVECTOR_H +CGAL_END_NAMESPACE +#endif // CGAL_IVECTOR_H diff --git a/Packages/Kernel_d/include/CGAL/Kernel_d/Linear_algebraCd.C b/Packages/Kernel_d/include/CGAL/Kernel_d/Linear_algebraCd.C new file mode 100644 index 00000000000..49d7e1fd92a --- /dev/null +++ b/Packages/Kernel_d/include/CGAL/Kernel_d/Linear_algebraCd.C @@ -0,0 +1,468 @@ +// ====================================================================== +// +// Copyright (c) 1999 The CGAL Consortium +// +// This software and related documentation is part of an INTERNAL release +// of the Computational Geometry Algorithms Library (CGAL). It is not +// intended for general use. +// +// ---------------------------------------------------------------------- +// +// release : 1.1 +// release_date : Jan 2001 +// +// file : include/CGAL/Linear_algebraCd.C +// package : Kernel_d +// revision : $Revision$ +// revision_date : $Date$ +// author(s) : Herve.Bronnimann@sophia.inria.fr +// author(s) : seel@mpi-sb.mpg.de +// coordinator : seel@mpi-sb.mpg.de +// +// ====================================================================== + +#ifndef CGAL_LINEAR_ALGEBRACD_C +#define CGAL_LINEAR_ALGEBRACD_C + +#include +#include + +CGAL_BEGIN_NAMESPACE + +template < class FT > +std::pair +Linear_algebraCd:: +transpose(std::pair p) +{ std::swap(p.first,p.second); + return p; +} + +template < class FT > +typename Linear_algebraCd::Matrix +Linear_algebraCd:: +transpose(const Matrix &M) +{ + Matrix P(transpose(M.dimension())); + int i, j; + for (i=0; i +inline // in order to facilitate the optimization with unused variables +void +Linear_algebraCd:: +Gaussian_elimination(const Matrix &M, + // return parameters + Matrix &L, Matrix &U, + std::vector &row_permutation, + std::vector &column_permutation, + FT &det, int &rank, Vector &c) +{ + // Method: we use Gaussian elimination with division at each step + // We do not use the variant by Bareiss (because we are on a field) + // We obtain L, U, q, such that LM' = U, and U is upper diagonal, + // where M' is M whose rows and columns are permuted + // Picked up on the way: + // det, rank, non-trivial solution c if system is degenerate (c^t M = 0) + // Preconditions: + CGAL_kernel_precondition( M.row_dimension() <= M.column_dimension() ); + // Temporaries + int i, j, k; + int dim = M.row_dimension(), cdim = M.column_dimension(); + // All the parameters are already initialized (as in C++) + int sign = 1; + // First create a copy of M into U, and set L and permutations to identity + U = M; + L = Matrix(dim,Matrix::Identity()); + for (i=0; i +bool +Linear_algebraCd:: +Triangular_system_solver(const Matrix &U, const Matrix& L, const Vector &b, + int rank, + // return parameters + Vector &x, FT &D) +{ + // METHOD: solve system Ux=b, returning solution (x/D) + // back substitution of x[rdim], x[rdim-1], etc. + // depends on "free" variables x[rdim+1], etc. x[cdim] + CGAL_kernel_assertion( U.row_dimension() == b.dimension()); + D = FT(1); + for (int i = rank; i < U.column_dimension(); ++i) + if ( b[i] != FT(0) ) + { x = L.row(i); return false; } + + x = Vector(U.column_dimension()); + for (int i = rank-1; i>=0; --i) { + x[i] = b[i]; + for (int j = i+1; j +inline // in order to facilitate the optimization with unused variables +void +Linear_algebraCd:: +Triangular_left_inverse(const Linear_algebraCd::Matrix &U, + // return parameters + Linear_algebraCd::Matrix &Uinv) +{ + int i, j, k; + CGAL_kernel_precondition( U.dimension() == transpose(Uinv.dimension()) ); + //DEBUG: std::cerr << "system : " << U << std::endl; + // std::fill(Uinv.begin(), Uinv.end(), FT(0)); + for (i=U.row_dimension()-1; i>=0; --i) { + Uinv[i][i] = FT(1)/U[i][i]; + for (j=i+1; j +bool +Linear_algebraCd:: +inverse(const Matrix &M, Matrix &I, FT &D, Vector &c) +{ + Matrix L(M.dimension()); + Matrix U(M.dimension()); + Matrix Uinv(M.column_dimension(),M.row_dimension()); + int rank; + std::vector rq, cq; + Gaussian_elimination(M, L, U, rq, cq, D, rank, c); + if (D == FT(0)) return false; // c holds the witness + // Otherwise, compute the inverse of U + Triangular_left_inverse(U,Uinv); + Uinv = Uinv * L; + // Don't forget to permute the rows of M back + + //DEBUG: std::cerr << "inverse before permutation : " << I << std::endl; + for (rank=0; rank +inline +typename Linear_algebraCd::Matrix +Linear_algebraCd:: +inverse(const Matrix &M, FT &D) +{ + CGAL_kernel_precondition( M.row_dimension() == M.column_dimension() ); + Matrix I(M.column_dimension(),M.row_dimension()); + Vector c(M.row_dimension()); + bool result = inverse(M,I,D,c); + CGAL_kernel_precondition( result ); + return I; +} + +template < class FT > +typename Linear_algebraCd::FT +Linear_algebraCd:: +determinant(const Matrix &M, Matrix &L, Matrix &U, + std::vector &q, Vector &c) +{ + FT det; int rank; + std::vector cq; + Gaussian_elimination(M, L, U, q, cq, det, rank, c); + return det; +} + +template < class FT > +inline +typename Linear_algebraCd::FT +Linear_algebraCd:: +determinant(const Matrix &M) +{ + Matrix L(M.dimension()); + Matrix U(M.dimension()); + std::vector q; + Vector c(M.column_dimension()); + return determinant(M,L,U,q,c); +} + +template < class FT > +inline +Sign +Linear_algebraCd:: +sign_of_determinant(const Matrix &M) +{ return CGAL_NTS sign(determinant(M)); } + +template < class FT > +bool +Linear_algebraCd:: +verify_determinant(const Matrix & /*M*/, + const FT & /*D*/, + const Matrix & /*L*/, + const Matrix & /*U*/, + const std::vector & /*q*/, + const Vector & /*c*/) +{ + // TODO: verify_determinant + return true; + CGAL_kernel_assertion( false ); + return false; +} + +template < class FT > +bool +Linear_algebraCd:: +linear_solver(const Matrix &M, const Vector &b, + Vector &x, FT &D, Vector &c) +{ + CGAL_kernel_precondition( M.row_dimension() <= M.column_dimension() ); + CGAL_kernel_precondition( b.dimension() == M.row_dimension() ); + Matrix L(M.dimension()); + Matrix U(M.dimension()); + int rank; + std::vector rq, cq; + + TRACEV(M); TRACEV(b); + Gaussian_elimination(M, L, U, rq, cq, D, rank, c); + TRACEV(U); TRACEV(L); + // Compute a solution by solving triangular system + // Since LM=U, and x is a solution of Mx=b, then Ux=Lb + // Temporary store the solution in c + if ( !Triangular_system_solver(U, L, L*b, rank, c, D) ) + return false; + + // Don't forget to permute the rows of M back + x = Vector(M.column_dimension()); + for (rank=0; rank +bool +Linear_algebraCd:: +linear_solver(const Matrix &M, const Vector &b, + Vector &x, FT &D, Matrix &spanning_vectors, + Vector &c) +{ + CGAL_kernel_precondition( M.row_dimension() <= M.column_dimension() ); + CGAL_kernel_precondition( b.dimension() == M.row_dimension() ); + Matrix L,U; + int rank; + std::vector dummy, var; + + TRACEV(M); TRACEV(b); + Gaussian_elimination(M, L, U, dummy, var, D, rank, c); + TRACEV(U); TRACEV(L); + // Compute a solution by solving triangular system + // Since LM=U, and x is a solution of Mx=b, then Ux=Lb + // Temporary store the solution in c + if ( !Triangular_system_solver(U, L, L*b, rank, c, D) ) + return false; + + // Don't forget to permute the rows of M back: + x = Vector(M.column_dimension()); + for (rank=0; rank 0) { + for(int l=0; l < defect; ++l) { + spanning_vectors(var[rank + l],l)=FT(1); + for(int i = rank - 1; i >= 0 ; i--) { + FT h = - U(i,rank+l); + for (int j= i + 1; j +bool +Linear_algebraCd:: +linear_solver(const Matrix &M, const Vector &b, Vector &x, FT &D) +{ Vector c; + return linear_solver(M, b, x, D, c); +} + +template < class FT > +bool +Linear_algebraCd:: +is_solvable(const Matrix &M, const Vector &b) +{ Vector x; FT D; + return linear_solver(M, b, x, D); +} + +template < class FT > +int +Linear_algebraCd:: +homogeneous_linear_solver(const Matrix &M, Matrix &spanning_vectors) +{ + Matrix L,U; Vector c, b(M.row_dimension()); + std::vector dummy,var; + FT D; int rank; + Gaussian_elimination(M, L, U, dummy, var, D, rank, c); + +#ifdef CGAL_LA_SELFTEST + Vector x; + Triangular_system_solver(U, L, b, rank, c, D); + TRACEV(M);TRACEV(U);TRACEV(b);TRACEV(rank);TRACEV(c);TRACEV(D); + x = Vector(M.column_dimension()); + for (int i=0; i 0) { + /* In the $l$-th spanning vector, $0 \le l < |defect|$ we set + variable |var[rank + l]| to $1 = |denom|/|denom|$ and then the + dependent variables as dictated by the $|rank| + l$ - th column of + |C|.*/ + + for(int l=0; l < defect; ++l) { + spanning_vectors(var[rank + l],l)=FT(1); + for(int i = rank - 1; i >= 0 ; i--) { + FT h = - U(i,rank+l); + for (int j= i + 1; j +bool +Linear_algebraCd:: +homogeneous_linear_solver(const Matrix &M, Vector &x) +{ + x = Vector(M.row_dimension()); + Matrix spanning_vectors; + int defect = homogeneous_linear_solver(M,spanning_vectors); + if (defect == 0) return false; + + /* return first column of |spanning_vectors| */ + for (int i = 0; i < spanning_vectors.row_dimension(); i++) + x[i] = spanning_vectors(i,0); + return true; + +} + + +template < class FT > +int +Linear_algebraCd:: +independent_columns(const Matrix &M, std::vector &q) +{ CGAL_kernel_precondition( M.row_dimension() <= M.column_dimension() ); + int rank; + std::vector dummy; + Matrix Dummy1,Dummy2; Vector Dummy3; FT null(0); + Gaussian_elimination(M, Dummy1, Dummy2, dummy, q, null, rank, Dummy3); + return rank; +} + +template < class FT > +int +Linear_algebraCd:: +rank(const Matrix &M) +{ std::vector q; + if ( M.row_dimension() > M.column_dimension() ) + return independent_columns(transpose(M),q); + return independent_columns(M,q); +} + +CGAL_END_NAMESPACE +#endif // CGAL_LINEAR_ALGEBRACD_C + + diff --git a/Packages/Kernel_d/include/CGAL/Kernel_d/Linear_algebraHd.C b/Packages/Kernel_d/include/CGAL/Kernel_d/Linear_algebraHd.C new file mode 100644 index 00000000000..acde27869b7 --- /dev/null +++ b/Packages/Kernel_d/include/CGAL/Kernel_d/Linear_algebraHd.C @@ -0,0 +1,1036 @@ +// ============================================================================ +// +// Copyright (c) 1997-2000 The CGAL Consortium +// +// This software and related documentation is part of an INTERNAL release +// of the Computational Geometry Algorithms Library (CGAL). It is not +// intended for general use. +// +// ---------------------------------------------------------------------------- +// +// release : $CGAL_Revision$ +// release_date : $CGAL_Date$ +// +// file : include/CGAL/Linear_algebraHd.C +// package : Kernel_d +// chapter : Kernel +// +// source : ddgeo/Linear_algebra.lw +// revision : $Revision$ +// revision_date : $Date$ +// +// author(s) : Michael Seel +// maintainer : Michael Seel +// coordinator : Susan Hert +// +// implementation: Higher dimensional geometry +// ============================================================================ +//--------------------------------------------------------------------- +// file generated by notangle from Linear_algebra.lw +// please debug or modify noweb file +// based on LEDA architecture by S. Naeher, C. Uhrig +// coding: K. Mehlhorn, M. Seel +// debugging and templatization: M. Seel +//--------------------------------------------------------------------- +#ifndef CGAL_LINEAR_ALGEBRAHD_C +#define CGAL_LINEAR_ALGEBRAHD_C +CGAL_BEGIN_NAMESPACE + +template +bool Linear_algebraHd<_RT,_ALLOC>:: +linear_solver(const Matrix& A, const Vector& b, + Vector& x, RT& D, + Matrix& spanning_vectors, Vector& c) +{ + bool solvable = true; + int i,j,k; // indices to step through the matrix + int rows = A.row_dimension(); + int cols = A.column_dimension(); + + /* at this point one might want to check whether the computation can + be carried out with doubles, see section \ref{ optimization }. + */ + + if (b.dimension() != rows) + ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); + + Matrix C(rows,cols + 1); + // the matrix in which we will calculate ($C = (A \vert b)$) + + /* copy |A| and |b| into |C| and L becomes the identity matrix */ + Matrix L(rows); // zero initialized + for(i=0; i var(cols); + // column $j$ of |C| represents the |var[j]| - th variable + // the array is indexed between |0| and |cols - 1| + + for(j=0; j= 0; i--) { + _RT h = C(i,cols) * D; + for (j = i + 1; j < rank; j++) { + h -= C(i,j)*x[var[j]]; + } + x[var[i]]= h / C(i,i); + } + #ifdef CGAL_LA_SELFTEST + CGAL_assertion( (M*x).is_zero() ); + #endif + + int defect = cols - rank; // dimension of kernel + spanning_vectors = Matrix(cols,defect); + + if (defect > 0) { + for(int l=0; l < defect; l++) { + spanning_vectors(var[rank + l],l)=D; + for(i = rank - 1; i >= 0 ; i--) { + _RT h = - C(i,rank + l)*D; + for ( j= i + 1; j +_RT Linear_algebraHd<_RT,_ALLOC>:: +determinant(const Matrix& A) +{ + if (A.row_dimension() != A.column_dimension()) + ERROR_HANDLER(1,"determinant: only square matrices are legal inputs."); + Vector b(A.row_dimension()); // zero - vector + int i,j,k; // indices to step through the matrix + int rows = A.row_dimension(); + int cols = A.column_dimension(); + + /* at this point one might want to check whether the computation can + be carried out with doubles, see section \ref{ optimization }. + */ + + if (b.dimension() != rows) + ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); + + Matrix C(rows,cols + 1); + // the matrix in which we will calculate ($C = (A \vert b)$) + + /* copy |A| and |b| into |C| and L becomes the identity matrix */ + Matrix L(rows); // zero initialized + for(i=0; i var(cols); + // column $j$ of |C| represents the |var[j]| - th variable + // the array is indexed between |0| and |cols - 1| + + for(j=0; j +_RT Linear_algebraHd<_RT,_ALLOC>:: +determinant(const Matrix& A, + Matrix& Ld, Matrix& Ud, + std::vector& q, Ivector<_RT,_ALLOC>& c) +{ + if (A.row_dimension() != A.column_dimension()) + ERROR_HANDLER(1,"determinant: only square matrices are legal inputs."); + Vector b(A.row_dimension()); // zero - vector + int i,j,k; // indices to step through the matrix + int rows = A.row_dimension(); + int cols = A.column_dimension(); + + /* at this point one might want to check whether the computation can + be carried out with doubles, see section \ref{ optimization }. + */ + + if (b.dimension() != rows) + ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); + + Matrix C(rows,cols + 1); + // the matrix in which we will calculate ($C = (A \vert b)$) + + /* copy |A| and |b| into |C| and L becomes the identity matrix */ + Matrix L(rows); // zero initialized + for(i=0; i var(cols); + // column $j$ of |C| represents the |var[j]| - th variable + // the array is indexed between |0| and |cols - 1| + + for(j=0; j +int Linear_algebraHd<_RT,_ALLOC>:: +sign_of_determinant(const Matrix& M) +{ return CGAL_NTS sign(determinant(M)); } + +template +bool Linear_algebraHd<_RT,_ALLOC>:: +verify_determinant(const Matrix& A, _RT D, + Matrix& L, Matrix& U, + const std::vector& q, Ivector<_RT,_ALLOC>& c) +{ + if ((int)q.size() != A.column_dimension()) + ERROR_HANDLER(1,"verify_determinant: \ + q should be a permutation array \ + with index range [0,A.column_dimension() - 1]."); + int n = A.row_dimension(); + int i,j; + if (D == 0) { /* we have $c^T \cdot A = 0$ */ + Ivector<_RT,_ALLOC> zero(n); + return (transpose(A) * c == zero); + } else { + /* we check the conditions on |L| and |U| */ + if (L(0,0) != 1) return false; + for (i = 0; i 0 && L(i,i) != U(i - 1,i - 1)) + return false; + + for (j = i + 1; j < n; j++) + if (L(i,j) != 0) + return false; + } + + /* check whether $L \cdot A \cdot Q = U$ */ + Matrix LA = L * A; + for (j = 0; j < n; j++) + if (LA.column(q[j]) != U.column(j)) + return false; + + /* compute sign |s| of |Q| */ + int sign = 1; + + /* we chase the cycles of |q|. An even length cycle contributes - 1 + and vice versa */ + + std::vector already_considered(n); + + for (i = 0; i < n; i++) + already_considered[i] = false; + + for (i = 0; i < n; i++) + already_considered[q[i]] = true; + + for (i = 0; i < n; i++) + if (! already_considered[i]) + ERROR_HANDLER(1,"verify_determinant:q is not a permutation."); + else + already_considered[i] = false; + + for (i = 0; i < n; i++) { + if (already_considered[i]) continue; + + /* we have found a new cycle with minimal element $i$. */ + int k = q[i]; + already_considered[i] =true; + + while (k != i) { + sign = - sign; + already_considered[k]= true; + k = q[k]; + } + } + return (D == RT(sign) * U(n - 1,n - 1)); + } +} + +template +int Linear_algebraHd<_RT,_ALLOC>:: +independent_columns(const Matrix& A, + std::vector& columns) +{ + Ivector<_RT,_ALLOC> b(A.row_dimension()); // zero - vector + int i,j,k; // indices to step through the matrix + int rows = A.row_dimension(); + int cols = A.column_dimension(); + + /* at this point one might want to check whether the computation can + be carried out with doubles, see section \ref{ optimization }. + */ + + if (b.dimension() != rows) + ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); + + Matrix C(rows,cols + 1); + // the matrix in which we will calculate ($C = (A \vert b)$) + + /* copy |A| and |b| into |C| and L becomes the identity matrix */ + Matrix L(rows); // zero initialized + for(i=0; i var(cols); + // column $j$ of |C| represents the |var[j]| - th variable + // the array is indexed between |0| and |cols - 1| + + for(j=0; j(rank); + for(i = 0; i < rank; i++) + columns[i] = var[i]; + return rank; +} + + +template +int Linear_algebraHd<_RT,_ALLOC>:: +rank(const Matrix& A) +{ + Ivector<_RT,_ALLOC> b(A.row_dimension()); // zero - vector + int i,j,k; // indices to step through the matrix + int rows = A.row_dimension(); + int cols = A.column_dimension(); + + /* at this point one might want to check whether the computation can + be carried out with doubles, see section \ref{ optimization }. + */ + + if (b.dimension() != rows) + ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); + + Matrix C(rows,cols + 1); + // the matrix in which we will calculate ($C = (A \vert b)$) + + /* copy |A| and |b| into |C| and L becomes the identity matrix */ + Matrix L(rows); // zero initialized + for(i=0; i var(cols); + // column $j$ of |C| represents the |var[j]| - th variable + // the array is indexed between |0| and |cols - 1| + + for(j=0; j +bool Linear_algebraHd<_RT,_ALLOC>:: +inverse(const Matrix& A, Matrix& inverse, + _RT& D, Ivector<_RT,_ALLOC>& c) +{ + if (A.row_dimension() != A.column_dimension()) + ERROR_HANDLER(1,"inverse: only square matrices are legal inputs."); + Ivector<_RT,_ALLOC> b(A.row_dimension()); // zero - vector + int i,j,k; // indices to step through the matrix + int rows = A.row_dimension(); + int cols = A.column_dimension(); + + /* at this point one might want to check whether the computation can + be carried out with doubles, see section \ref{ optimization }. + */ + + if (b.dimension() != rows) + ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); + + Matrix C(rows,cols + 1); + // the matrix in which we will calculate ($C = (A \vert b)$) + + /* copy |A| and |b| into |C| and L becomes the identity matrix */ + Matrix L(rows); // zero initialized + for(i=0; i var(cols); + // column $j$ of |C| represents the |var[j]| - th variable + // the array is indexed between |0| and |cols - 1| + + for(j=0; j= 0; j--) { + h = L (j,i) * D; + for (int l = j + 1; l +int Linear_algebraHd<_RT,_ALLOC>:: +homogeneous_linear_solver(const Matrix &A, + Matrix& spanning_vectors) +/* returns the dimension of the solution space of the homogeneous system + $Ax = 0$. The columns of spanning\_vectors span the solution space. */ +{ + Ivector<_RT,_ALLOC> b(A.row_dimension()); // zero - vector + _RT D; + int i,j,k; // indices to step through the matrix + int rows = A.row_dimension(); + int cols = A.column_dimension(); + + /* at this point one might want to check whether the computation can + be carried out with doubles, see section \ref{ optimization }. + */ + + if (b.dimension() != rows) + ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); + + Matrix C(rows,cols + 1); + // the matrix in which we will calculate ($C = (A \vert b)$) + + /* copy |A| and |b| into |C| and L becomes the identity matrix */ + Matrix L(rows); // zero initialized + for(i=0; i var(cols); + // column $j$ of |C| represents the |var[j]| - th variable + // the array is indexed between |0| and |cols - 1| + + for(j=0; j x; + x = Vector(cols); + D = denom; + for(i = rank - 1; i >= 0; i--) { + _RT h = C(i,cols) * D; + for (j = i + 1; j < rank; j++) { + h -= C(i,j)*x[var[j]]; + } + x[var[i]]= h / C(i,i); + } + #ifdef CGAL_LA_SELFTEST + CGAL_assertion( (M*x).is_zero() ); + #endif + + int defect = cols - rank; // dimension of kernel + spanning_vectors = Matrix(cols,defect); + + if (defect > 0) { + for(int l=0; l < defect; l++) { + spanning_vectors(var[rank + l],l)=D; + for(i = rank - 1; i >= 0 ; i--) { + _RT h = - C(i,rank + l)*D; + for ( j= i + 1; j +bool Linear_algebraHd<_RT,_ALLOC>:: +homogeneous_linear_solver(const Matrix& A, Vector& x) +/* returns true if the homogeneous system $Ax = 0$ has a non - trivial + solution and false otherwise. */ +{ + x = Vector(A.row_dimension()); + Matrix spanning_vectors; + int dimension = homogeneous_linear_solver(A,spanning_vectors); + if (dimension == 0) return false; + + /* return first column of |spanning_vectors| */ + for (int i = 0; i < spanning_vectors.row_dimension(); i++) + x[i] = spanning_vectors(i,0); + return true; +} + + +template +typename Linear_algebraHd<_RT,_ALLOC>::Matrix +Linear_algebraHd<_RT,_ALLOC>::transpose(const Matrix& M) +{ + int d1 = M.row_dimension(); + int d2 = M.column_dimension(); + Matrix result(d2,d1); + for(int i = 0; i < d2; i++) + for(int j = 0; j < d1; j++) + result(i,j) = M(j,i); + return result; +} + + +CGAL_END_NAMESPACE +#endif //CGAL_LINEAR_ALGEBRAHD_C + + diff --git a/Packages/Kernel_d/include/CGAL/Kernel_d/Vector_d.h b/Packages/Kernel_d/include/CGAL/Kernel_d/Vector_d.h new file mode 100644 index 00000000000..74a1c62c3fc --- /dev/null +++ b/Packages/Kernel_d/include/CGAL/Kernel_d/Vector_d.h @@ -0,0 +1,77 @@ +#ifndef CGAL_VECTOR_D_H +#define CGAL_VECTOR_D_H + +CGAL_BEGIN_NAMESPACE + +template +class Vector_d : public pR::Vector_d_base +{ public: + typedef typename pR::Vector_d_base Base; + typedef Vector_d Self; + typedef pR R; + typedef typename R::RT RT; + typedef typename R::FT FT; + typedef typename R::LA LA; + + Vector_d(int d=0) : Base(d) {} + Vector_d(int d, Null_vector v) : Base(d,v) {} + Vector_d(int a, int b, int c = 1) : Base(a,b,c) {} + Vector_d(const RT& a, const RT& b, const RT& c = 1) : + Base(a,b,c) {} + Vector_d(int a, int b, int c, int d) : Base(a,b,c,d) {} + Vector_d(const RT& a, const RT& b, const RT& c, const RT& d) : + Base(a,b,c,d) {} + Vector_d(typename Base::Base_vector, int d, int i) : + Base(typename Base::Base_vector(), d,i) {} + + template + Vector_d (int d, InputIterator first, InputIterator last) + : Base (d, first, last) {} + template + Vector_d (int d, InputIterator first, InputIterator last, const RT& D) + : Base (d, first, last, D) {} + Vector_d(const Self& v) : Base(v) {} + Vector_d(const Base& v) : Base(v) {} + + Direction_d direction() const { return Base::direction(); } + + FT operator* (const Self& w) const + { return Base::operator*(w); } + Self operator+(const Self& w) const + { return Base::operator+(w); } + Self operator-(const Self& w) const + { return Base::operator-(w); } + Self operator-() const + { return Base::operator-(); } + + template + Self operator/(const NT& n) const { return Base::operator/(n); } + + Self& operator+=(const Self& w) + { return static_cast(Base::operator+=(w)); } + Self& operator-=(const Self& w) + { return static_cast(Base::operator-=(w)); } + template + Self& operator*=(const NT& n) + { return static_cast(Base::operator*=(n)); } + template + Self& operator/=(const NT& n) + { return static_cast(Base::operator/=(n)); } + + bool operator==(const Self& w) const + { return Base::operator==(w); } + bool operator!=(const Self& w) const + { return Base::operator!=(w); } + +}; + +template Point_d +operator+ (const Origin& o, const Vector_d& v) +{ return Point_d( o + static_cast::Base>(v) ); } + +template +Vector_d operator*(const NT& n, const Vector_d& v) +{ return Vector_d( n * static_cast::Base>(v) ); } + +CGAL_END_NAMESPACE +#endif //CGAL_VECTOR_D_H diff --git a/Packages/Kernel_d/include/CGAL/Kernel_d/d_utils.h b/Packages/Kernel_d/include/CGAL/Kernel_d/d_utils.h new file mode 100644 index 00000000000..75fe7ae3313 --- /dev/null +++ b/Packages/Kernel_d/include/CGAL/Kernel_d/d_utils.h @@ -0,0 +1,67 @@ +// ====================================================================== +// +// Copyright (c) 1999 The CGAL Consortium +// +// This software and related documentation is part of an INTERNAL release +// of the Computational Geometry Algorithms Library (CGAL). It is not +// intended for general use. +// +// ---------------------------------------------------------------------- +// +// release : 1.1 +// release_date : %-e Dec 1999 +// +// file : include/CGAL/Cartesian/d_utils.h +// package : Cd (1.1) +// revision : $Revision$ +// revision_date : $Date$ +// author : Herve Brönnimann +// coordinator : INRIA Sophia-Antipolis (Herve.Bronnimann@sophia.inria.fr) +// +// ====================================================================== + +#ifndef CGAL_CARTESIAN_D_UTILS_H +#define CGAL_CARTESIAN_D_UTILS_H + +#include +#include +#include + +CGAL_BEGIN_NAMESPACE + +// Small object utility for printing objects +// with separators depending on the ostream mode +template < class Object > +struct print_d +{ + char * _separator; + std::ostream*_os; + bool _print_sep; + + print_d(std::ostream *os) : _os(os), _print_sep(false) + { + if (_os->iword(IO::mode)==IO::ASCII) _separator = " "; + else if (_os->iword(IO::mode)==IO::BINARY) _separator = ""; + else _separator = ", "; + } + void reset() + { + _print_sep = false; + } + void operator()(const Object &x) { + if (_print_sep && _os->iword(IO::mode) != IO::BINARY) + (*_os) << _separator; + (*_os) << x; + _print_sep = true; + } + void operator()(const int &x) { + if (_print_sep && _os->iword(IO::mode) != IO::BINARY) + (*_os) << _separator; + (*_os) << x; + _print_sep = true; + } +}; + +CGAL_END_NAMESPACE + +#endif // CGAL_CARTESIAN_D_UTILS_H diff --git a/Packages/Kernel_d/include/CGAL/Linear_algebra.h b/Packages/Kernel_d/include/CGAL/Linear_algebra.h deleted file mode 100644 index 92b4d225eb9..00000000000 --- a/Packages/Kernel_d/include/CGAL/Linear_algebra.h +++ /dev/null @@ -1,1323 +0,0 @@ -// ============================================================================ -// -// Copyright (c) 1997-2000 The CGAL Consortium -// -// This software and related documentation is part of an INTERNAL release -// of the Computational Geometry Algorithms Library (CGAL). It is not -// intended for general use. -// -// ---------------------------------------------------------------------------- -// -// release : $CGAL_Revision$ -// release_date : $CGAL_Date$ -// -// file : include/CGAL/Linear_algebra.h -// package : Kernel_d -// chapter : Kernel -// -// source : ddgeo/Linear_algebra.lw -// revision : $Revision$ -// revision_date : $Date$ -// -// author(s) : Michael Seel -// maintainer : Michael Seel -// coordinator : Susan Hert -// -// implementation: Higher dimensional geometry -// ============================================================================ -//--------------------------------------------------------------------- -// file generated by notangle from Linear_algebra.lw -// please debug or modify noweb file -// based on LEDA architecture by S. Naeher, C. Uhrig -// coding: K. Mehlhorn, M. Seel -// debugging and templatization: M. Seel -//--------------------------------------------------------------------- - - -#ifndef LINEAR_ALGEBRA_H -#define LINEAR_ALGEBRA_H - -#include -#include - -// #define CGAL_LA_SELFTEST -LA_BEGIN_NAMESPACE - -/*{\Moptions outfile=Linear_algebra.man}*/ -/*{\Manpage {Linear_algebra} {RT} {Linear Algebra on RT} {LA}}*/ - -template -class Linear_algebra -{ -/*{\Mdefinition -The data type |\Mname| encapsulates two classes |Imatrix|, |Ivector| -and many functions of basic linear algebra. It is parametrized by a -number type |RT|. An instance of data type |Imatrix| is a matrix of -variables of type |RT|, the so called ring type. Accordingly, -|Ivector| implements vectors of variables of type |RT|. The arithmetic -type |RT| is required to behave like integers in the mathematical -sense. - -All functions compute the exact result, i.e., there is no rounding -error. Most functions of linear algebra are \emph{checkable}, i.e., -the programs can be asked for a proof that their output is -correct. For example, if the linear system solver declares a linear -system $A x = b$ unsolvable it also returns a vector $c$ such that -$c^T A = 0$ and $c^T b \neq 0$. All internal correctness checks can -be switched on by the flag [[CGAL_LA_SELFTEST]].}*/ - -public: - -/*{\Mtypes 5.5}*/ - -typedef _RT RT; -/*{\Mtypemember the ring type of the components.}*/ - -typedef Ivector<_RT,_ALLOC> Vector; -/*{\Mtypemember the vector type.}*/ - -typedef Imatrix<_RT,_ALLOC> Matrix; -/*{\Mtypemember the matrix type.}*/ - -typedef _ALLOC allocator_type; -/*{\Mtypemember the allocator used for memory management. |\Mname| is -an abbreviation for |Linear_algebra >|. Thus -|allocator_type| defaults to the standard allocator offered by the STL.}*/ - -/*{\Moperations 2 1}*/ - -static Matrix transpose(const Matrix& M); -/*{\Mstatic returns $M^T$ ($m\times n$ - matrix). }*/ - -static bool inverse(const Matrix& M, Matrix& I, RT& D, Vector& c); -/*{\Mstatic determines whether |M| has an inverse. It also computes -either the inverse as $(1/D) \cdot |I|$ or when no inverse -exists, a vector $c$ such that $c^T \cdot M = 0 $. }*/ - -static Matrix inverse(const Matrix& M, RT& D) -/*{\Mstatic returns the inverse matrix of |M|. More precisely, $1/D$ - times the matrix returned is the inverse of |M|.\\ - \precond |determinant(M) != 0|. }*/ -{ - Matrix result; - Vector c; - if (!inverse(M,result,D,c)) - ERROR_HANDLER(1,"Imatrix::inverse: matrix is singular."); - return result; -} - -static RT determinant (const Matrix& M, Matrix& L, Matrix& U, - std::vector& q, Vector& c); -/*{\Mstatic returns the determinant $D$ of |M| and sufficient information - to verify that the value of the determinant is correct. If - the determinant is zero then $c$ is a vector such that - $c^T \cdot M = 0$. If the determinant is non-zero then $L$ - and $U$ are lower and upper diagonal matrices respectively - and $q$ encodes a permutation matrix $Q$ with $Q(i,j) = 1$ - iff $i = q(j)$ such that $L \cdot M \cdot Q = U$, - $L(0,0) = 1$, $L(i,i) = U(i - 1,i - 1)$ for all $i$, - $1 \le i < n$, and $D = s \cdot U(n - 1,n - 1)$ where $s$ is - the determinant of $Q$. \precond |M| is square. }*/ - -static bool verify_determinant (const Matrix& M, RT D, Matrix& L, Matrix& U, - const std::vector& q, Vector& c); -/*{\Mstatic verifies the conditions stated above. }*/ - -static RT determinant (const Matrix& M); -/*{\Mstatic returns the determinant of |M|. - \precond |M| is square. }*/ - -static int sign_of_determinant (const Matrix& M); -/*{\Mstatic returns the sign of the determinant of |M|. - \precond |M| is square. }*/ - -static bool linear_solver(const Matrix& M, const Vector& b, - Vector& x, RT& D, - Matrix& spanning_vectors, - Vector& c); -/*{\Mstatic determines the complete solution space of the linear system - $M\cdot x = b$. If the system is unsolvable then - $c^T \cdot M = 0$ and $c^T \cdot b \not= 0$. - If the system is solvable then $(1/D) x$ is a solution, and - the columns of |spanning_vectors| are a maximal set of linearly - independent solutions to the corresponding homogeneous system. - \precond |M.row_dimension() = b.dimension()|. }*/ - -static bool linear_solver(const Matrix& M, const Vector& b, - Vector& x, RT& D, - Vector& c) -/*{\Mstatic determines whether the linear system $M\cdot x = b$ is - solvable. If yes, then $(1/D) x$ is a solution, if not then - $c^T \cdot M = 0$ and $c^T \cdot b \not= 0$. - \precond |M.row_dimension() = b.dimension()|. }*/ -{ - Matrix spanning_vectors; - return linear_solver(M,b,x,D,spanning_vectors,c); -} - -static bool linear_solver(const Matrix& M, const Vector& b, - Vector& x, RT& D) -/*{\Mstatic as above, but without the witness $c$ - \precond |M.row_dimension() = b.dimension()|. }*/ -{ - Matrix spanning_vectors; Vector c; - return linear_solver(M,b,x,D,spanning_vectors,c); -} - -static bool is_solvable(const Matrix& M, const Vector& b) -/*{\Mstatic determines whether the system $M \cdot x = b$ is solvable \\ - \precond |M.row_dimension() = b.dimension()|. }*/ -{ - Vector x; RT D; Matrix spanning_vectors; Vector c; - return linear_solver(M,b,x,D,spanning_vectors,c); -} - -static bool homogeneous_linear_solver (const Matrix& M, Vector& x); -/*{\Mstatic determines whether the homogeneous linear system - $M\cdot x = 0$ has a non - trivial solution. If - yes, then $x$ is such a solution. }*/ - -static int homogeneous_linear_solver (const Matrix& M, Matrix& spanning_vecs); -/*{\Mstatic determines the solution space of the homogeneous linear -system $M\cdot x = 0$. It returns the dimension of the solution space. -Moreover the columns of |spanning_vecs| span the solution space. }*/ - -static void independent_columns (const Matrix& M, std::vector& columns); -/*{\Mstatic returns the indices of a maximal subset of independent -columns of |M|.}*/ - -static int rank (const Matrix & M); -/*{\Mstatic returns the rank of matrix |M| }*/ - -/*{\Mimplementation The datatype |\Mname| is a wrapper class for the -linear algebra functionality on matrices and vectors. Operations -|determinant|, |inverse|, |linear_solver|, and |rank| take time -$O(n^3)$, and all other operations take time $O(nm)$. These time -bounds ignore the cost for multiprecision arithmetic operations. - -All functions on integer matrices compute the exact result, i.e., -there is no rounding error. The implemenation follows a proposal of -J. Edmonds (J. Edmonds, Systems of distinct representatives and linear -algebra, Journal of Research of the Bureau of National Standards, (B), -71, 241 - 245). Most functions of linear algebra are { \em checkable -}, i.e., the programs can be asked for a proof that their output is -correct. For example, if the linear system solver declares a linear -system $A x = b$ unsolvable it also returns a vector $c$ such that -$c^T A = 0$ and $c^T b \not= 0$.}*/ - -}; - - -template -bool Linear_algebra<_RT,_ALLOC>:: -linear_solver(const Imatrix<_RT,_ALLOC>& A, const Vector& b, - Vector& x, RT& D, - Imatrix<_RT,_ALLOC>& spanning_vectors, Vector& c) -{ - bool solvable = true; - - int i,j,k; // indices to step through the matrix - int rows = A.row_dimension(); - int cols = A.column_dimension(); - - /* at this point one might want to check whether the computation can - be carried out with doubles, see section \ref{ optimization }. - */ - - if (b.dimension() != rows) - ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); - - Imatrix<_RT,_ALLOC> C(rows,cols + 1); - // the matrix in which we will calculate ($C = (A \vert b)$) - - /* copy |A| and |b| into |C| and L becomes the identity matrix */ - Imatrix<_RT,_ALLOC> L(rows); // zero initialized - for(i=0; i var(cols); - // column $j$ of |C| represents the |var[j]| - th variable - // the array is indexed between |0| and |cols - 1| - - for(j=0; j= 0; i--) { - _RT h =C(i,cols) * D; - for (j = i + 1; j 0) { - - /* In the $l$ - th spanning vector, $0 \le l < |dimension|$ we set - variable |var[rank + l]| to $1 = |denom|/|denom|$ and then the - dependent variables as dictated by the $|rank| + l$ - th column of - |C|.*/ - - for(int l=0; l < dimension; l++) { - spanning_vectors(var[rank + l],l)=D; - for(i = rank - 1; i >= 0 ; i--) { - _RT h = - C(i,rank + l)* D; - for ( j= i + 1; j -_RT Linear_algebra<_RT,_ALLOC>:: -determinant(const Imatrix<_RT,_ALLOC>& A) -{ - if (A.row_dimension() != A.column_dimension()) - ERROR_HANDLER(1,"determinant: only square matrices are legal inputs."); - Vector b(A.row_dimension()); // zero - vector - - int i,j,k; // indices to step through the matrix - int rows = A.row_dimension(); - int cols = A.column_dimension(); - - /* at this point one might want to check whether the computation can - be carried out with doubles, see section \ref{ optimization }. - */ - - if (b.dimension() != rows) - ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); - - Imatrix<_RT,_ALLOC> C(rows,cols + 1); - // the matrix in which we will calculate ($C = (A \vert b)$) - - /* copy |A| and |b| into |C| and L becomes the identity matrix */ - Imatrix<_RT,_ALLOC> L(rows); // zero initialized - for(i=0; i var(cols); - // column $j$ of |C| represents the |var[j]| - th variable - // the array is indexed between |0| and |cols - 1| - - for(j=0; j -_RT Linear_algebra<_RT,_ALLOC>:: -determinant(const Imatrix<_RT,_ALLOC>& A, - Imatrix<_RT,_ALLOC>& Ld, Imatrix<_RT,_ALLOC>& Ud, - std::vector& q, Ivector<_RT,_ALLOC>& c) -{ - if (A.row_dimension() != A.column_dimension()) - ERROR_HANDLER(1,"determinant: only square matrices are legal inputs."); - Vector b(A.row_dimension()); // zero - vector - - int i,j,k; // indices to step through the matrix - int rows = A.row_dimension(); - int cols = A.column_dimension(); - - /* at this point one might want to check whether the computation can - be carried out with doubles, see section \ref{ optimization }. - */ - - if (b.dimension() != rows) - ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); - - Imatrix<_RT,_ALLOC> C(rows,cols + 1); - // the matrix in which we will calculate ($C = (A \vert b)$) - - /* copy |A| and |b| into |C| and L becomes the identity matrix */ - Imatrix<_RT,_ALLOC> L(rows); // zero initialized - for(i=0; i var(cols); - // column $j$ of |C| represents the |var[j]| - th variable - // the array is indexed between |0| and |cols - 1| - - for(j=0; j -int Linear_algebra<_RT,_ALLOC>:: -sign_of_determinant(const Imatrix<_RT,_ALLOC>& M) -{ return CGAL_NTS sign(determinant(M)); } - -template -bool Linear_algebra<_RT,_ALLOC>:: -verify_determinant(const Imatrix<_RT,_ALLOC>& A, _RT D, - Imatrix<_RT,_ALLOC>& L, Imatrix<_RT,_ALLOC>& U, - const std::vector& q, Ivector<_RT,_ALLOC>& c) -{ - if ((int)q.size() != A.column_dimension()) - ERROR_HANDLER(1,"verify_determinant: \ - q should be a permutation array \ - with index range [0,A.column_dimension() - 1]."); - int n = A.row_dimension(); - int i,j; - if (D == 0) { /* we have $c^T \cdot A = 0$ */ - Ivector<_RT,_ALLOC> zero(n); - return (transpose(A) * c == zero); - } else { - /* we check the conditions on |L| and |U| */ - if (L(0,0) != 1) return false; - for (i = 0; i 0 && L(i,i) != U(i - 1,i - 1)) - return false; - - for (j = i + 1; j < n; j++) - if (L(i,j) != 0) - return false; - } - - /* check whether $L \cdot A \cdot Q = U$ */ - Imatrix<_RT,_ALLOC> LA = L * A; - for (j = 0; j < n; j++) - if (LA.column(q[j]) != U.column(j)) - return false; - - /* compute sign |s| of |Q| */ - int sign = 1; - - /* we chase the cycles of |q|. An even length cycle contributes - 1 - and vice versa */ - - std::vector already_considered(n); - - for (i = 0; i < n; i++) - already_considered[i] = false; - - for (i = 0; i < n; i++) - already_considered[q[i]] = true; - - for (i = 0; i < n; i++) - if (! already_considered[i]) - ERROR_HANDLER(1,"verify_determinant:q is not a permutation."); - else - already_considered[i] = false; - - for (i = 0; i < n; i++) { - if (already_considered[i]) continue; - - /* we have found a new cycle with minimal element $i$. */ - int k = q[i]; - already_considered[i] =true; - - while (k != i) { - sign = - sign; - already_considered[k]= true; - k = q[k]; - } - } - return (D == RT(sign) * U(n - 1,n - 1)); - } -} - -template -void Linear_algebra<_RT,_ALLOC>:: -independent_columns(const Imatrix<_RT,_ALLOC>& A, - std::vector& columns) -{ - Ivector<_RT,_ALLOC> b(A.row_dimension()); // zero - vector - - int i,j,k; // indices to step through the matrix - int rows = A.row_dimension(); - int cols = A.column_dimension(); - - /* at this point one might want to check whether the computation can - be carried out with doubles, see section \ref{ optimization }. - */ - - if (b.dimension() != rows) - ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); - - Imatrix<_RT,_ALLOC> C(rows,cols + 1); - // the matrix in which we will calculate ($C = (A \vert b)$) - - /* copy |A| and |b| into |C| and L becomes the identity matrix */ - Imatrix<_RT,_ALLOC> L(rows); // zero initialized - for(i=0; i var(cols); - // column $j$ of |C| represents the |var[j]| - th variable - // the array is indexed between |0| and |cols - 1| - - for(j=0; j(rank); - for(i = 0; i < rank; i++) - columns[i] = var[i]; -} - - -template -int Linear_algebra<_RT,_ALLOC>:: -rank(const Imatrix<_RT,_ALLOC>& A) -{ - Ivector<_RT,_ALLOC> b(A.row_dimension()); // zero - vector - - int i,j,k; // indices to step through the matrix - int rows = A.row_dimension(); - int cols = A.column_dimension(); - - /* at this point one might want to check whether the computation can - be carried out with doubles, see section \ref{ optimization }. - */ - - if (b.dimension() != rows) - ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); - - Imatrix<_RT,_ALLOC> C(rows,cols + 1); - // the matrix in which we will calculate ($C = (A \vert b)$) - - /* copy |A| and |b| into |C| and L becomes the identity matrix */ - Imatrix<_RT,_ALLOC> L(rows); // zero initialized - for(i=0; i var(cols); - // column $j$ of |C| represents the |var[j]| - th variable - // the array is indexed between |0| and |cols - 1| - - for(j=0; j -bool Linear_algebra<_RT,_ALLOC>:: -inverse(const Imatrix<_RT,_ALLOC>& A, Imatrix<_RT,_ALLOC>& inverse, - _RT& D, Ivector<_RT,_ALLOC>& c) -{ - if (A.row_dimension() != A.column_dimension()) - ERROR_HANDLER(1,"inverse: only square matrices are legal inputs."); - Ivector<_RT,_ALLOC> b(A.row_dimension()); // zero - vector - - int i,j,k; // indices to step through the matrix - int rows = A.row_dimension(); - int cols = A.column_dimension(); - - /* at this point one might want to check whether the computation can - be carried out with doubles, see section \ref{ optimization }. - */ - - if (b.dimension() != rows) - ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); - - Imatrix<_RT,_ALLOC> C(rows,cols + 1); - // the matrix in which we will calculate ($C = (A \vert b)$) - - /* copy |A| and |b| into |C| and L becomes the identity matrix */ - Imatrix<_RT,_ALLOC> L(rows); // zero initialized - for(i=0; i var(cols); - // column $j$ of |C| represents the |var[j]| - th variable - // the array is indexed between |0| and |cols - 1| - - for(j=0; j(rows); //square - _RT h; - for(i = 0; i = 0; j--) - { - h = L (j,i) * D; - for (int l = j + 1; l -int Linear_algebra<_RT,_ALLOC>:: -homogeneous_linear_solver(const Imatrix<_RT,_ALLOC> &A, - Imatrix<_RT,_ALLOC>& spanning_vectors) -/* returns the dimension of the solution space of the homogeneous system - $Ax = 0$. The columns of spanning\_vectors span the solution space. */ -{ - Ivector<_RT,_ALLOC> b(A.row_dimension()); // zero - vector - _RT D; - - int i,j,k; // indices to step through the matrix - int rows = A.row_dimension(); - int cols = A.column_dimension(); - - /* at this point one might want to check whether the computation can - be carried out with doubles, see section \ref{ optimization }. - */ - - if (b.dimension() != rows) - ERROR_HANDLER(1,"linear_solver: b has wrong dimension"); - - Imatrix<_RT,_ALLOC> C(rows,cols + 1); - // the matrix in which we will calculate ($C = (A \vert b)$) - - /* copy |A| and |b| into |C| and L becomes the identity matrix */ - Imatrix<_RT,_ALLOC> L(rows); // zero initialized - for(i=0; i var(cols); - // column $j$ of |C| represents the |var[j]| - th variable - // the array is indexed between |0| and |cols - 1| - - for(j=0; j x; - x = Vector(cols); - D = denom; - for(i=rank - 1 ; i>= 0; i--) { - _RT h =C(i,cols) * D; - for (j = i + 1; j 0) { - - /* In the $l$ - th spanning vector, $0 \le l < |dimension|$ we set - variable |var[rank + l]| to $1 = |denom|/|denom|$ and then the - dependent variables as dictated by the $|rank| + l$ - th column of - |C|.*/ - - for(int l=0; l < dimension; l++) { - spanning_vectors(var[rank + l],l)=D; - for(i = rank - 1; i >= 0 ; i--) { - _RT h = - C(i,rank + l)* D; - for ( j= i + 1; j -bool Linear_algebra<_RT,_ALLOC>:: -homogeneous_linear_solver(const Imatrix<_RT,_ALLOC>& A, Ivector<_RT,_ALLOC>& x) -/* returns true if the homogeneous system $Ax = 0$ has a non - trivial - solution and false otherwise. */ -{ - Imatrix<_RT,_ALLOC> spanning_vectors; - int dimension = homogeneous_linear_solver(A,spanning_vectors); - - if (dimension == 0) - return false; - - /* return first column of |spanning_vectors| */ - for (int i = 0; i < spanning_vectors.row_dimension() ; i++) - x[i] = spanning_vectors(i,0); - return true; -} - - - -template -Imatrix<_RT,_ALLOC> Linear_algebra<_RT,_ALLOC>:: -transpose(const Imatrix<_RT,_ALLOC>& M) -{ - int d1 = M.row_dimension(); - int d2 = M.column_dimension(); - Imatrix<_RT,_ALLOC> result(d2,d1); - for(int i = 0; i < d2; i++) - for(int j = 0; j < d1; j++) - result(i,j) = M(j,i); - return result; -} - - -LA_END_NAMESPACE -#endif // LINALG_ALGEBRA_H - diff --git a/Packages/Kernel_d/include/CGAL/Linear_algebraCd.h b/Packages/Kernel_d/include/CGAL/Linear_algebraCd.h new file mode 100644 index 00000000000..1a741ef58a3 --- /dev/null +++ b/Packages/Kernel_d/include/CGAL/Linear_algebraCd.h @@ -0,0 +1,122 @@ +// ====================================================================== +// +// Copyright (c) 1999 The CGAL Consortium +// +// This software and related documentation is part of an INTERNAL release +// of the Computational Geometry Algorithms Library (CGAL). It is not +// intended for general use. +// +// ---------------------------------------------------------------------- +// +// release : 1.1 +// release_date : Jan 2001 +// +// file : include/CGAL/Linear_algebra_c.h +// package : Kernel_d +// revision : $Revision$ +// revision_date : $Date$ +// author(s) : Herve.Bronnimann@sophia.inria.fr +// author(s) : seel@mpi-sb.mpg.de +// coordinator : seel@mpi-sb.mpg.de +// +// ====================================================================== + +#ifndef CGAL_LINEAR_ALGEBRACD_H +#define CGAL_LINEAR_ALGEBRACD_H + +#include +#include +#include +#include +#undef _DEBUG +#define _DEBUG 13 +#include + +CGAL_BEGIN_NAMESPACE + +template < class _FT > +class Linear_algebraCd +{ +public: + typedef _FT FT; + typedef _FT RT; + typedef Linear_algebraCd Self; + typedef Ivector Vector; + typedef Imatrix Matrix; + typedef const FT* const_iterator; + typedef FT* iterator; + + Linear_algebraCd() {} + +protected: + // Major routines for Linear_algebra_d + static + void Gaussian_elimination(const Matrix &M, + Matrix &L, Matrix &U, + std::vector &row_permutation, + std::vector &column_permutation, + FT &det, int &rank, Vector &c); + static + bool Triangular_system_solver(const Matrix &U, const Matrix &L, + const Vector &b, int rank, Vector &x, FT &det); + static + void Triangular_left_inverse(const Matrix &U, + Matrix &Uinv); + +public: + static + std::pair transpose(std::pair dim); + static + Matrix transpose(const Matrix &M); + + static + bool inverse(const Matrix &M, Matrix &I, FT &D, Vector &c); + + static + Matrix inverse(const Matrix &M, RT &D); + + static + FT determinant(const Matrix &M, Matrix &L, Matrix &U, + std::vector &q, Vector &c); + static + FT determinant(const Matrix &M); + + static + Sign sign_of_determinant(const Matrix &M); + + static + bool verify_determinant(const Matrix &M, const RT &D, + const Matrix &L, const Matrix &U, + const std::vector &q, const Vector &c); + static + bool linear_solver(const Matrix &M, const Vector &b, + Vector &x, RT &D, Matrix &spanning_vectors, Vector &c); + static + bool linear_solver(const Matrix &M, const Vector &b, + Vector &x, RT &D, Vector &c); + static + bool linear_solver(const Matrix &M, const Vector &b, + Vector &x, RT &D); + static + bool is_solvable(const Matrix &M, const Vector &b); + + static + bool homogeneous_linear_solver(const Matrix &M, Vector &x); + static + int homogeneous_linear_solver(const Matrix &M, + Matrix &spanning_vectors); + static + int rank(const Matrix &M); + + static + int independent_columns(const Matrix& M, std::vector& columns); + +}; + +CGAL_END_NAMESPACE + +#include + +#endif // CGAL_LINEAR_ALGEBRACD_H + + diff --git a/Packages/Kernel_d/include/CGAL/Linear_algebraHd.h b/Packages/Kernel_d/include/CGAL/Linear_algebraHd.h new file mode 100644 index 00000000000..d8dfe41b81c --- /dev/null +++ b/Packages/Kernel_d/include/CGAL/Linear_algebraHd.h @@ -0,0 +1,213 @@ +// ============================================================================ +// +// Copyright (c) 1997-2000 The CGAL Consortium +// +// This software and related documentation is part of an INTERNAL release +// of the Computational Geometry Algorithms Library (CGAL). It is not +// intended for general use. +// +// ---------------------------------------------------------------------------- +// +// release : $CGAL_Revision$ +// release_date : $CGAL_Date$ +// +// file : include/CGAL/Linear_algebraHd.h +// package : Kernel_d +// chapter : Kernel +// +// source : ddgeo/Linear_algebra.lw +// revision : $Revision$ +// revision_date : $Date$ +// +// author(s) : Michael Seel +// maintainer : Michael Seel +// coordinator : Susan Hert +// +// implementation: Higher dimensional geometry +// ============================================================================ +//--------------------------------------------------------------------- +// file generated by notangle from Linear_algebra.lw +// please debug or modify noweb file +// based on LEDA architecture by S. Naeher, C. Uhrig +// coding: K. Mehlhorn, M. Seel +// debugging and templatization: M. Seel +//--------------------------------------------------------------------- +#ifndef CGAL_LINEAR_ALGEBRAHD_H +#define CGAL_LINEAR_ALGEBRAHD_H + +#include +#include + +// #define CGAL_LA_SELFTEST +CGAL_BEGIN_NAMESPACE + +/*{\Moptions outfile=Linear_algebra.man}*/ +/*{\Manpage {Linear_algebraHd}{RT}{Linear Algebra on RT}{LA}}*/ + +template +class Linear_algebraHd +{ +/*{\Mdefinition +The data type |\Mname| encapsulates two classes |Matrix|, |Vector| +and many functions of basic linear algebra. It is parametrized by a +number type |RT|. An instance of data type |Matrix| is a matrix of +variables of type |RT|, the so called ring type. Accordingly, +|Vector| implements vectors of variables of type |RT|. The arithmetic +type |RT| is required to behave like integers in the mathematical +sense. The manual pages of |Vector| and |Matrix| follow below. + +All functions compute the exact result, i.e., there is no rounding +error. Most functions of linear algebra are \emph{checkable}, i.e., +the programs can be asked for a proof that their output is +correct. For example, if the linear system solver declares a linear +system $A x = b$ unsolvable it also returns a vector $c$ such that +$c^T A = 0$ and $c^T b \neq 0$. All internal correctness checks can +be switched on by the flag [[CGAL_LA_SELFTEST]].}*/ + +public: + +/*{\Mtypes 5.5}*/ + +typedef _RT RT; +/*{\Mtypemember the ring type of the components.}*/ + +typedef CGAL::Ivector<_RT,_ALLOC> Vector; +/*{\Mtypemember the vector type.}*/ + +typedef CGAL::Imatrix<_RT,_ALLOC> Matrix; +/*{\Mtypemember the matrix type.}*/ + +typedef _ALLOC allocator_type; +/*{\Mtypemember the allocator used for memory management. |\Mname| is +an abbreviation for |Linear_algebraHd >|. Thus +|allocator_type| defaults to the standard allocator offered by the STL.}*/ + +/*{\Moperations 2 1}*/ + +static Matrix transpose(const Matrix& M); +/*{\Mstatic returns $M^T$ ($m\times n$ - matrix). }*/ + +static bool inverse(const Matrix& M, Matrix& I, RT& D, Vector& c); +/*{\Mstatic determines whether |M| has an inverse. It also computes +either the inverse as $(1/D) \cdot |I|$ or when no inverse +exists, a vector $c$ such that $c^T \cdot M = 0 $. }*/ + +static Matrix inverse(const Matrix& M, RT& D) +/*{\Mstatic returns the inverse matrix of |M|. More precisely, $1/D$ + times the matrix returned is the inverse of |M|.\\ + \precond |determinant(M) != 0|. }*/ +{ + Matrix result; + Vector c; + if (!inverse(M,result,D,c)) + ERROR_HANDLER(1,"Imatrix::inverse: matrix is singular."); + return result; +} + +static RT determinant (const Matrix& M, Matrix& L, Matrix& U, + std::vector& q, Vector& c); +/*{\Mstatic returns the determinant $D$ of |M| and sufficient information + to verify that the value of the determinant is correct. If + the determinant is zero then $c$ is a vector such that + $c^T \cdot M = 0$. If the determinant is non-zero then $L$ + and $U$ are lower and upper diagonal matrices respectively + and $q$ encodes a permutation matrix $Q$ with $Q(i,j) = 1$ + iff $i = q(j)$ such that $L \cdot M \cdot Q = U$, + $L(0,0) = 1$, $L(i,i) = U(i - 1,i - 1)$ for all $i$, + $1 \le i < n$, and $D = s \cdot U(n - 1,n - 1)$ where $s$ is + the determinant of $Q$. \precond |M| is square. }*/ + +static bool verify_determinant (const Matrix& M, RT D, Matrix& L, Matrix& U, + const std::vector& q, Vector& c); +/*{\Mstatic verifies the conditions stated above. }*/ + +static RT determinant (const Matrix& M); +/*{\Mstatic returns the determinant of |M|. + \precond |M| is square. }*/ + +static int sign_of_determinant (const Matrix& M); +/*{\Mstatic returns the sign of the determinant of |M|. + \precond |M| is square. }*/ + +static bool linear_solver(const Matrix& M, const Vector& b, + Vector& x, RT& D, + Matrix& spanning_vectors, + Vector& c); +/*{\Mstatic determines the complete solution space of the linear system + $M\cdot x = b$. If the system is unsolvable then + $c^T \cdot M = 0$ and $c^T \cdot b \not= 0$. + If the system is solvable then $(1/D) x$ is a solution, and + the columns of |spanning_vectors| are a maximal set of linearly + independent solutions to the corresponding homogeneous system. + \precond |M.row_dimension() = b.dimension()|. }*/ + +static bool linear_solver(const Matrix& M, const Vector& b, + Vector& x, RT& D, + Vector& c) +/*{\Mstatic determines whether the linear system $M\cdot x = b$ is + solvable. If yes, then $(1/D) x$ is a solution, if not then + $c^T \cdot M = 0$ and $c^T \cdot b \not= 0$. + \precond |M.row_dimension() = b.dimension()|. }*/ +{ + Matrix spanning_vectors; + return linear_solver(M,b,x,D,spanning_vectors,c); +} + +static bool linear_solver(const Matrix& M, const Vector& b, + Vector& x, RT& D) +/*{\Mstatic as above, but without the witness $c$ + \precond |M.row_dimension() = b.dimension()|. }*/ +{ + Matrix spanning_vectors; Vector c; + return linear_solver(M,b,x,D,spanning_vectors,c); +} + +static bool is_solvable(const Matrix& M, const Vector& b) +/*{\Mstatic determines whether the system $M \cdot x = b$ is solvable \\ + \precond |M.row_dimension() = b.dimension()|. }*/ +{ + Vector x; RT D; Matrix spanning_vectors; Vector c; + return linear_solver(M,b,x,D,spanning_vectors,c); +} + +static bool homogeneous_linear_solver (const Matrix& M, Vector& x); +/*{\Mstatic determines whether the homogeneous linear system + $M\cdot x = 0$ has a non - trivial solution. If + yes, then $x$ is such a solution. }*/ + +static int homogeneous_linear_solver (const Matrix& M, Matrix& spanning_vecs); +/*{\Mstatic determines the solution space of the homogeneous linear +system $M\cdot x = 0$. It returns the dimension of the solution space. +Moreover the columns of |spanning_vecs| span the solution space. }*/ + +static int independent_columns (const Matrix& M, std::vector& columns); +/*{\Mstatic returns the indices of a maximal subset of independent +columns of |M|.}*/ + +static int rank (const Matrix & M); +/*{\Mstatic returns the rank of matrix |M| }*/ + +/*{\Mimplementation The datatype |\Mname| is a wrapper class for the +linear algebra functionality on matrices and vectors. Operations +|determinant|, |inverse|, |linear_solver|, and |rank| take time +$O(n^3)$, and all other operations take time $O(nm)$. These time +bounds ignore the cost for multiprecision arithmetic operations. + +All functions on integer matrices compute the exact result, i.e., +there is no rounding error. The implemenation follows a proposal of +J. Edmonds (J. Edmonds, Systems of distinct representatives and linear +algebra, Journal of Research of the Bureau of National Standards, (B), +71, 241 - 245). Most functions of linear algebra are { \em checkable +}, i.e., the programs can be asked for a proof that their output is +correct. For example, if the linear system solver declares a linear +system $A x = b$ unsolvable it also returns a vector $c$ such that +$c^T A = 0$ and $c^T b \not= 0$.}*/ + +}; // Linear_algebraHd + + +CGAL_END_NAMESPACE + +#include +#endif // CGAL_LINALG_ALGEBRAHD_H + diff --git a/Packages/Kernel_d/test/Kernel_d/Linear_algebra-test.C b/Packages/Kernel_d/test/Kernel_d/Linear_algebra-test.C index 666545efd9a..dd26178497e 100644 --- a/Packages/Kernel_d/test/Kernel_d/Linear_algebra-test.C +++ b/Packages/Kernel_d/test/Kernel_d/Linear_algebra-test.C @@ -1,7 +1,9 @@ #include -#include +#include +#include #include #include +#include #define VEC_DIM 10 #define MAT_DIM 10 @@ -13,173 +15,356 @@ typedef leda_integer RT; #include typedef CGAL::Gmpz RT; #else -#include typedef double RT; #endif #endif -typedef CGAL::Linear_algebra LA; -typedef LA::Matrix Matrix; -typedef LA::Vector Vector; +#ifdef CGAL_USE_LEDA +#include +typedef leda_real FT; +#else +typedef double FT; +#endif int main(int argc, char* argv[]) { CGAL_TEST_START; - { - int vec_dim; - if (argc == 2) vec_dim = atoi(argv[1]); - else vec_dim = VEC_DIM; + { + typedef RT NT; + typedef CGAL::Linear_algebraHd LA; + typedef LA::Matrix Matrix; + typedef LA::Vector Vector; + bool IOTEST = true; + { + int vec_dim; + if (argc == 2) vec_dim = atoi(argv[1]); + else vec_dim = VEC_DIM; - /* some construction and access ops */ + /* some construction and access ops */ - Vector v0(vec_dim), v1(vec_dim), v2(vec_dim); - int F[] = { 1,2,3,4,5 }; - Vector v11(1,2), v12(1,2,3), v13(1,2,3,4), v14(F,F+5), - v15(v13.begin(),v13.end()), - v16(vec_dim,Vector::Initialize(),1); - CGAL_TEST(v13==v15); - CGAL_TEST(v0==v1); - for (int i = 0; i < vec_dim; i++) { - v1[i] = i; - v2[i] = vec_dim-i; - } - CGAL_TEST(v0!=v1); - Vector v3(v1); - CGAL_TEST(v3==v1); - - v1 += 2*v2; - v1 -= v2; - v1 = -v1; - CGAL_TEST((v1*v1 == vec_dim*vec_dim*vec_dim)); - /* squared length of (v1 + v2) = v1*v1 - should be equal to dim*dim*dim */ - RT res1 = (v1*v1 - v2*v2); - RT res2 = ((v1-v2)*(v1+v2)); - CGAL_TEST(res1==res2); - /* we test the third binomial formula */ - - CGAL_IO_TEST(v1,v2); - } - - { - int mat_dim; - if (argc == 2) mat_dim = atoi(argv[1]); - else mat_dim = MAT_DIM; - - /* some construction and access ops */ - Matrix A(mat_dim,mat_dim), B(mat_dim), - I(mat_dim,Matrix::Identity()), - One(mat_dim,mat_dim,Matrix::Initialize(),2); - CGAL_TEST(A==B); - CGAL_TEST(One*I==One); - std::vector F(mat_dim); - - int i,j; - for (i = 0; i < mat_dim; i++) { - Vector v(mat_dim); - for (j = 0; j < mat_dim; j++) { - A(i,j) = i; - B(i,j) = j; - v[j] = j; + Vector v0(vec_dim), v1(vec_dim), v2(vec_dim); + int F[] = { 1,2,3,4,5 }; + Vector v11(1,2), v12(1,2,3), v13(1,2,3,4), v14(F,F+5), + v15(v13.begin(),v13.end()), + v16(vec_dim,Vector::Initialize(),1); + CGAL_TEST(v13==v15); + CGAL_TEST(v0==v1); + for (int i = 0; i < vec_dim; i++) { + v1[i] = i; + v2[i] = vec_dim-i; } - F[i] = v; - } - Matrix C(F), D(A), E(F.begin(),F.end()); - CGAL_TEST(C==F && E==F); - /* A = (0,0,0,0...) - (1,1,1,1...) - (2,2,2,2...) - ... - B = (0,1,2,3...) - (0,1,2,3...) - (0,1,2,3...) - ... - C = A = D = E; - */ + CGAL_TEST(v0!=v1); + Vector v3(v1); + CGAL_TEST(v3==v1); - CGAL_TEST(A==C); - CGAL_TEST(A!=B); - CGAL_TEST(A==D); - CGAL_TEST(A.row_dimension()==mat_dim && A.column_dimension()==mat_dim); - CGAL_TEST(A.row(1)*B.column(1)==mat_dim); - CGAL_IO_TEST(A,C); + v1 += 2*v2; + v1 -= v2; + v1 = -v1; + CGAL_TEST((v1*v1 == vec_dim*vec_dim*vec_dim)); + /* squared length of (v1 + v2) = v1*v1 + should be equal to dim*dim*dim */ + NT res1 = (v1*v1 - v2*v2); + NT res2 = ((v1-v2)*(v1+v2)); + CGAL_TEST(res1==res2); + /* we test the third binomial formula */ - /* some basic arithmetic testing */ - C += A; - C -= 3*A; - C = -C; - CGAL_TEST(C==A); - - /* row sum test: */ - Vector ones(mat_dim), row_sum_vec(mat_dim); - for (i = 0; i < mat_dim; i++) { - ones[i] = 1; - row_sum_vec[i] = i*mat_dim; - } - CGAL_TEST(A*ones==row_sum_vec); - - C = A+B; - /* matrix operations + ,* and |transpose|, |rank| */ - CGAL_TEST(C==LA::transpose(C)); - C = C-B; - CGAL_TEST(C==A); - - C = I+A; - CGAL_TEST(LA::rank(C)==mat_dim); - - /* matrix operations 2* , |determinant| */ - C = 2 * I; - C = C * 2; - Matrix L,U; - Vector c; - std::vector q; - RT det = LA::determinant(C, L, U, q, c); // det must be 2^{2*mat_dim} - CGAL_TEST(log(det)==2*mat_dim); - CGAL_TEST(LA::verify_determinant(C, det, L, U, q, c)); - CGAL_TEST(det == LA::determinant(C)); - CGAL_TEST(sign(det) == LA::sign_of_determinant(C)); - LA::independent_columns(C,q); - - /* a random linear solver task: */ - Vector b(mat_dim),x(mat_dim),e; - RT denom; - for (i = 0; i < mat_dim; i++) { - for (j = 0; j < mat_dim; j++) - E(i,j) = CGAL::default_random.get_int( - mat_dim,mat_dim); - b[i] = CGAL::default_random.get_int( - mat_dim,mat_dim); + if (IOTEST) CGAL_IO_TEST(v1,v2); + CGAL::set_pretty_mode ( cerr ); } - // linear solver - if (LA::linear_solver(E, b, x, denom, A, e)) { - CGAL_TEST(E*x == b*denom); - LA::linear_solver(E,b,x,denom,e); - CGAL_TEST(E*x == b*denom); - LA::linear_solver(E,b,x,denom); - CGAL_TEST(E*x == b*denom); - CGAL_TEST(LA::is_solvable(E,b)); - } - else { - Vector null(mat_dim); - CGAL_TEST(LA::transpose(E)*e==null && b*e!=0); - CGAL_TEST(!LA::is_solvable(E,b)); + { + int mat_dim; + if (argc == 2) mat_dim = atoi(argv[1]); + else mat_dim = MAT_DIM; + + /* some construction and access ops */ + Matrix A(mat_dim,mat_dim), B(mat_dim), + I(mat_dim,Matrix::Identity()), + One(mat_dim,mat_dim,Matrix::Initialize(),2); + CGAL_TEST(A==B); + CGAL_TEST(One*I==One); + std::vector F(mat_dim); + + int i,j; + for (i = 0; i < mat_dim; i++) { + Vector v(mat_dim); + for (j = 0; j < mat_dim; j++) { + A(i,j) = i; + B(i,j) = j; + v[j] = j; + } + F[i] = v; + } + Matrix C(F), D(A), E(F.begin(),F.end()); + CGAL_TEST(C==F && E==F); + /* A = (0,0,0,0...) + (1,1,1,1...) + (2,2,2,2...) + ... + B = (0,1,2,3...) + (0,1,2,3...) + (0,1,2,3...) + ... + C = A = D = E; + */ + + CGAL_TEST(A==C); + CGAL_TEST(A!=B); + CGAL_TEST(A==D); + CGAL_TEST(A.row_dimension()==mat_dim && A.column_dimension()==mat_dim); + CGAL_TEST(A.row(1)*B.column(1)==mat_dim); + + /* some basic arithmetic testing */ + C = A; C += A; + C -= 3*A; + C = -C; + CGAL_TEST(C==A); + + /* row sum test: */ + Vector ones(mat_dim), row_sum_vec(mat_dim); + for (i = 0; i < mat_dim; i++) { + ones[i] = 1; + row_sum_vec[i] = i*mat_dim; + } + CGAL_TEST(A*ones==row_sum_vec); + + C = A+B; + /* matrix operations + ,* and |transpose|, |rank| */ + CGAL_TEST(C==LA::transpose(C)); + C = C-B; + CGAL_TEST(C==A); + C = I+A; + CGAL_TEST(LA::rank(C)==mat_dim); + + /* matrix operations 2* , |determinant| */ + C = 2 * I; + C = C * 2; + Matrix L,U; + Vector c; + std::vector q; + NT det = LA::determinant(C, L, U, q, c); // det must be 2^{2*mat_dim} + NT pot = 1; for (int i=0; i 0.0; f = f-=0.1) { + Matrix E = C; + for (int i = 0; i< mat_dim; ++i) + for (int j = 0; j < mat_dim; ++j) + if ( CGAL::default_random.get_double() > f ) E(i,j)=0; + + if (LA::linear_solver(E, b, x, denom, A, e)) { + CGAL_TEST(E*x == b*denom); + LA::linear_solver(E,b,x,denom,e); + CGAL_TEST(E*x == b*denom); + LA::linear_solver(E,b,x,denom); + CGAL_TEST(E*x == b*denom); + CGAL_TEST(LA::is_solvable(E,b)); + } + else { + Vector null(mat_dim); + CGAL_TEST(LA::transpose(E)*e==null && b*e!=0); + CGAL_TEST(!LA::is_solvable(E,b)); + } + + Matrix SV; + if (LA::homogeneous_linear_solver(E,x)) { + CGAL_TEST(E*x==Vector(mat_dim)); + CGAL_TEST(LA::homogeneous_linear_solver(E,SV)==LA::rank(SV)); + } + + // inverse + if (LA::inverse(E,D,denom,c)) { + CGAL_TEST(E*D==denom*Matrix(mat_dim,Matrix::Identity())); + CGAL_TEST(D==LA::inverse(E,denom)); + } else { + CGAL_TEST(LA::transpose(E)*c==Vector(mat_dim)); + } + LA::independent_columns(E,q); + } } - Matrix SV; - if (LA::homogeneous_linear_solver(E,x)) { - CGAL_TEST(E*x==Vector(mat_dim)); - CGAL_TEST(LA::homogeneous_linear_solver(E,SV)==LA::rank(SV)); - } - // inverse - if (LA::inverse(E,D,denom,c)) { - CGAL_TEST(E*D==denom*Matrix(mat_dim,Matrix::Identity())); - CGAL_TEST(D==LA::inverse(E,denom)); - } else { - CGAL_TEST(LA::transpose(E)*c==Vector(mat_dim)); - } } + { + typedef FT NT; + typedef CGAL::Linear_algebraCd LA; + typedef LA::Matrix Matrix; + typedef LA::Vector Vector; + bool IOTEST = false; + { + int vec_dim; + if (argc == 2) vec_dim = atoi(argv[1]); + else vec_dim = VEC_DIM; + + /* some construction and access ops */ + + Vector v0(vec_dim), v1(vec_dim), v2(vec_dim); + int F[] = { 1,2,3,4,5 }; + Vector v11(1,2), v12(1,2,3), v13(1,2,3,4), v14(F,F+5), + v15(v13.begin(),v13.end()), + v16(vec_dim,Vector::Initialize(),1); + CGAL_TEST(v13==v15); + CGAL_TEST(v0==v1); + for (int i = 0; i < vec_dim; i++) { + v1[i] = i; + v2[i] = vec_dim-i; + } + CGAL_TEST(v0!=v1); + Vector v3(v1); + CGAL_TEST(v3==v1); + + v1 += 2*v2; + v1 -= v2; + v1 = -v1; + CGAL_TEST((v1*v1 == vec_dim*vec_dim*vec_dim)); + /* squared length of (v1 + v2) = v1*v1 + should be equal to dim*dim*dim */ + NT res1 = (v1*v1 - v2*v2); + NT res2 = ((v1-v2)*(v1+v2)); + CGAL_TEST(res1==res2); + /* we test the third binomial formula */ + + if (IOTEST) CGAL_IO_TEST(v1,v2); + CGAL::set_pretty_mode ( cerr ); + } + + { + int mat_dim; + if (argc == 2) mat_dim = atoi(argv[1]); + else mat_dim = MAT_DIM; + + /* some construction and access ops */ + Matrix A(mat_dim,mat_dim), B(mat_dim), + I(mat_dim,Matrix::Identity()), + One(mat_dim,mat_dim,Matrix::Initialize(),2); + CGAL_TEST(A==B); + CGAL_TEST(One*I==One); + std::vector F(mat_dim); + + int i,j; + for (i = 0; i < mat_dim; i++) { + Vector v(mat_dim); + for (j = 0; j < mat_dim; j++) { + A(i,j) = i; + B(i,j) = j; + v[j] = j; + } + F[i] = v; + } + Matrix C(F), D(A), E(F.begin(),F.end()); + CGAL_TEST(C==F && E==F); + /* A = (0,0,0,0...) + (1,1,1,1...) + (2,2,2,2...) + ... + B = (0,1,2,3...) + (0,1,2,3...) + (0,1,2,3...) + ... + C = A = D = E; + */ + + CGAL_TEST(A==C); + CGAL_TEST(A!=B); + CGAL_TEST(A==D); + CGAL_TEST(A.row_dimension()==mat_dim && A.column_dimension()==mat_dim); + CGAL_TEST(A.row(1)*B.column(1)==mat_dim); + + /* some basic arithmetic testing */ + C = A; C += A; + C -= 3*A; + C = -C; + CGAL_TEST(C==A); + + /* row sum test: */ + Vector ones(mat_dim), row_sum_vec(mat_dim); + for (i = 0; i < mat_dim; i++) { + ones[i] = 1; + row_sum_vec[i] = i*mat_dim; + } + CGAL_TEST(A*ones==row_sum_vec); + + C = A+B; + /* matrix operations + ,* and |transpose|, |rank| */ + CGAL_TEST(C==LA::transpose(C)); + C = C-B; + CGAL_TEST(C==A); + C = I+A; + CGAL_TEST(LA::rank(C)==mat_dim); + + /* matrix operations 2* , |determinant| */ + C = 2 * I; + C = C * 2; + Matrix L,U; + Vector c; + std::vector q; + NT det = LA::determinant(C, L, U, q, c); // det must be 2^{2*mat_dim} + NT pot = 1; for (int i=0; i 0.0; f = f-=0.1) { + Matrix E = C; + for (int i = 0; i< mat_dim; ++i) + for (int j = 0; j < mat_dim; ++j) + if ( CGAL::default_random.get_double() > f ) E(i,j)=0; + + if (LA::linear_solver(E, b, x, denom, A, e)) { + CGAL_TEST(E*x == b*denom); + LA::linear_solver(E,b,x,denom,e); + CGAL_TEST(E*x == b*denom); + LA::linear_solver(E,b,x,denom); + CGAL_TEST(E*x == b*denom); + CGAL_TEST(LA::is_solvable(E,b)); + } + else { + Vector null(mat_dim); + CGAL_TEST(LA::transpose(E)*e==null && b*e!=0); + CGAL_TEST(!LA::is_solvable(E,b)); + } + + Matrix SV; + if (LA::homogeneous_linear_solver(E,x)) { + CGAL_TEST(E*x==Vector(mat_dim)); + CGAL_TEST(LA::homogeneous_linear_solver(E,SV)==LA::rank(SV)); + } + + // inverse + if (LA::inverse(E,D,denom,c)) { + CGAL_TEST(E*D==denom*Matrix(mat_dim,Matrix::Identity())); + CGAL_TEST(D==LA::inverse(E,denom)); + } else { + CGAL_TEST(LA::transpose(E)*c==Vector(mat_dim)); + } + LA::independent_columns(E,q); + } + } + } CGAL_TEST_END; } diff --git a/Packages/Kernel_d/test/Kernel_d/Linear_algebra-test.cmd b/Packages/Kernel_d/test/Kernel_d/Linear_algebra-test.cmd new file mode 100644 index 00000000000..7ed6ff82de6 --- /dev/null +++ b/Packages/Kernel_d/test/Kernel_d/Linear_algebra-test.cmd @@ -0,0 +1 @@ +5