debug de lapack...

This commit is contained in:
Marc Pouget 2006-05-03 17:46:35 +00:00
parent b9ef01ba52
commit 4501a57ae3
2 changed files with 36 additions and 11 deletions

View File

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

View File

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