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

Added pure cblas implementation of multiplyAdd functions. This does not have autotools support yet, so to enable, goto cpp/Makefile.am, comment out the block concerning GSL/ATLAS, and uncomment the above section enabling just blas.

release/4.3a0
Alex Cunningham 2010-02-24 18:12:48 +00:00
parent f8c4e1fe27
commit 9955ea20bd
2 changed files with 65 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
cpp/Makefile.am
View File

@ -293,6 +293,12 @@ include_HEADERS = $(headers)
AM_CXXFLAGS += -I..
AM_LDFLAGS = -L../CppUnitLite -lCppUnitLite $(BOOST_LDFLAGS) $(boost_serialization) #-L. -lgtsam
# TESTING: adding cblas implementation with atlas
#AM_CXXFLAGS += -DCBLAS
#libgtsam_la_CPPFLAGS += -DCBLAS
#AM_LDFLAGS += -lcblas -latlas
#libgtsam_la_LDFLAGS += -lcblas -latlas
#######################
# GSL/ATLAS inclusion
#######################

77
cpp/Matrix.cpp
View File

@ -10,6 +10,11 @@
#include <list>
#include <fstream>
// TESTING
#ifdef CBLAS
#include <cblas.h>
#endif
#ifdef GSL
#include <gsl/gsl_blas.h> // needed for gsl blas
#include <gsl/gsl_linalg.h>
@ -139,7 +144,24 @@ bool assert_equal(const Matrix& expected, const Matrix& actual, double tol) {
/* ************************************************************************* */
void multiplyAdd(double alpha, const Matrix& A, const Vector& x, Vector& e) {
#ifdef GSL
#if defined CBLAS
// get sizes
const size_t m = A.size1(), n = A.size2();
// get pointers
const double * Aptr = A.data().begin();
const double * Xptr = x.data().begin();
double * Eptr = e.data().begin();
// fill in parameters
const double beta = 1.0;
const int incx = 1, incy = 1, ida = n;
// execute blas call
cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, Aptr, ida, Xptr, incx, beta, Eptr, incy);
#elif defined GSL
gsl_vector_const_view xg = gsl_vector_const_view_array(x.data().begin(), x.size());
gsl_vector_view eg = gsl_vector_view_array(e.data().begin(), e.size());
gsl_matrix_const_view Ag = gsl_matrix_const_view_array(A.data().begin(), A.size1(), A.size2());
@ -159,15 +181,16 @@ void multiplyAdd(double alpha, const Matrix& A, const Vector& x, Vector& e) {
/* ************************************************************************* */
void multiplyAdd(const Matrix& A, const Vector& x, Vector& e) {
multiplyAdd(1.0, A, x, e);
// ublas e += prod(A,x) is terribly slow
size_t m = A.size1(), n = A.size2();
double * ei = e.data().begin();
const double * aij = A.data().begin();
for (int i = 0; i < m; i++, ei++) {
const double * xj = x.data().begin();
for (int j = 0; j < n; j++, aij++, xj++)
(*ei) += (*aij) * (*xj);
}
// size_t m = A.size1(), n = A.size2();
// double * ei = e.data().begin();
// const double * aij = A.data().begin();
// for (int i = 0; i < m; i++, ei++) {
// const double * xj = x.data().begin();
// for (int j = 0; j < n; j++, aij++, xj++)
// (*ei) += (*aij) * (*xj);
// }
}
/* ************************************************************************* */
@ -181,7 +204,24 @@ Vector operator^(const Matrix& A, const Vector & v) {
/* ************************************************************************* */
void transposeMultiplyAdd(double alpha, const Matrix& A, const Vector& e, Vector& x) {
#ifdef GSL
#if defined CBLAS
// get sizes
const size_t m = A.size1(), n = A.size2();
// get pointers
const double * Aptr = A.data().begin();
const double * Xptr = e.data().begin();
double * Eptr = x.data().begin();
// fill in parameters
const double beta = 1.0;
const int incx = 1, incy = 1, ida = n;
// execute blas call
cblas_dgemv(CblasRowMajor, CblasTrans, m, n, alpha, Aptr, ida, Xptr, incx, beta, Eptr, incy);
#elif defined GSL
gsl_vector_const_view eg = gsl_vector_const_view_array(e.data().begin(), e.size());
gsl_vector_view xg = gsl_vector_view_array(x.data().begin(), x.size());
gsl_matrix_const_view Ag = gsl_matrix_const_view_array(A.data().begin(), A.size1(), A.size2());
@ -202,16 +242,17 @@ void transposeMultiplyAdd(double alpha, const Matrix& A, const Vector& e, Vector
/* ************************************************************************* */
void transposeMultiplyAdd(const Matrix& A, const Vector& e, Vector& x) {
transposeMultiplyAdd(1.0, A, e, x);
// ublas x += prod(trans(A),e) is terribly slow
// TODO: use BLAS
size_t m = A.size1(), n = A.size2();
double * xj = x.data().begin();
for (int j = 0; j < n; j++,xj++) {
const double * ei = e.data().begin();
const double * aij = A.data().begin() + j;
for (int i = 0; i < m; i++, aij+=n, ei++)
(*xj) += (*aij) * (*ei);
}
// size_t m = A.size1(), n = A.size2();
// double * xj = x.data().begin();
// for (int j = 0; j < n; j++,xj++) {
// const double * ei = e.data().begin();
// const double * aij = A.data().begin() + j;
// for (int i = 0; i < m; i++, aij+=n, ei++)
// (*xj) += (*aij) * (*ei);
// }
}
/* ************************************************************************* */