From 5010033d88ee9cab2df4fbdd7b023e9d80256c98 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 17 Jan 2011 11:11:22 +0100 Subject: [PATCH] do not stop the factorization if one pivot is exactly 0, and return the index of the first zero pivot if any --- Eigen/src/LU/PartialPivLU.h | 82 +++++++++++++++++-------------------- 1 file changed, 37 insertions(+), 45 deletions(-) diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 908c00c5c..8694abae7 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -225,7 +225,7 @@ PartialPivLU::PartialPivLU(const MatrixType& matrix) namespace internal { /** \internal This is the blocked version of fullpivlu_unblocked() */ -template +template struct partial_lu_impl { // FIXME add a stride to Map, so that the following mapping becomes easier, @@ -247,51 +247,50 @@ struct partial_lu_impl * of columns of the matrix \a lu, and an integer \a nb_transpositions * which returns the actual number of transpositions. * - * \returns false if some pivot is exactly zero, in which case the matrix is left with - * undefined coefficients (to avoid generating inf/nan values). Returns true - * otherwise. + * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise. */ - static bool unblocked_lu(MatrixType& lu, Index* row_transpositions, Index& nb_transpositions) + static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions) { const Index rows = lu.rows(); - const Index size = std::min(lu.rows(),lu.cols()); + const Index cols = lu.cols(); + const Index size = std::min(rows,cols); nb_transpositions = 0; + int first_zero_pivot = -1; for(Index k = 0; k < size; ++k) { + Index rrows = rows-k-1; + Index rcols = cols-k-1; + Index row_of_biggest_in_col; RealScalar biggest_in_corner = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col); row_of_biggest_in_col += k; - if(biggest_in_corner == 0) // the pivot is exactly zero: the matrix is singular - { - // end quickly, avoid generating inf/nan values. Although in this unblocked_lu case - // the result is still valid, there's no need to boast about it because - // the blocked_lu code can't guarantee the same. - // before exiting, make sure to initialize the still uninitialized row_transpositions - // in a sane state without destroying what we already have. - for(Index i = k; i < size; i++) - row_transpositions[i] = i; - return false; - } - row_transpositions[k] = row_of_biggest_in_col; - if(k != row_of_biggest_in_col) + if(biggest_in_corner != 0) { - lu.row(k).swap(lu.row(row_of_biggest_in_col)); - ++nb_transpositions; + if(k != row_of_biggest_in_col) + { + lu.row(k).swap(lu.row(row_of_biggest_in_col)); + ++nb_transpositions; + } + + // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k) + // overflow but not the actual quotient? + lu.col(k).tail(rrows) /= lu.coeff(k,k); + } + else if(first_zero_pivot==-1) + { + // the pivot is exactly zero, we record the index of the first pivot which is exactly 0, + // and continue the factorization such we still have A = PLU + first_zero_pivot = k; } if(k > > */ - static bool blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, Index* row_transpositions, Index& nb_transpositions, Index maxBlockSize=256) + static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256) { MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols); MatrixType lu(lu1,0,0,rows,cols); @@ -334,6 +331,7 @@ struct partial_lu_impl } nb_transpositions = 0; + int first_zero_pivot = -1; for(Index k = 0; k < size; k+=blockSize) { Index bs = std::min(size-k,blockSize); // actual size of the block @@ -351,21 +349,15 @@ struct partial_lu_impl BlockType A21(lu,k+bs,k,trows,bs); BlockType A22(lu,k+bs,k+bs,trows,tsize); - Index nb_transpositions_in_panel; + PivIndex nb_transpositions_in_panel; // recursively call the blocked LU algorithm on [A11^T A21^T]^T // with a very small blocking size: - if(!blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride, - row_transpositions+k, nb_transpositions_in_panel, 16)) - { - // end quickly with undefined coefficients, just avoid generating inf/nan values. - // before exiting, make sure to initialize the still uninitialized row_transpositions - // in a sane state without destroying what we already have. - for(Index i=k; i=0 && first_zero_pivot==-1) + first_zero_pivot = k+ret; + nb_transpositions += nb_transpositions_in_panel; // update permutations and apply them to A_0 for(Index i=k; i