Split the computation of the ILUT into two steps
diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
index c0a9b12..d813ea8 100644
--- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
+++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
@@ -35,9 +35,10 @@
* approximation of Ax=b (regardless of b)
* \param iters On input the max number of iteration, on output the number of performed iterations.
* \param tol_error On input the tolerance error, on output an estimation of the relative error.
+ * \return false in the case of numerical issue, for example a break down of BiCGSTAB.
*/
template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
-void bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
+bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
const Preconditioner& precond, int& iters,
typename Dest::RealScalar& tol_error)
{
@@ -46,7 +47,6 @@
typedef typename Dest::RealScalar RealScalar;
typedef typename Dest::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,1> VectorType;
-
RealScalar tol = tol_error;
int maxIters = iters;
@@ -70,10 +70,11 @@
while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters )
{
+// std::cout<<i<<" : Relative residual norm " << sqrt(r.squaredNorm()/r0_sqnorm)<<"\n";
Scalar rho_old = rho;
rho = r0.dot(r);
- eigen_assert((rho != Scalar(0)) && "BiCGSTAB BROKE DOWN !!!");
+ if (rho == Scalar(0)) return false; /* New search directions cannot be found */
Scalar beta = (rho/rho_old) * (alpha / w);
p = r + beta * (p - w * v);
@@ -94,6 +95,7 @@
}
tol_error = sqrt(r.squaredNorm()/r0_sqnorm);
iters = i;
+ return true;
}
}
@@ -214,17 +216,18 @@
template<typename Rhs,typename Dest>
void _solveWithGuess(const Rhs& b, Dest& x) const
{
+ bool ok;
for(int j=0; j<b.cols(); ++j)
{
m_iterations = Base::maxIterations();
m_error = Base::m_tolerance;
typename Dest::ColXpr xj(x,j);
- internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
+ ok = internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
}
-
+ if (ok == false) m_info = NumericalIssue;
+ else m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
m_isInitialized = true;
- m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
}
/** \internal */
diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
index bfa6e29..ce451ae 100644
--- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
+++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
@@ -63,220 +63,237 @@
public:
typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
- IncompleteLUT() : m_droptol(NumTraits<Scalar>::dummy_precision()),m_fillfactor(50),m_isInitialized(false) {};
+ IncompleteLUT() : m_droptol(NumTraits<Scalar>::dummy_precision()),m_fillfactor(10),m_isInitialized(false),m_analysisIsOk(false),m_factorizationIsOk(false) {};
template<typename MatrixType>
IncompleteLUT(const MatrixType& mat, RealScalar droptol, int fillfactor)
- : m_droptol(droptol),m_fillfactor(fillfactor),m_isInitialized(false)
+ : m_droptol(droptol),m_fillfactor(fillfactor),m_isInitialized(false),m_analysisIsOk(false),m_factorizationIsOk(false)
{
- compute(mat);
+ eigen_assert(fillfactor != 0);
+ compute(mat);
}
Index rows() const { return m_lu.rows(); }
Index cols() const { return m_lu.cols(); }
-
-/**
- * Compute an incomplete LU factorization with dual threshold on the matrix mat
- * No pivoting is done in this version
- *
- **/
-template<typename MatrixType>
-IncompleteLUT<Scalar>& compute(const MatrixType& amat)
-{
- int n = amat.cols(); /* Size of the matrix */
- m_lu.resize(n,n);
- int fill_in; /* Number of largest elements to keep in each row */
- int nnzL, nnzU; /* Number of largest nonzero elements to keep in the L and the U part of the current row */
- /* Declare Working vectors and variables */
- int sizeu; /* number of nonzero elements in the upper part of the current row */
- int sizel; /* number of nonzero elements in the lower part of the current row */
- Vector u(n) ; /* real values of the row -- maximum size is n -- */
- VectorXi ju(n); /*column position of the values in u -- maximum size is n*/
- VectorXi jr(n); /* Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1*/
- int j, k, ii, jj, jpos, minrow, len;
- Scalar fact, prod;
- RealScalar rownorm;
-
- /* Compute the Fill-reducing permutation */
- SparseMatrix<Scalar,ColMajor, Index> mat1 = amat;
- SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose();
- SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 * mat1; /* Symmetrize the pattern */
- AtA.prune(keep_diag());
- internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P); /* Then compute the AMD ordering... */
-
- m_Pinv = m_P.inverse(); /* ... and the inverse permutation */
-
-// m_Pinv.indices().setLinSpaced(0,n);
-// m_P.indices().setLinSpaced(0,n);
-
- SparseMatrix<Scalar,RowMajor, Index> mat;
- mat = amat.twistedBy(m_Pinv);
- /* Initialization */
- fact = 0;
- jr.fill(-1);
- ju.fill(0);
- u.fill(0);
- fill_in = static_cast<int> (amat.nonZeros()*m_fillfactor)/n+1;
- if (fill_in > n) fill_in = n;
- nnzL = fill_in/2;
- nnzU = nnzL;
- m_lu.reserve(n * (nnzL + nnzU + 1));
- for (int ii = 0; ii < n; ii++)
- { /* global loop over the rows of the sparse matrix */
-
- /* Copy the lower and the upper part of the row i of mat in the working vector u */
- sizeu = 1;
- sizel = 0;
- ju(ii) = ii;
- u(ii) = 0;
- jr(ii) = ii;
- rownorm = 0;
-
- typename FactorType::InnerIterator j_it(mat, ii); /* Iterate through the current row ii */
- for (; j_it; ++j_it)
+ template<typename MatrixType>
+ void analyzePattern(const MatrixType& amat)
{
- k = j_it.index();
- if (k < ii)
- { /* Copy the lower part */
- ju(sizel) = k;
- u(sizel) = j_it.value();
- jr(k) = sizel;
- ++sizel;
- }
- else if (k == ii)
- {
- u(ii) = j_it.value();
- }
- else
- { /* Copy the upper part */
- jpos = ii + sizeu;
- ju(jpos) = k;
- u(jpos) = j_it.value();
- jr(k) = jpos;
- ++sizeu;
- }
- rownorm += internal::abs2(j_it.value());
- } /* end copy of the row */
- /* detect possible zero row */
- if (rownorm == 0) eigen_internal_assert(false);
- rownorm = std::sqrt(rownorm); /* Take the 2-norm of the current row as a relative tolerance */
-
- /* Now, eliminate the previous nonzero rows */
- jj = 0; len = 0;
- while (jj < sizel)
- { /* In order to eliminate in the correct order, we must select first the smallest column index among ju(jj:sizel) */
-
- minrow = ju.segment(jj,sizel-jj).minCoeff(&k); /* k est relatif au segment */
- k += jj;
- if (minrow != ju(jj)) { /* swap the two locations */
- j = ju(jj);
- std::swap(ju(jj), ju(k));
- jr(minrow) = jj; jr(j) = k;
- std::swap(u(jj), u(k));
- }
- /* Reset this location to zero */
- jr(minrow) = -1;
+ /* Compute the Fill-reducing permutation */
+ SparseMatrix<Scalar,ColMajor, Index> mat1 = amat;
+ SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose();
+ SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 * mat1; /* Symmetrize the pattern */
+ AtA.prune(keep_diag());
+ internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P); /* Then compute the AMD ordering... */
- /* Start elimination */
- typename FactorType::InnerIterator ki_it(m_lu, minrow);
- while (ki_it && ki_it.index() < minrow) ++ki_it;
- if(ki_it && ki_it.col()==minrow) fact = u(jj) / ki_it.value();
- else { eigen_internal_assert(false); }
- if( std::abs(fact) <= m_droptol ) {
- jj++;
- continue ; /* This element is been dropped */
- }
- /* linear combination of the current row ii and the row minrow */
- ++ki_it;
- for (; ki_it; ++ki_it) {
- prod = fact * ki_it.value();
- j = ki_it.index();
- jpos = jr(j);
- if (j >= ii) { /* Dealing with the upper part */
- if (jpos == -1) { /* Fill-in element */
- int newpos = ii + sizeu;
- ju(newpos) = j;
- u(newpos) = - prod;
- jr(j) = newpos;
- sizeu++;
- if (sizeu > n) { eigen_internal_assert(false);}
+ m_Pinv = m_P.inverse(); /* ... and the inverse permutation */
+ m_analysisIsOk = true;
+ }
+
+ template<typename MatrixType>
+ void factorize(const MatrixType& amat)
+ {
+ eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix");
+ int n = amat.cols(); /* Size of the matrix */
+ m_lu.resize(n,n);
+ int fill_in; /* Number of largest elements to keep in each row */
+ int nnzL, nnzU; /* Number of largest nonzero elements to keep in the L and the U part of the current row */
+ /* Declare Working vectors and variables */
+ int sizeu; /* number of nonzero elements in the upper part of the current row */
+ int sizel; /* number of nonzero elements in the lower part of the current row */
+ Vector u(n) ; /* real values of the row -- maximum size is n -- */
+ VectorXi ju(n); /*column position of the values in u -- maximum size is n*/
+ VectorXi jr(n); /* Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1*/
+ int j, k, ii, jj, jpos, minrow, len;
+ Scalar fact, prod;
+ RealScalar rownorm;
+
+ /* Apply the fill-reducing permutation */
+
+ eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
+ SparseMatrix<Scalar,RowMajor, Index> mat;
+ mat = amat.twistedBy(m_Pinv);
+
+ /* Initialization */
+ fact = 0;
+ jr.fill(-1);
+ ju.fill(0);
+ u.fill(0);
+ fill_in = static_cast<int> (amat.nonZeros()*m_fillfactor)/n+1;
+ if (fill_in > n) fill_in = n;
+ nnzL = fill_in/2;
+ nnzU = nnzL;
+ m_lu.reserve(n * (nnzL + nnzU + 1));
+ for (int ii = 0; ii < n; ii++)
+ { /* global loop over the rows of the sparse matrix */
+
+ /* Copy the lower and the upper part of the row i of mat in the working vector u */
+ sizeu = 1;
+ sizel = 0;
+ ju(ii) = ii;
+ u(ii) = 0;
+ jr(ii) = ii;
+ rownorm = 0;
+
+ typename FactorType::InnerIterator j_it(mat, ii); /* Iterate through the current row ii */
+ for (; j_it; ++j_it)
+ {
+ k = j_it.index();
+ if (k < ii)
+ { /* Copy the lower part */
+ ju(sizel) = k;
+ u(sizel) = j_it.value();
+ jr(k) = sizel;
+ ++sizel;
}
- else { /* Not a fill_in element */
- u(jpos) -= prod;
+ else if (k == ii)
+ {
+ u(ii) = j_it.value();
+ }
+ else
+ { /* Copy the upper part */
+ jpos = ii + sizeu;
+ ju(jpos) = k;
+ u(jpos) = j_it.value();
+ jr(k) = jpos;
+ ++sizeu;
+ }
+ rownorm += internal::abs2(j_it.value());
+ } /* end copy of the row */
+ /* detect possible zero row */
+ if (rownorm == 0) eigen_internal_assert(false);
+ rownorm = std::sqrt(rownorm); /* Take the 2-norm of the current row as a relative tolerance */
+
+ /* Now, eliminate the previous nonzero rows */
+ jj = 0; len = 0;
+ while (jj < sizel)
+ { /* In order to eliminate in the correct order, we must select first the smallest column index among ju(jj:sizel) */
+
+ minrow = ju.segment(jj,sizel-jj).minCoeff(&k); /* k est relatif au segment */
+ k += jj;
+ if (minrow != ju(jj)) { /* swap the two locations */
+ j = ju(jj);
+ std::swap(ju(jj), ju(k));
+ jr(minrow) = jj; jr(j) = k;
+ std::swap(u(jj), u(k));
+ }
+ /* Reset this location to zero */
+ jr(minrow) = -1;
+
+ /* Start elimination */
+ typename FactorType::InnerIterator ki_it(m_lu, minrow);
+ while (ki_it && ki_it.index() < minrow) ++ki_it;
+ if(ki_it && ki_it.col()==minrow) fact = u(jj) / ki_it.value();
+ else { eigen_internal_assert(false); }
+ if( std::abs(fact) <= m_droptol ) {
+ jj++;
+ continue ; /* This element is been dropped */
+ }
+ /* linear combination of the current row ii and the row minrow */
+ ++ki_it;
+ for (; ki_it; ++ki_it) {
+ prod = fact * ki_it.value();
+ j = ki_it.index();
+ jpos = jr(j);
+ if (j >= ii) { /* Dealing with the upper part */
+ if (jpos == -1) { /* Fill-in element */
+ int newpos = ii + sizeu;
+ ju(newpos) = j;
+ u(newpos) = - prod;
+ jr(j) = newpos;
+ sizeu++;
+ if (sizeu > n) { eigen_internal_assert(false);}
+ }
+ else { /* Not a fill_in element */
+ u(jpos) -= prod;
+ }
+ }
+ else { /* Dealing with the lower part */
+ if (jpos == -1) { /* Fill-in element */
+ ju(sizel) = j;
+ jr(j) = sizel;
+ u(sizel) = - prod;
+ sizel++;
+ if(sizel > n) { eigen_internal_assert(false);}
+ }
+ else {
+ u(jpos) -= prod;
+ }
+ }
+ }
+ /* Store the pivot element */
+ u(len) = fact;
+ ju(len) = minrow;
+ ++len;
+
+ jj++;
+ } /* End While loop -- end of the elimination on the row ii*/
+ /* Reset the upper part of the pointer jr to zero */
+ for (k = 0; k <sizeu; k++){
+ jr(ju(ii+k)) = -1;
+ }
+ /* Sort the L-part of the row --use Quicksplit()*/
+ sizel = len;
+ len = std::min(sizel, nnzL );
+ typename Vector::SegmentReturnType ul(u.segment(0, len));
+ typename VectorXi::SegmentReturnType jul(ju.segment(0, len));
+ QuickSplit(ul, jul, len);
+
+
+ /* Store the largest m_fill elements of the L part */
+ m_lu.startVec(ii);
+ for (k = 0; k < len; k++){
+ m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
+ }
+
+ /* Store the diagonal element */
+ if (u(ii) == Scalar(0))
+ u(ii) = std::sqrt(m_droptol ) * rownorm ; /* NOTE This is used to avoid a zero pivot, because we are doing an incomplete factorization */
+ m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
+ /* Sort the U-part of the row -- Use Quicksplit() */
+ len = 0;
+ for (k = 1; k < sizeu; k++) { /* First, drop any element that is below a relative tolerance */
+ if ( std::abs(u(ii+k)) > m_droptol * rownorm ) {
+ ++len;
+ u(ii + len) = u(ii + k);
+ ju(ii + len) = ju(ii + k);
}
}
- else { /* Dealing with the lower part */
- if (jpos == -1) { /* Fill-in element */
- ju(sizel) = j;
- jr(j) = sizel;
- u(sizel) = - prod;
- sizel++;
- if(sizel > n) { eigen_internal_assert(false);}
- }
- else {
- u(jpos) -= prod;
- }
+ sizeu = len + 1; /* To take into account the diagonal element */
+ len = std::min(sizeu, nnzU);
+ typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
+ typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
+ QuickSplit(uu, juu, len);
+ /* Store the largest <fill> elements of the U part */
+ for (k = ii + 1; k < ii + len; k++){
+ m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
}
- }
- /* Store the pivot element */
- u(len) = fact;
- ju(len) = minrow;
- ++len;
+ } /* End global for-loop */
+ m_lu.finalize();
+ m_lu.makeCompressed(); /* NOTE To save the extra space */
- jj++;
- } /* End While loop -- end of the elimination on the row ii*/
- /* Reset the upper part of the pointer jr to zero */
- for (k = 0; k <sizeu; k++){
- jr(ju(ii+k)) = -1;
- }
- /* Sort the L-part of the row --use Quicksplit()*/
- sizel = len;
- len = std::min(sizel, nnzL );
- typename Vector::SegmentReturnType ul(u.segment(0, len));
- typename VectorXi::SegmentReturnType jul(ju.segment(0, len));
- QuickSplit(ul, jul, len);
-
-
- /* Store the largest m_fill elements of the L part */
- m_lu.startVec(ii);
- for (k = 0; k < len; k++){
- m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
+ m_factorizationIsOk = true;
}
- /* Store the diagonal element */
- if (u(ii) == Scalar(0))
- u(ii) = std::sqrt(m_droptol ) * rownorm ; /* NOTE This is used to avoid a zero pivot, because we are doing an incomplete factorization */
- m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
- /* Sort the U-part of the row -- Use Quicksplit() */
- len = 0;
- for (k = 1; k < sizeu; k++) { /* First, drop any element that is below a relative tolerance */
- if ( std::abs(u(ii+k)) > m_droptol * rownorm ) {
- ++len;
- u(ii + len) = u(ii + k);
- ju(ii + len) = ju(ii + k);
- }
+ /**
+ * Compute an incomplete LU factorization with dual threshold on the matrix mat
+ * No pivoting is done in this version
+ *
+ **/
+ template<typename MatrixType>
+ IncompleteLUT<Scalar>& compute(const MatrixType& amat)
+ {
+ analyzePattern(amat);
+ factorize(amat);
+ eigen_assert(m_factorizationIsOk == true);
+ m_isInitialized = true;
+ return *this;
}
- sizeu = len + 1; /* To take into account the diagonal element */
- len = std::min(sizeu, nnzU);
- typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
- typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
- QuickSplit(uu, juu, len);
- /* Store the largest <fill> elements of the U part */
- for (k = ii + 1; k < ii + len; k++){
- m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
- }
- } /* End global for-loop */
- m_lu.finalize();
- m_lu.makeCompressed(); /* NOTE To save the extra space */
- m_isInitialized = true;
- return *this;
-}
void setDroptol(RealScalar droptol);
- void setFill(int fillfactor);
+ void setFillfactor(int fillfactor);
@@ -302,8 +319,10 @@
protected:
FactorType m_lu;
RealScalar m_droptol;
- int m_fillfactor;
- bool m_isInitialized;
+ int m_fillfactor;
+ bool m_factorizationIsOk;
+ bool m_analysisIsOk;
+ bool m_isInitialized;
template <typename VectorV, typename VectorI>
int QuickSplit(VectorV &row, VectorI &ind, int ncut);
PermutationMatrix<Dynamic,Dynamic,Index> m_P; /* Fill-reducing permutation */
@@ -333,7 +352,7 @@
* \param fillfactor This is used to compute the number @p fill_in of largest elements to keep on each row.
**/
template<typename Scalar>
-void IncompleteLUT<Scalar>::setFill(int fillfactor)
+void IncompleteLUT<Scalar>::setFillfactor(int fillfactor)
{
this->m_fillfactor = fillfactor;
}
