GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-08 17:59:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

SimplicialCholesky: the shift offset must be real, and fix a comparison issue for complexes

Browse Source
This commit is contained in:
Gael Guennebaud 2012-01-26 10:34:45 +01:00
parent d9f5840f7a
commit 65d5311c68
1 changed files with 8 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
Eigen/src/SparseCholesky/SimplicialCholesky.h
View File

@ -181,7 +181,7 @@ class SimplicialCholeskyBase
* *
* \returns a reference to \c *this. * \returns a reference to \c *this.
*/ */
Derived& setShift(const Scalar& offset, const RealScalar& scale = 1) Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
{ {
m_shiftOffset = offset; m_shiftOffset = offset;
m_shiftScale = scale; m_shiftScale = scale;
@ -281,7 +281,7 @@ class SimplicialCholeskyBase
PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation
PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation
Scalar m_shiftOffset; RealScalar m_shiftOffset;
RealScalar m_shiftScale; RealScalar m_shiftScale;
}; };
@ -739,7 +739,7 @@ void SimplicialCholeskyBase<Derived>::factorize(const MatrixType& a)
/* compute numerical values kth row of L (a sparse triangular solve) */ /* compute numerical values kth row of L (a sparse triangular solve) */
Scalar d = y[k] * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k) RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k)
y[k] = 0.0; y[k] = 0.0;
for(; top < size; ++top) for(; top < size; ++top)
{ {
@ -758,7 +758,7 @@ void SimplicialCholeskyBase<Derived>::factorize(const MatrixType& a)
Index p; Index p;
for(p = Lp[i] + (DoLDLt ? 0 : 1); p < p2; ++p) for(p = Lp[i] + (DoLDLt ? 0 : 1); p < p2; ++p)
y[Li[p]] -= internal::conj(Lx[p]) * yi; y[Li[p]] -= internal::conj(Lx[p]) * yi;
d -= l_ki * internal::conj(yi); d -= internal::real(l_ki * internal::conj(yi));
Li[p] = k; /* store L(k,i) in column form of L */ Li[p] = k; /* store L(k,i) in column form of L */
Lx[p] = l_ki; Lx[p] = l_ki;
++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
@ -766,7 +766,7 @@ void SimplicialCholeskyBase<Derived>::factorize(const MatrixType& a)
if(DoLDLt) if(DoLDLt)
{ {
m_diag[k] = d; m_diag[k] = d;
if(d == Scalar(0)) if(d == RealScalar(0))
{ {
ok = false; /* failure, D(k,k) is zero */ ok = false; /* failure, D(k,k) is zero */
break; break;
@ -774,9 +774,9 @@ void SimplicialCholeskyBase<Derived>::factorize(const MatrixType& a)
} }
else else
{ {
Index p = Lp[k]+m_nonZerosPerCol[k]++; Index p = Lp[k] + m_nonZerosPerCol[k]++;
Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */ Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */
if(d <= Scalar(0)) { if(d <= RealScalar(0)) {
ok = false; /* failure, matrix is not positive definite */ ok = false; /* failure, matrix is not positive definite */
break; break;
} }