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

1.0 versions of BLAS level 2

release/4.3a0
Frank Dellaert 2010-02-20 23:44:07 +00:00
parent ac316f2d56
commit 5b6af9e4bd
2 changed files with 37 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
cpp/Matrix.cpp
View File

@ -156,6 +156,19 @@ void multiplyAdd(double alpha, const Matrix& A, const Vector& x, Vector& e) {
#endif
}
/* ************************************************************************* */
void multiplyAdd(const Matrix& A, const Vector& x, Vector& 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);
}
}
/* ************************************************************************* */
Vector operator^(const Matrix& A, const Vector & v) {
if (A.size1()!=v.size()) throw std::invalid_argument(
@ -186,6 +199,20 @@ void transposeMultiplyAdd(double alpha, const Matrix& A, const Vector& e, Vector
#endif
}
/* ************************************************************************* */
void transposeMultiplyAdd(const Matrix& A, const Vector& e, Vector& 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);
}
}
/* ************************************************************************* */
void transposeMultiplyAdd(double alpha, const Matrix& A, const Vector& e, SubVector x) {
// ublas x += prod(trans(A),e) is terribly slow

10
cpp/Matrix.h
View File

@ -91,6 +91,11 @@ inline Vector operator*(const Matrix& A, const Vector & v) { return prod(A,v);}
*/
void multiplyAdd(double alpha, const Matrix& A, const Vector& x, Vector& e);
/**
* BLAS Level-2 style e <- e + A*x
*/
void multiplyAdd(const Matrix& A, const Vector& x, Vector& e);
/**
* overload ^ for trans(A)*v
* We transpose the vectors for speed.
@ -102,6 +107,11 @@ Vector operator^(const Matrix& A, const Vector & v);
*/
void transposeMultiplyAdd(double alpha, const Matrix& A, const Vector& e, Vector& x);
/**
* BLAS Level-2 style x <- x + A'*e
*/
void transposeMultiplyAdd(const Matrix& A, const Vector& e, Vector& x);
/**
* BLAS Level-2 style x <- x + alpha*A'*e
*/