This commit is contained in:
Marc Pouget 2006-06-28 21:42:42 +00:00
parent 05c4292350
commit fbae35c971
4 changed files with 22 additions and 185 deletions

View File

@ -153,7 +153,7 @@ dump_4ogl(std::ostream& out_stream, const FT scale)
template <class DataKernel>
std::ostream&
operator<<(std::ostream& out_stream, const Monge_form<DataKernel>& monge)
operator<<(std::ostream& out_stream, Monge_form<DataKernel>& monge)
{
monge.dump_verbose(out_stream);
return out_stream;

View File

@ -7,21 +7,19 @@
#include <stdlib.h>
#include <vector>
#include "../../include/CGAL/Monge_via_jet_fitting.h"
#include "include/LinAlg_lapack.h"
#include <CGAL/Monge_via_jet_fitting.h>
typedef double DFT;
typedef CGAL::Cartesian<DFT> Data_Kernel;
typedef Data_Kernel::Point_3 DPoint;
typedef Data_Kernel::Vector_3 DVector;
typedef CGAL::Monge_rep<Data_Kernel> My_Monge_rep;
typedef double LFT;
typedef CGAL::Cartesian<LFT> Local_Kernel;
typedef CGAL::Monge_info<Local_Kernel> My_Monge_info;
typedef CGAL::Monge_via_jet_fitting<Data_Kernel, Local_Kernel, Lapack> My_Monge_via_jet_fitting;
typedef CGAL::Monge_via_jet_fitting<Data_Kernel> 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";
}

View File

@ -1,110 +0,0 @@
#ifndef _LAPACK_H_
#define _LAPACK_H_
#include <stdlib.h>
#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

View File

@ -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