From 4501a57ae3dc87edc9cee1cde7f631ee1f1ed161 Mon Sep 17 00:00:00 2001 From: Marc Pouget Date: Wed, 3 May 2006 17:46:35 +0000 Subject: [PATCH] debug de lapack... --- .../examples/Jet_fitting_3/LinAlg_lapack.h | 43 ++++++++++++++----- .../include/CGAL/Monge_via_jet_fitting.h | 4 +- 2 files changed, 36 insertions(+), 11 deletions(-) diff --git a/Jet_fitting_3/examples/Jet_fitting_3/LinAlg_lapack.h b/Jet_fitting_3/examples/Jet_fitting_3/LinAlg_lapack.h index 5d3758b1eb2..d1b9731b202 100644 --- a/Jet_fitting_3/examples/Jet_fitting_3/LinAlg_lapack.h +++ b/Jet_fitting_3/examples/Jet_fitting_3/LinAlg_lapack.h @@ -66,13 +66,27 @@ void GSL::eigen_symm_algo(Matrix& S, double* eval, Matrix& evec) liwork = 10*n, info; double vl=0, vu=0, abstol = 0; - integer* isuppz = (integer*) malloc(n*sizeof(integer)); - double* work = (double*) malloc(lwork*sizeof(double)); - integer* iwork = (integer*) malloc(liwork*sizeof(integer)); +/* //c style */ +/* integer* isuppz = (integer*) malloc(n*sizeof(integer)); */ +/* double* work = (double*) malloc(lwork*sizeof(double)); */ +/* integer* iwork = (integer*) malloc(liwork*sizeof(integer)); */ + //c++ style + integer* isuppz; + isuppz = new integer [n]; + double* work = new double [lwork]; + integer* iwork = new integer [liwork]; 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); */ +/* //C++ style */ + delete [] isuppz; + delete [] work; + delete [] iwork; } void GSL::solve_ls_svd_algo(Matrix& M, double* B, double &cond_nb) @@ -85,18 +99,27 @@ void GSL::solve_ls_svd_algo(Matrix& M, double* B, double &cond_nb) rank, lwork = 5*m, info; - double* sing_values = (double*) malloc(n*sizeof(double)); - double* work = (double*) malloc(lwork*sizeof(double)); +/* //c style */ +/* double* sing_values = (double*) malloc(n*sizeof(double)); */ +/* double* work = (double*) malloc(lwork*sizeof(double)); */ + //c++ style + double* sing_values = new double [n]; + double* work = new double [lwork]; + double rcond = -1; - dgelss_(&m, &n, &nrhs, M.matrix(), &lda, B, &ldb, sing_values, &rcond, &rank, work, &lwork, &info); + 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); +/* //clean up */ +/* free(sing_values); */ +/* free(work); */ +//c++ style + delete [] sing_values; + delete [] work; } #endif 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 ae7a8ed0bc2..69e6c494468 100644 --- a/Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h +++ b/Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h @@ -349,7 +349,7 @@ compute_PCA(Range_Iterator begin, Range_Iterator end, Cov.set_elt(1,1,yy); Cov.set_elt(1,2,yz); Cov.set_elt(2,0,xz); - Cov.set_elt(2,2,yz); + Cov.set_elt(2,1,yz); Cov.set_elt(2,2,zz); // solve for eigenvalues and eigenvectors. // eigen values are sorted in ascending order, @@ -360,9 +360,11 @@ compute_PCA(Range_Iterator begin, Range_Iterator end, monge_info.pca_eigen_vals()[0] = eval[2];//implicit cast LAFT->LFT LVector temp_vectn(evec.get_elt(0,2),evec.get_elt(1,2),evec.get_elt(2,2)); monge_info.pca_eigen_vecs()[0] = temp_vectn; + monge_info.pca_eigen_vals()[1] = eval[1]; LVector temp_vect1(evec.get_elt(0,1),evec.get_elt(1,1),evec.get_elt(2,1)); monge_info.pca_eigen_vecs()[1] = temp_vect1; + monge_info.pca_eigen_vals()[2] = eval[0]; LVector temp_vect2(evec.get_elt(0,0),evec.get_elt(1,0),evec.get_elt(2,0)); monge_info.pca_eigen_vecs()[2] = temp_vect2;