aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h178
-rw-r--r--test/incomplete_cholesky.cpp30
2 files changed, 131 insertions, 77 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h
index 18e9d1466..babc14d9e 100644
--- a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h
+++ b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h
@@ -240,7 +240,7 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
else
m_scale(j) = 1;
- // FIXME disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
+ // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
// Scale and compute the shift for the matrix
RealScalar mindiag = NumTraits<RealScalar>::highest();
@@ -251,96 +251,122 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
}
+
+ FactorType L_save = m_L;
RealScalar shift = 0;
if(mindiag <= RealScalar(0.))
shift = m_initialShift - mindiag;
- // Apply the shift to the diagonal elements of the matrix
- for (Index j = 0; j < n; j++)
- vals[colPtr[j]] += shift;
-
- // jki version of the Cholesky factorization
- for (Index j=0; j < n; ++j)
- {
- // Left-looking factorization of the j-th column
- // First, load the j-th column into col_vals
- Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
- col_nnz = 0;
- for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
- {
- StorageIndex l = rowIdx[i];
- col_vals(col_nnz) = vals[i];
- col_irow(col_nnz) = l;
- col_pattern(l) = col_nnz;
- col_nnz++;
- }
+ m_info = NumericalIssue;
+
+ // Try to perform the incomplete factorization using the current shift
+ int iter = 0;
+ do
+ {
+ // Apply the shift to the diagonal elements of the matrix
+ for (Index j = 0; j < n; j++)
+ vals[colPtr[j]] += shift;
+
+ // jki version of the Cholesky factorization
+ Index j=0;
+ for (; j < n; ++j)
{
- typename std::list<StorageIndex>::iterator k;
- // Browse all previous columns that will update column j
- for(k = listCol[j].begin(); k != listCol[j].end(); k++)
+ // Left-looking factorization of the j-th column
+ // First, load the j-th column into col_vals
+ Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
+ col_nnz = 0;
+ for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
{
- Index jk = firstElt(*k); // First element to use in the column
- eigen_internal_assert(rowIdx[jk]==j);
- Scalar v_j_jk = numext::conj(vals[jk]);
-
- jk += 1;
- for (Index i = jk; i < colPtr[*k+1]; i++)
+ StorageIndex l = rowIdx[i];
+ col_vals(col_nnz) = vals[i];
+ col_irow(col_nnz) = l;
+ col_pattern(l) = col_nnz;
+ col_nnz++;
+ }
+ {
+ typename std::list<StorageIndex>::iterator k;
+ // Browse all previous columns that will update column j
+ for(k = listCol[j].begin(); k != listCol[j].end(); k++)
{
- StorageIndex l = rowIdx[i];
- if(col_pattern[l]<0)
+ Index jk = firstElt(*k); // First element to use in the column
+ eigen_internal_assert(rowIdx[jk]==j);
+ Scalar v_j_jk = numext::conj(vals[jk]);
+
+ jk += 1;
+ for (Index i = jk; i < colPtr[*k+1]; i++)
{
- col_vals(col_nnz) = vals[i] * v_j_jk;
- col_irow[col_nnz] = l;
- col_pattern(l) = col_nnz;
- col_nnz++;
+ StorageIndex l = rowIdx[i];
+ if(col_pattern[l]<0)
+ {
+ col_vals(col_nnz) = vals[i] * v_j_jk;
+ col_irow[col_nnz] = l;
+ col_pattern(l) = col_nnz;
+ col_nnz++;
+ }
+ else
+ col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
}
- else
- col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
+ updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
}
- updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
}
+
+ // Scale the current column
+ if(numext::real(diag) <= 0)
+ {
+ if(++iter>=10)
+ return;
+
+ // increase shift
+ shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
+ // restore m_L, col_pattern, and listCol
+ vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
+ rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
+ colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
+ col_pattern.fill(-1);
+ for(Index i=0; i<n; ++i)
+ listCol[i].clear();
+
+ break;
+ }
+
+ RealScalar rdiag = sqrt(numext::real(diag));
+ vals[colPtr[j]] = rdiag;
+ for (Index k = 0; k<col_nnz; ++k)
+ {
+ Index i = col_irow[k];
+ //Scale
+ col_vals(k) /= rdiag;
+ //Update the remaining diagonals with col_vals
+ vals[colPtr[i]] -= numext::abs2(col_vals(k));
+ }
+ // Select the largest p elements
+ // p is the original number of elements in the column (without the diagonal)
+ Index p = colPtr[j+1] - colPtr[j] - 1 ;
+ Ref<VectorSx> cvals = col_vals.head(col_nnz);
+ Ref<VectorIx> cirow = col_irow.head(col_nnz);
+ internal::QuickSplit(cvals,cirow, p);
+ // Insert the largest p elements in the matrix
+ Index cpt = 0;
+ for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
+ {
+ vals[i] = col_vals(cpt);
+ rowIdx[i] = col_irow(cpt);
+ // restore col_pattern:
+ col_pattern(col_irow(cpt)) = -1;
+ cpt++;
+ }
+ // Get the first smallest row index and put it after the diagonal element
+ Index jk = colPtr(j)+1;
+ updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
}
-
- // Scale the current column
- if(numext::real(diag) <= 0)
- {
- m_info = NumericalIssue;
- return;
- }
-
- RealScalar rdiag = sqrt(numext::real(diag));
- vals[colPtr[j]] = rdiag;
- for (Index k = 0; k<col_nnz; ++k)
- {
- Index i = col_irow[k];
- //Scale
- col_vals(k) /= rdiag;
- //Update the remaining diagonals with col_vals
- vals[colPtr[i]] -= numext::abs2(col_vals(k));
- }
- // Select the largest p elements
- // p is the original number of elements in the column (without the diagonal)
- Index p = colPtr[j+1] - colPtr[j] - 1 ;
- Ref<VectorSx> cvals = col_vals.head(col_nnz);
- Ref<VectorIx> cirow = col_irow.head(col_nnz);
- internal::QuickSplit(cvals,cirow, p);
- // Insert the largest p elements in the matrix
- Index cpt = 0;
- for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
+
+ if(j==n)
{
- vals[i] = col_vals(cpt);
- rowIdx[i] = col_irow(cpt);
- // restore col_pattern:
- col_pattern(col_irow(cpt)) = -1;
- cpt++;
+ m_factorizationIsOk = true;
+ m_info = Success;
}
- // Get the first smallest row index and put it after the diagonal element
- Index jk = colPtr(j)+1;
- updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
- }
- m_factorizationIsOk = true;
- m_info = Success;
+ } while(m_info!=Success);
}
template<typename Scalar, int _UpLo, typename OrderingType>
diff --git a/test/incomplete_cholesky.cpp b/test/incomplete_cholesky.cpp
index 7acad9872..59ffe9259 100644
--- a/test/incomplete_cholesky.cpp
+++ b/test/incomplete_cholesky.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2015-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -34,4 +34,32 @@ void test_incomplete_cholesky()
CALL_SUBTEST_1(( test_incomplete_cholesky_T<double,int>() ));
CALL_SUBTEST_2(( test_incomplete_cholesky_T<std::complex<double>, int>() ));
CALL_SUBTEST_3(( test_incomplete_cholesky_T<double,long int>() ));
+
+#ifdef EIGEN_TEST_PART_1
+ // regression for bug 1150
+ for(int N = 1; N<20; ++N)
+ {
+ Eigen::MatrixXd b( N, N );
+ b.setOnes();
+
+ Eigen::SparseMatrix<double> m( N, N );
+ m.reserve(Eigen::VectorXi::Constant(N,4));
+ for( int i = 0; i < N; ++i )
+ {
+ m.insert( i, i ) = 1;
+ m.coeffRef( i, i / 2 ) = 2;
+ m.coeffRef( i, i / 3 ) = 2;
+ m.coeffRef( i, i / 4 ) = 2;
+ }
+
+ Eigen::SparseMatrix<double> A;
+ A = m * m.transpose();
+
+ Eigen::ConjugateGradient<Eigen::SparseMatrix<double>,
+ Eigen::Lower | Eigen::Upper,
+ Eigen::IncompleteCholesky<double> > solver( A );
+ VERIFY(solver.preconditioner().info() == Eigen::Success);
+ VERIFY(solver.info() == Eigen::Success);
+ }
+#endif
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-05 01:22:10 +0000