Change Cholesky decomposed matrix from 'Upper' to 'Lower' in order to avoid confusion.
parent
c5b4d731cc
commit
e9b0f7b98f
|
@ -95,7 +95,7 @@ void BlockJacobiPreconditioner::solve(const Vector& y, Vector &x) const {
|
||||||
|
|
||||||
const Eigen::Map<const Eigen::MatrixXd> R(ptr, d, d);
|
const Eigen::Map<const Eigen::MatrixXd> R(ptr, d, d);
|
||||||
Eigen::Map<Eigen::VectorXd> b(dst, d, 1);
|
Eigen::Map<Eigen::VectorXd> b(dst, d, 1);
|
||||||
R.triangularView<Eigen::Upper>().solveInPlace(b);
|
R.triangularView<Eigen::Lower>().solveInPlace(b);
|
||||||
|
|
||||||
dst += d;
|
dst += d;
|
||||||
ptr += sz;
|
ptr += sz;
|
||||||
|
@ -158,12 +158,12 @@ void BlockJacobiPreconditioner::build(
|
||||||
double *ptr = buffer_;
|
double *ptr = buffer_;
|
||||||
for ( size_t i = 0 ; i < n ; ++i ) {
|
for ( size_t i = 0 ; i < n ; ++i ) {
|
||||||
/* use eigen to decompose Di */
|
/* use eigen to decompose Di */
|
||||||
/* It is same as R = chol(M) in MATLAB where M is full preconditioner */
|
/* It is same as L = chol(M,'lower') in MATLAB where M is full preconditioner */
|
||||||
const Matrix R = blocks[i].llt().matrixL().transpose();
|
const Matrix L = blocks[i].llt().matrixL();
|
||||||
|
|
||||||
/* store the data in the buffer */
|
/* store the data in the buffer */
|
||||||
size_t sz = dims_[i]*dims_[i] ;
|
size_t sz = dims_[i]*dims_[i] ;
|
||||||
std::copy(R.data(), R.data() + sz, ptr);
|
std::copy(L.data(), L.data() + sz, ptr);
|
||||||
|
|
||||||
/* advance the pointer */
|
/* advance the pointer */
|
||||||
ptr += sz;
|
ptr += sz;
|
||||||
|
|
Loading…
Reference in New Issue