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

Integrated blas into gtsam using autotools

release/4.3a0
Alex Cunningham 2010-03-15 18:17:43 +00:00
parent 3a5aeeeb0d
commit a1273a91fc
5 changed files with 91 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
configure.ac
View File

@ -26,30 +26,44 @@ AC_ARG_ENABLE([debug],
AM_CONDITIONAL([DEBUG], [test x$debug = xtrue]) AM_CONDITIONAL([DEBUG], [test x$debug = xtrue])
# search for gsl # search for gsl
AM_PATH_GSL(1.1) #AM_PATH_GSL(1.1)
# enable using GSL for linalg # enable using GSL for linalg
AC_ARG_ENABLE([gsl], #AC_ARG_ENABLE([gsl],
[ --enable-gsl Enable the GSL library], # [ --enable-gsl Enable the GSL library],
[case "${enableval}" in # [case "${enableval}" in
yes) gsl=true ;; # yes) gsl=true ;;
no) gsl=false ;; # no) gsl=false ;;
*) AC_MSG_ERROR([bad value ${enableval} for --enable-gsl]) ;; # *) AC_MSG_ERROR([bad value ${enableval} for --enable-gsl]) ;;
esac],[gsl=false]) # esac],[gsl=false])
#
AM_CONDITIONAL([GSL], [test x$gsl = xtrue]) #AM_CONDITIONAL([GSL], [test x$gsl = xtrue])
# enable using ATLAS for BLAS # enable using ATLAS for BLAS
AC_ARG_ENABLE([atlas], #AC_ARG_ENABLE([atlas],
[ --enable-atlas Enable ATLAS optimized BLAS], # [ --enable-atlas Enable ATLAS optimized BLAS],
[case "${enableval}" in # [case "${enableval}" in
yes) atlas=true ;; # yes) atlas=true ;;
no) atlas=false ;; # no) atlas=false ;;
*) AC_MSG_ERROR([bad value ${enableval} for --enable-atlas]) ;; # *) AC_MSG_ERROR([bad value ${enableval} for --enable-atlas]) ;;
esac],[atlas=false]) # esac],[atlas=false])
#
#AM_CONDITIONAL([ATLAS], [test x$atlas = xtrue])
AM_CONDITIONAL([ATLAS], [test x$atlas = xtrue]) # search for a blas implementation
AX_BLAS()
# enable BLAS with general purpose script
AC_ARG_ENABLE([blas],
[ --enable-blas Enable external BLAS library],
[case "${enableval}" in
yes) blas=true ;;
no) blas=false ;;
*) AC_MSG_ERROR([bad value ${enableval} for --enable-blas]) ;;
esac],[blas=false])
AM_CONDITIONAL([USE_BLAS], [test x$blas = xtrue])
# Checks for programs. # Checks for programs.
AC_PROG_CXX AC_PROG_CXX

27
cpp/Makefile.am
View File

@ -293,27 +293,14 @@ include_HEADERS = $(headers)
AM_CXXFLAGS += -I.. AM_CXXFLAGS += -I..
AM_LDFLAGS = -L../CppUnitLite -lCppUnitLite $(BOOST_LDFLAGS) $(boost_serialization) #-L. -lgtsam AM_LDFLAGS = -L../CppUnitLite -lCppUnitLite $(BOOST_LDFLAGS) $(boost_serialization) #-L. -lgtsam
# TESTING: adding cblas implementation with atlas # adding cblas implementation with atlas
#AM_CXXFLAGS += -DCBLAS if USE_BLAS
#libgtsam_la_CPPFLAGS += -DCBLAS AM_CXXFLAGS += -DCBLAS
#AM_LDFLAGS += -L/opt/local/lib -lcblas -latlas -llapack -lptcblas -lptf77blas -lf77blas libgtsam_la_CPPFLAGS += -DCBLAS
#libgtsam_la_LDFLAGS += -L/opt/local/lib -lcblas -latlas -llapack -lptcblas -lptf77blas -lf77blas AM_LDFLAGS += $(BLAS_LIBS) $(LIBS) $(FLIBS) # -lcblas -latlas
libgtsam_la_LDFLAGS += $(BLAS_LIBS) $(LIBS) $(FLIBS) # -lcblas -latlas
endif
#######################
# GSL/ATLAS inclusion
#######################
#GSL using ATLAS
#if GSL
#AM_CXXFLAGS += -DGSL $(GSL_CFLAGS)
#libgtsam_la_CPPFLAGS += -DGSL $(GSL_CFLAGS)
#if ATLAS
#AM_LDFLAGS += $(GSL_LIBS_NO_CBLAS) -lcblas -latlas
#libgtsam_la_LDFLAGS += $(GSL_LIBS_NO_CBLAS) -lcblas -latlas
#else
#AM_LDFLAGS += $(GSL_LIBS)
#libgtsam_la_LDFLAGS += $(GSL_LIBS)
#endif
#endif
TESTS = $(check_PROGRAMS) TESTS = $(check_PROGRAMS)
CXXLINK = $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=link \ CXXLINK = $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=link \

63
cpp/Matrix.cpp
View File

@ -34,10 +34,6 @@
#include "Vector.h" #include "Vector.h"
#include "svdcmp.h" #include "svdcmp.h"
// use for switching quickly between GSL and nonGSL versions without reconfigure
//#define REVERTGSL
using namespace std; using namespace std;
namespace ublas = boost::numeric::ublas; namespace ublas = boost::numeric::ublas;
@ -158,7 +154,6 @@ bool assert_equal(const std::list<Matrix>& As, const std::list<Matrix>& Bs, doub
/* ************************************************************************* */ /* ************************************************************************* */
void multiplyAdd(double alpha, const Matrix& A, const Vector& x, Vector& e) { void multiplyAdd(double alpha, const Matrix& A, const Vector& x, Vector& e) {
#if defined CBLAS #if defined CBLAS
// get sizes // get sizes
const size_t m = A.size1(), n = A.size2(); const size_t m = A.size1(), n = A.size2();
@ -227,15 +222,15 @@ void transposeMultiplyAdd(double alpha, const Matrix& A, const Vector& e, Vector
// get pointers // get pointers
const double * Aptr = A.data().begin(); const double * Aptr = A.data().begin();
const double * Xptr = e.data().begin(); const double * Eptr = e.data().begin();
double * Eptr = x.data().begin(); double * Xptr = x.data().begin();
// fill in parameters // fill in parameters
const double beta = 1.0; const double beta = 1.0;
const int incx = 1, incy = 1, ida = n; const int incx = 1, incy = 1, ida = n;
// execute blas call // execute blas call
cblas_dgemv(CblasRowMajor, CblasTrans, m, n, alpha, Aptr, ida, Xptr, incx, beta, Eptr, incy); cblas_dgemv(CblasRowMajor, CblasTrans, m, n, alpha, Aptr, ida, Eptr, incx, beta, Xptr, incy);
#elif defined GSL #elif defined GSL
gsl_vector_const_view eg = gsl_vector_const_view_array(e.data().begin(), e.size()); gsl_vector_const_view eg = gsl_vector_const_view_array(e.data().begin(), e.size());
@ -482,28 +477,60 @@ inline void householder_update_manual(Matrix &A, int j, double beta, const Vecto
} }
void householder_update(Matrix &A, int j, double beta, const Vector& vjm) { void householder_update(Matrix &A, int j, double beta, const Vector& vjm) {
// TODO: SWAP IN ATLAS VERSION OF THE SYSTEM const size_t m = A.size1(), n = A.size2();
#if defined CBLAS
// CBLAS version not working, using manual approach
householder_update_manual(A,j,beta,vjm);
// // straight atlas version // // straight atlas version
// const size_t m = A.size1(), n = A.size2(), mj = m-j; // const size_t mj = m-j;
// //
// // find pointers to the data // // find pointers to the data
// const double * vptr = vjm.data().begin(); // mj long // const double * vptr = vjm.data().begin(); // mj long
// double * Aptr = A.data().begin() + n*j; // mj x n - note that this starts at row j // double * Aptr = A.data().begin() + n*j; // mj x n - note that this starts at row j
// //
// // first step: get w = beta*trans(A(j:m,:))*vjm // // first step: get w = beta*trans(A(j:m,:))*vjm
// Vector w(n); // Vector w = zero(n);
// double * wptr = w.data().begin(); // double * wptr = w.data().begin();
// //
// // DEBUG: create a duplicate version of the problem to solve simultaneously
// Matrix aA(A); Vector avjm(vjm);
//
// // execute w generation // // execute w generation
// cblas_dgemv(CblasRowMajor, CblasTrans, mj, n, beta, Aptr, n, vptr, 1, 0.0, wptr, 1); // cblas_dgemv(CblasRowMajor, CblasTrans, mj, n, beta, Aptr, n, vptr, 1, 0.0, wptr, 1);
// //
// // Execute w generation with alternate code
// Vector aw(n);
// for( size_t c = 0; c < n; c++) {
// aw(c) = 0.0;
// // dangerous as relies on row-major scheme
// const double *a = &aA(j,c), * const v = &avjm(0);
// for( size_t r=j, s=0 ; r < m ; r++, s++, a+=n )
// // w(c) += A(r,c) * vjm(r-j)
// aw(c) += (*a) * v[s];
// aw(c) *= beta;
// }
//
// print(w, "CBLAS w");
// print(aw, "Alternate w");
//
// // second step: rank 1 update A(j:m,:) = v(j:m)*w' + A(j:m,:) // // second step: rank 1 update A(j:m,:) = v(j:m)*w' + A(j:m,:)
// cblas_dger(CblasRowMajor, mj, n, 1.0, vptr, 1, wptr, 1, Aptr, n); // cblas_dger(CblasRowMajor, mj, n, 1.0, vptr, 1, wptr, 1, Aptr, n); // not correct
//
// // Execute second step using alternate code
// for( size_t c = 0 ; c < n; c++) {
// double wc = aw(c);
// double *a = &aA(j,c); const double * const v =&avjm(0);
// for( size_t r=j, s=0 ; r < m ; r++, s++, a+=n )
// // A(r,c) -= vjm(r-j) * wjn(c-j);
// (*a) -= v[s] * wc;
// }
//
// // copy in alternate results, which should be correct
// A = aA;
#elif defined GSL
#ifdef GSL
#ifndef REVERTGSL
const size_t m = A.size1(), n = A.size2();
// use GSL version // use GSL version
gsl_vector_const_view v = gsl_vector_const_view_array(vjm.data().begin(), m-j); gsl_vector_const_view v = gsl_vector_const_view_array(vjm.data().begin(), m-j);
gsl_matrix_view Ag = gsl_matrix_view_array(A.data().begin(), m, n); gsl_matrix_view Ag = gsl_matrix_view_array(A.data().begin(), m, n);
@ -512,10 +539,6 @@ void householder_update(Matrix &A, int j, double beta, const Vector& vjm) {
#else #else
householder_update_manual(A,j,beta,vjm); householder_update_manual(A,j,beta,vjm);
#endif #endif
#else
householder_update_manual(A,j,beta,vjm);
#endif
} }
/* ************************************************************************* */ /* ************************************************************************* */

12
cpp/NonlinearOptimizer-inl.h
View File

@ -149,11 +149,12 @@ namespace gtsam {
// update config // update config
shared_config newConfig(new C(expmap(*config_,delta))); // TODO: updateConfig shared_config newConfig(new C(expmap(*config_,delta))); // TODO: updateConfig
if (verbosity >= TRYCONFIG) // if (verbosity >= TRYCONFIG)
newConfig->print("config"); // newConfig->print("config");
// create new optimization state with more adventurous lambda // create new optimization state with more adventurous lambda
NonlinearOptimizer next(graph_, newConfig, solver_, lambda_ / factor); NonlinearOptimizer next(graph_, newConfig, solver_, lambda_ / factor);
cout << "next error = " << next.error_ << endl;
if(lambdaMode >= CAUTIOUS) { if(lambdaMode >= CAUTIOUS) {
throw runtime_error("CAUTIOUS mode not working yet, please use BOUNDED."); throw runtime_error("CAUTIOUS mode not working yet, please use BOUNDED.");
@ -171,8 +172,9 @@ namespace gtsam {
// Either we're not cautious, or we are but the adventerous lambda is better than the same one. // Either we're not cautious, or we are but the adventerous lambda is better than the same one.
return next; return next;
} else if (lambda_ > 1e+10) // if lambda gets too big, something is broken
} else { throw runtime_error("Lambda has grown too large!");
else {
// A more adventerous lambda was worse. If we're cautious, try the same lambda. // A more adventerous lambda was worse. If we're cautious, try the same lambda.
if(lambdaMode == CAUTIOUS) { if(lambdaMode == CAUTIOUS) {
@ -184,6 +186,7 @@ namespace gtsam {
// Either we're not cautious, or the same lambda was worse than the current error. // Either we're not cautious, or the same lambda was worse than the current error.
// The more adventerous lambda was worse too, so make lambda more conservative // The more adventerous lambda was worse too, so make lambda more conservative
// and keep the same config. // and keep the same config.
// TODO: can we avoid copying the config ? // TODO: can we avoid copying the config ?
if(lambdaMode >= BOUNDED && lambda_ >= 1.0e5) { if(lambdaMode >= BOUNDED && lambda_ >= 1.0e5) {
return NonlinearOptimizer(graph_, newConfig, solver_, lambda_);; return NonlinearOptimizer(graph_, newConfig, solver_, lambda_);;
@ -216,6 +219,7 @@ namespace gtsam {
linear->print("linear"); linear->print("linear");
// try lambda steps with successively larger lambda until we achieve descent // try lambda steps with successively larger lambda until we achieve descent
cout << "Trying Lambda for the first time" << endl;
return try_lambda(*linear, verbosity, lambdaFactor, lambdaMode); return try_lambda(*linear, verbosity, lambdaFactor, lambdaMode);
} }

2
myconfigure
View File

@ -1 +1 @@
./configure --prefix=$HOME --with-toolbox=$HOME/toolbox/ --with-boost=/opt/local/include/ --enable-gsl=no --enable-atlas=no CXXFLAGS=" -g -O2 -march=core2 -DNDEBUG" --disable-static ./configure --prefix=$HOME --with-toolbox=$HOME/toolbox/ --with-boost=/opt/local/include/ --enable-blas --with-blas=atlas CXXFLAGS=" -g -O2 -march=core2 -DNDEBUG" --disable-static