diff --git a/Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h b/Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h index 6b327dd408f..a1ce2131951 100644 --- a/Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h +++ b/Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h @@ -153,7 +153,7 @@ dump_4ogl(std::ostream& out_stream, const FT scale) template std::ostream& -operator<<(std::ostream& out_stream, const Monge_form& monge) +operator<<(std::ostream& out_stream, Monge_form& monge) { monge.dump_verbose(out_stream); return out_stream; diff --git a/Jet_fitting_3/test/Jet_fitting_3/blind_1pt.C b/Jet_fitting_3/test/Jet_fitting_3/blind_1pt.C index 3d4615d763c..09bedb86ddf 100644 --- a/Jet_fitting_3/test/Jet_fitting_3/blind_1pt.C +++ b/Jet_fitting_3/test/Jet_fitting_3/blind_1pt.C @@ -7,21 +7,19 @@ #include #include -#include "../../include/CGAL/Monge_via_jet_fitting.h" -#include "include/LinAlg_lapack.h" +#include typedef double DFT; typedef CGAL::Cartesian Data_Kernel; typedef Data_Kernel::Point_3 DPoint; typedef Data_Kernel::Vector_3 DVector; -typedef CGAL::Monge_rep My_Monge_rep; typedef double LFT; typedef CGAL::Cartesian Local_Kernel; -typedef CGAL::Monge_info My_Monge_info; -typedef CGAL::Monge_via_jet_fitting My_Monge_via_jet_fitting; - - +typedef CGAL::Monge_via_jet_fitting My_Monge_via_jet_fitting; +typedef My_Monge_via_jet_fitting::Monge_form My_Monge_form; +typedef My_Monge_via_jet_fitting::Monge_form_condition_numbers My_Monge_form_condition_numbers; + int main() { //open the input file @@ -49,17 +47,23 @@ int main() // fct parameters int d_fitting = 4; int d_monge = 4; - My_Monge_rep monge_rep; - My_Monge_info monge_info; + My_Monge_form monge_form; + My_Monge_form_condition_numbers monge_form_condition_numbers; //run the main fct My_Monge_via_jet_fitting do_it(in_points.begin(), in_points.end(), d_fitting, d_monge, - monge_rep, monge_info); + monge_form, monge_form_condition_numbers); - monge_rep.comply_wrt_given_normal( -monge_rep.n() ); + monge_form.comply_wrt_given_normal( -monge_form.n() ); //OUTPUT on std::cout - monge_rep.dump_verbose( std::cout ); - monge_rep.dump_4ogl( std::cout, 1 ); - monge_info.dump_verbose( std::cout ); + std::cout << monge_form + << monge_form_condition_numbers; + + monge_form.dump_4ogl( std::cout, 1 ); + double precision = 0.01; + assert(monge_form.coefficients()[0] >= -0.2 - precision); + assert(monge_form.coefficients()[0] <= -0.2 + precision); + assert(monge_form.coefficients()[1] >= -0.4 - precision); + assert(monge_form.coefficients()[1] <= -0.4 + precision); std::cout << "success\n"; } diff --git a/Jet_fitting_3/test/Jet_fitting_3/include/LinAlg_lapack.h b/Jet_fitting_3/test/Jet_fitting_3/include/LinAlg_lapack.h deleted file mode 100644 index 1f9d6cc9dc3..00000000000 --- a/Jet_fitting_3/test/Jet_fitting_3/include/LinAlg_lapack.h +++ /dev/null @@ -1,110 +0,0 @@ -#ifndef _LAPACK_H_ -#define _LAPACK_H_ - -#include -#include "blaswrap.h" -#include "f2c.h" -extern "C" { -#include "clapack.h" -} - -////////////////////////class Lapack_matrix///////////////////// -//in Lapack, matrices are one-dimensional arrays -// and elements are column-major ordered -// this class is a wrapper defining set and get in the usual way -class Lapack_matrix{ -protected: - double* m_matrix; -public: - size_t nb_rows; - size_t nb_columns; - //contructor - // initializes all the elements of the matrix to zero. - Lapack_matrix(size_t n1, size_t n2) { - m_matrix = (double*) calloc (n1*n2, sizeof(double)); - nb_rows = n1; - nb_columns = n2; - } - - //destructor - // ~Lapack_matrix();// { free m_matrix; } - - //access - const double* matrix() const { return m_matrix; } - double* matrix() { return m_matrix; } - - void set_elt(size_t i, size_t j, const double value) { m_matrix[j*nb_rows+i] = value; } - double get_elt(size_t i, size_t j) { return m_matrix[j*nb_rows+i]; } -}; - -//Lapack_matrix::~Lapack_matrix() { free m_matrix; } bug! - -////////////////////////class Lapack///////////////////// -class Lapack{ -public: - typedef Lapack_matrix Matrix; - // solve for eigenvalues and eigenvectors. - // eigen values are sorted in ascending order, - // eigen vectors are sorted in accordance. - static - void eigen_symm_algo(Matrix& S, double* eval, Matrix& evec); - //solve MX=B using SVD and give the condition number of M - static - void solve_ls_svd_algo(Matrix& M, double* B, double &cond_nb); -}; - -void Lapack::eigen_symm_algo(Matrix& S, double* eval, Matrix& evec) -{ - char jobz = 'V'; - char range = 'A'; - char uplo = 'U'; - integer n = S.nb_rows, - lda = S.nb_rows, - il = 0, iu = 0, - m = n, - lwork = 26*n, - liwork = 10*n, - info; - double vl=0, vu=0, abstol = 0; - - integer* isuppz = (integer*) malloc(2*n*sizeof(integer)); - double* work = (double*) malloc(lwork*sizeof(double)); - integer* iwork = (integer*) malloc(liwork*sizeof(integer)); - - dsyevr_(&jobz, &range, &uplo, &n, S.matrix(), &lda, &vl, &vu, &il, &iu, &abstol, &m, - eval, evec.matrix(), &n, isuppz, work, &lwork, iwork, &liwork, &info); - - //clean up - free(isuppz); - free(work); - free(iwork); -} - -void Lapack::solve_ls_svd_algo(Matrix& M, double* B, double &cond_nb) -{ - integer m = M.nb_rows, - n = M.nb_columns, - nrhs = 1, - lda = m, - ldb = m, - rank, - lwork = 5*m, - info; - //c style - double* sing_values = (double*) malloc(n*sizeof(double)); - double* work = (double*) malloc(lwork*sizeof(double)); - - double rcond = -1; - - dgelss_(&m, &n, &nrhs, M.matrix(), &lda, B, &ldb, sing_values, - &rcond, &rank, work, &lwork, &info); - assert(info==0); - - cond_nb = sing_values[0]/sing_values[n-1]; - - //clean up - free(sing_values); - free(work); -} - -#endif diff --git a/Jet_fitting_3/test/Jet_fitting_3/makefile b/Jet_fitting_3/test/Jet_fitting_3/makefile index de4bc9e4e7e..af78eee1667 100644 --- a/Jet_fitting_3/test/Jet_fitting_3/makefile +++ b/Jet_fitting_3/test/Jet_fitting_3/makefile @@ -48,15 +48,12 @@ all: \ blind_1pt blind_1pt: blind_1pt.o - $(CGAL_CXX) $(LIBPATH) -o blind_1pt.exe blind_1pt.o \ + $(CGAL_CXX) $(LIBPATH) -o blind_1pt blind_1pt.o \ $(LDFLAGS) -rmexe: - rm *.exe - clean: \ blind_1pt.clean - rmexe + depend: @@ -84,58 +81,4 @@ blind_1pt.o: /usr/include/bits/endian.h /usr/include/sys/select.h blind_1pt.o: /usr/include/bits/select.h /usr/include/bits/sigset.h blind_1pt.o: /usr/include/bits/time.h /usr/include/sys/sysmacros.h blind_1pt.o: /usr/include/bits/pthreadtypes.h /usr/include/bits/sched.h -blind_1pt.o: /usr/include/alloca.h ../../include/CGAL/Monge_via_jet_fitting.h -blind_1pt.o: ../../include/CGAL/jet_fitting_3_assertions.h -blind_1pt.o: /usr/include/math.h /usr/include/bits/huge_val.h -blind_1pt.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h -blind_1pt.o: include/GSL.h /usr/include/gsl/gsl_linalg.h -blind_1pt.o: /usr/include/gsl/gsl_mode.h /usr/include/gsl/gsl_permutation.h -blind_1pt.o: /usr/include/gsl/gsl_types.h /usr/include/gsl/gsl_errno.h -blind_1pt.o: /usr/include/errno.h /usr/include/bits/errno.h -blind_1pt.o: /usr/include/linux/errno.h /usr/include/asm/errno.h -blind_1pt.o: /usr/include/asm-i386/errno.h /usr/include/gsl/gsl_check_range.h -blind_1pt.o: /usr/include/gsl/gsl_vector.h -blind_1pt.o: /usr/include/gsl/gsl_vector_complex_long_double.h -blind_1pt.o: /usr/include/gsl/gsl_complex.h -blind_1pt.o: /usr/include/gsl/gsl_vector_long_double.h -blind_1pt.o: /usr/include/gsl/gsl_block_long_double.h -blind_1pt.o: /usr/include/gsl/gsl_vector_complex.h -blind_1pt.o: /usr/include/gsl/gsl_block_complex_long_double.h -blind_1pt.o: /usr/include/gsl/gsl_vector_complex_double.h -blind_1pt.o: /usr/include/gsl/gsl_vector_double.h -blind_1pt.o: /usr/include/gsl/gsl_block_double.h -blind_1pt.o: /usr/include/gsl/gsl_block_complex_double.h -blind_1pt.o: /usr/include/gsl/gsl_vector_complex_float.h -blind_1pt.o: /usr/include/gsl/gsl_vector_float.h -blind_1pt.o: /usr/include/gsl/gsl_block_float.h -blind_1pt.o: /usr/include/gsl/gsl_block_complex_float.h -blind_1pt.o: /usr/include/gsl/gsl_vector_ulong.h -blind_1pt.o: /usr/include/gsl/gsl_block_ulong.h -blind_1pt.o: /usr/include/gsl/gsl_vector_long.h -blind_1pt.o: /usr/include/gsl/gsl_block_long.h -blind_1pt.o: /usr/include/gsl/gsl_vector_uint.h -blind_1pt.o: /usr/include/gsl/gsl_block_uint.h -blind_1pt.o: /usr/include/gsl/gsl_vector_int.h -blind_1pt.o: /usr/include/gsl/gsl_block_int.h -blind_1pt.o: /usr/include/gsl/gsl_vector_ushort.h -blind_1pt.o: /usr/include/gsl/gsl_block_ushort.h -blind_1pt.o: /usr/include/gsl/gsl_vector_short.h -blind_1pt.o: /usr/include/gsl/gsl_block_short.h -blind_1pt.o: /usr/include/gsl/gsl_vector_uchar.h -blind_1pt.o: /usr/include/gsl/gsl_block_uchar.h -blind_1pt.o: /usr/include/gsl/gsl_vector_char.h -blind_1pt.o: /usr/include/gsl/gsl_block_char.h /usr/include/gsl/gsl_matrix.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_complex_long_double.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_complex_double.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_complex_float.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_long_double.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_double.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_float.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_ulong.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_long.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_uint.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_int.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_ushort.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_short.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_uchar.h -blind_1pt.o: /usr/include/gsl/gsl_matrix_char.h /usr/include/gsl/gsl_eigen.h +blind_1pt.o: /usr/include/alloca.h