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

the lapack version of householder

release/4.3a0
Kai Ni 2010-06-30 15:57:12 +00:00
parent 1b3a312abd
commit 70adcb26c9
6 changed files with 96 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
cpp/GaussianConditional.cpp
View File

@ -63,11 +63,13 @@ bool GaussianConditional::equals(const Conditional &c, double tol) const {
// check if the size of the parents_ map is the same
if (parents_.size() != p->parents_.size()) return false;
// check if R_ is equal
if (!(equal_with_abs_tol(R_, p->R_, tol))) return false;
bool equal_R1 = equal_with_abs_tol(R_, p->R_, tol);
bool equal_R2 = equal_with_abs_tol(R_*(-1), p->R_, tol);
bool equal_d1 = ::equal_with_abs_tol(d_, p->d_, tol);
bool equal_d2 = ::equal_with_abs_tol(d_*(-1), p->d_, tol);
// check if d_ is equal
if (!(::equal_with_abs_tol(d_, p->d_, tol))) return false;
// check if R_ and d_ are equal up to a sign
if (!((equal_R1 && equal_d1) || (equal_R2 && equal_d2))) return false;
// check if sigmas are equal
if (!(::equal_with_abs_tol(sigmas_, p->sigmas_, tol))) return false;
@ -77,7 +79,9 @@ bool GaussianConditional::equals(const Conditional &c, double tol) const {
for (it = parents_.begin(); it != parents_.end(); it++) {
Parents::const_iterator it2 = p->parents_.find(it->first);
if (it2 != p->parents_.end()) {
if (!(equal_with_abs_tol(it->second, it2->second, tol))) return false;
if (!(equal_with_abs_tol(it->second, it2->second, tol)) &&
!(equal_with_abs_tol(it->second*(-1), it2->second, tol)))
return false;
} else
return false;
}

37
cpp/Matrix.cpp
View File

@ -13,6 +13,7 @@
#ifdef GT_USE_CBLAS
#ifdef YA_BLAS
#include <vecLib/cblas.h>
#include <vecLib/clapack.h>
#else
#include <cblas.h>
#endif
@ -650,6 +651,42 @@ void householder(Matrix &A, size_t k) {
A(i,j) = 0.0;
}
/* ************************************************************************* */
/** in-place householder */
/* ************************************************************************* */
#ifdef GT_USE_CBLAS
void householder(Matrix &A) {
int m = A.size1();
int n = A.size2();
// convert from row major to column major
double a[m*n]; int k = 0;
for(int j=0; j<n; j++)
for(int i=0; i<m; i++, k++)
a[k] = A(i,j);
double tau[n];
double work_optimal_size;
int lwork = -1;
int info;
dgeqrf_(&m, &n, a, &m, tau, &work_optimal_size, &lwork, &info);
lwork = (int)work_optimal_size;
double work[lwork];
dgeqrf_(&m, &n, a, &m, tau, work, &lwork, &info);
int k0 = 0;
int j0;
memset(A.data().begin(), 0, m*n*sizeof(double));
for(int j=0; j<n; j++, k0+=m) {
k = k0;
j0 = min(j+1,m);
for(int i=0; i<j0; i++, k++)
A(i,j) = a[k];
}
}
#endif
/* ************************************************************************* */
Vector backSubstituteUpper(const Matrix& U, const Vector& b, bool unit) {
size_t n = U.size2();

9
cpp/Matrix.h
View File

@ -250,6 +250,15 @@ void householder_(Matrix& A, size_t k);
*/
void householder(Matrix& A, size_t k);
/**
* Householder tranformation, zeros below diagonal
* @param k number of columns to zero out below diagonal
* @return nothing: in place !!!
*/
#ifdef GT_USE_CBLAS
void householder(Matrix &A);
#endif
/**
* backSubstitute U*x=b
* @param U an upper triangular matrix

5
cpp/NoiseModel.cpp
View File

@ -118,7 +118,12 @@ SharedDiagonal Gaussian::QR(Matrix& Ab) const {
if (verbose) gtsam::print(Ab, "Ab after whitening");
// Perform in-place Householder
#ifdef GT_USE_CBLAS
householder(Ab);
#else
householder(Ab, maxRank);
#endif
if (verbose) gtsam::print(Ab, "Ab before householder");
return Unit::Create(maxRank);

5
cpp/testGaussianFactorGraph.cpp
View File

@ -729,7 +729,10 @@ TEST( GaussianFactorGraph, elimination )
Matrix expected = Matrix_(2,2,
0.707107, -0.353553,
0.0, 0.612372);
CHECK(assert_equal(expected,R,1e-6));
Matrix expected2 = Matrix_(2,2,
0.707107, -0.353553,
0.0, -0.612372);
CHECK(equal_with_abs_tol(expected, R, 1e-6) || equal_with_abs_tol(expected2, R, 1e-6));
}
/* ************************************************************************* */

37
cpp/testMatrix.cpp
View File

@ -613,13 +613,16 @@ TEST( matrix, backsubtitution )
/* ************************************************************************* */
TEST( matrix, houseHolder )
{
double data[] = { -5, 0, 5, 0, 0, 0, -1, 00, -5, 0, 5, 0, 0, 1.5, 10, 0, 0,
0, -10, 0, 2, 00, 10, 0, 0, 0, -10, -1 };
double data[] = { -5, 0, 5, 0, 0, 0, -1,
00,-5, 0, 5, 0, 0, 1.5,
10, 0, 0, 0,-10,0, 2,
00, 10,0, 0, 0, -10, -1 };
// check in-place householder, with v vectors below diagonal
double data1[] = { 11.1803, 0, -2.2361, 0, -8.9443, 0, 2.236, 0, 11.1803,
0, -2.2361, 0, -8.9443, -1.565, -0.618034, 0, 4.4721, 0, -4.4721,
0, 0, 0, -0.618034, 0, 4.4721, 0, -4.4721, 0.894 };
double data1[] = { 11.1803, 0, -2.2361, 0, -8.9443, 0, 2.236,
0, 11.1803, 0, -2.2361, 0, -8.9443, -1.565,
-0.618034, 0, 4.4721, 0, -4.4721, 0, 0,
0, -0.618034, 0, 4.4721, 0, -4.4721, 0.894 };
Matrix expected1 = Matrix_(4, 7, data1);
Matrix A1 = Matrix_(4, 7, data);
householder_(A1, 3);
@ -634,6 +637,30 @@ TEST( matrix, houseHolder )
householder(A2, 3);
CHECK(assert_equal(expected, A2, 1e-3));
}
/* ************************************************************************* */
// unit tests for housholder transformation
/* ************************************************************************* */
#ifdef GT_USE_CBLAS
TEST( matrix, houseHolder2 )
{
double data[] = { -5, 0, 5, 0, 0, 0, -1,
00,-5, 0, 5, 0, 0, 1.5,
10, 0, 0, 0,-10,0, 2,
00, 10,0, 0, 0, -10, -1 };
// check in-place householder, with v vectors below diagonal
double data1[] = { 11.1803, 0, -2.2361, 0, -8.9443, 0, 2.236,
0, 11.1803, 0, -2.2361, 0, -8.9443, -1.565,
0, 0, 4.4721, 0, -4.4721, 0, 0,
0, 0, 0, 4.4721, 0, -4.4721, 0.894 };
Matrix expected1 = Matrix_(4, 7, data1);
Matrix A1 = Matrix_(4, 7, data);
householder(A1);
CHECK(assert_equal(expected1, A1, 1e-3));
}
#endif
/* ************************************************************************* */
// unit test for qr factorization (and hence householder)
// This behaves the same as QR in matlab: [Q,R] = qr(A), except for signs