12345qiupeng
/
gtsam
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Changed problematic blas call to ublas - was sometimes producing incorrect results but I do not know why :-(

release/4.3a0
Richard Roberts 2011-02-04 00:00:31 +00:00
parent 654f4da4e5
commit 62aa7a6681
1 changed files with 4 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
gtsam/base/cholesky.cpp
View File

@ -21,6 +21,7 @@
#include <gtsam/base/timing.h> #include <gtsam/base/timing.h>
#include <boost/numeric/ublas/matrix_proxy.hpp> #include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/blas.hpp> #include <boost/numeric/ublas/blas.hpp>
#include <boost/format.hpp> #include <boost/format.hpp>
@ -254,7 +255,8 @@ void choleskyPartial(MatrixColMajor& ABC, size_t nFrontal) {
// Views of R, S, and L submatrices. // Views of R, S, and L submatrices.
ublas::matrix_range<MatrixColMajor> R(ublas::project(ABC, ublas::range(0,nFrontal), ublas::range(0,nFrontal))); ublas::matrix_range<MatrixColMajor> R(ublas::project(ABC, ublas::range(0,nFrontal), ublas::range(0,nFrontal)));
ublas::matrix_range<MatrixColMajor> S(ublas::project(ABC, ublas::range(0,nFrontal), ublas::range(nFrontal,n))); ublas::matrix_range<MatrixColMajor> S(ublas::project(ABC, ublas::range(0,nFrontal), ublas::range(nFrontal,n)));
ublas::matrix_range<MatrixColMajor> L(ublas::project(ABC, ublas::range(nFrontal,n), ublas::range(nFrontal,n))); ublas::matrix_range<MatrixColMajor> Lf(ublas::project(ABC, ublas::range(nFrontal,n), ublas::range(nFrontal,n)));
ublas::symmetric_adaptor<typeof(Lf), ublas::upper> L(Lf);
// Compute S = inv(R') * B // Compute S = inv(R') * B
tic(2, "compute S"); tic(2, "compute S");
@ -267,7 +269,7 @@ void choleskyPartial(MatrixColMajor& ABC, size_t nFrontal) {
tic(3, "compute L"); tic(3, "compute L");
if(debug) gtsam::print(L, "C: "); if(debug) gtsam::print(L, "C: ");
if(L.size2() > 0) if(L.size2() > 0)
cblas_dsyrk(CblasColMajor, CblasUpper, CblasTrans, L.size1(), S.size1(), -1.0, &S(0,0), n, 1.0, &L(0,0), n); L -= ublas::prod(ublas::trans(S), S);
if(debug) gtsam::print(L, "L: "); if(debug) gtsam::print(L, "L: ");
toc(3, "compute L"); toc(3, "compute L");
} }