From 3f3bc6d8622dcc69599c1191e6d0e70aa1436bd3 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 9 Jan 2024 18:18:05 +0000 Subject: [PATCH] Improve documentation of SparseLU --- Eigen/src/SparseLU/SparseLU.h | 145 ++++++++++++++++++++++++++-------- 1 file changed, 112 insertions(+), 33 deletions(-) diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index aee3d947d..29be01a27 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -99,13 +99,34 @@ class SparseLUTransposeView : public SparseSolverBase A; - * SparseLU, COLAMDOrdering > solver; - * // fill A and b; - * // Compute the ordering permutation vector from the structural pattern of A + * SparseLU, COLAMDOrdering > solver; + * // Fill A and b. + * // Compute the ordering permutation vector from the structural pattern of A. * solver.analyzePattern(A); - * // Compute the numerical factorization + * // Compute the numerical factorization. * solver.factorize(A); - * //Use the factors to solve the linear system + * // Use the factors to solve the linear system. + * x = solver.solve(b); + * \endcode + * + * We can directly call compute() instead of analyzePattern() and factorize() + * \code + * VectorXd x(n), b(n); + * SparseMatrix A; + * SparseLU, COLAMDOrdering > solver; + * // Fill A and b. + * solver.compute(A); + * // Use the factors to solve the linear system. + * x = solver.solve(b); + * \endcode + * + * Or give the matrix to the constructor SparseLU(const MatrixType& matrix) + * \code + * VectorXd x(n), b(n); + * SparseMatrix A; + * // Fill A and b. + * SparseLU, COLAMDOrdering > solver(A); + * // Use the factors to solve the linear system. * x = solver.solve(b); * \endcode * @@ -150,10 +171,18 @@ class SparseLU : public SparseSolverBase>, enum { ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; public: + /** \brief Basic constructor of the solver. + * + * Construct a SparseLU. As no matrix is given as argument, compute() should be called afterward with a matrix. + */ SparseLU() : m_lastError(""), m_Ustore(0, 0, 0, 0, 0, 0), m_symmetricmode(false), m_diagpivotthresh(1.0), m_detPermR(1) { initperfvalues(); } + /** \brief Constructor of the solver already based on a specific matrix. + * + * Construct a SparseLU. compute() is already called with the given matrix. + */ explicit SparseLU(const MatrixType& matrix) : m_lastError(""), m_Ustore(0, 0, 0, 0, 0, 0), m_symmetricmode(false), m_diagpivotthresh(1.0), m_detPermR(1) { initperfvalues(); @@ -168,9 +197,15 @@ class SparseLU : public SparseSolverBase>, void factorize(const MatrixType& matrix); void simplicialfactorize(const MatrixType& matrix); - /** + /** \brief Analyze and factorize the matrix so the solver is ready to solve. + * * Compute the symbolic and numeric factorization of the input sparse matrix. - * The input matrix should be in column-major storage. + * The input matrix should be in column-major storage, otherwise analyzePattern() + * will do a heavy copy. + * + * Call analyzePattern() followed by factorize() + * + * \sa analyzePattern(), factorize() */ void compute(const MatrixType& matrix) { // Analyze @@ -179,7 +214,9 @@ class SparseLU : public SparseSolverBase>, factorize(matrix); } - /** \returns an expression of the transposed of the factored matrix. + /** \brief Return a solver for the transposed matrix. + * + * \returns an expression of the transposed of the factored matrix. * * A typical usage is to solve for the transposed problem A^T x = b: * \code @@ -196,7 +233,9 @@ class SparseLU : public SparseSolverBase>, return transposeView; } - /** \returns an expression of the adjoint of the factored matrix + /** \brief Return a solver for the adjointed matrix. + * + * \returns an expression of the adjoint of the factored matrix * * A typical usage is to solve for the adjoint problem A' x = b: * \code @@ -215,19 +254,28 @@ class SparseLU : public SparseSolverBase>, return adjointView; } + /** \brief Give the number of rows. + */ inline Index rows() const { return m_mat.rows(); } + /** \brief Give the numver of columns. + */ inline Index cols() const { return m_mat.cols(); } - /** Indicate that the pattern of the input matrix is symmetric */ + /** \brief Let you set that the pattern of the input matrix is symmetric + */ void isSymmetric(bool sym) { m_symmetricmode = sym; } - /** \returns an expression of the matrix L, internally stored as supernodes + /** \brief Give the matrixL + * + * \returns an expression of the matrix L, internally stored as supernodes * The only operation available with this expression is the triangular solve * \code * y = b; matrixL().solveInPlace(y); * \endcode */ SparseLUMatrixLReturnType matrixL() const { return SparseLUMatrixLReturnType(m_Lstore); } - /** \returns an expression of the matrix U, + /** \brief Give the MatrixU + * + * \returns an expression of the matrix U, * The only operation available with this expression is the triangular solve * \code * y = b; matrixU().solveInPlace(y); @@ -237,12 +285,14 @@ class SparseLU : public SparseSolverBase>, return SparseLUMatrixUReturnType>>(m_Lstore, m_Ustore); } - /** + /** \brief Give the row matrix permutation. + * * \returns a reference to the row matrix permutation \f$ P_r \f$ such that \f$P_r A P_c^T = L U\f$ * \sa colsPermutation() */ inline const PermutationType& rowsPermutation() const { return m_perm_r; } - /** + /** \brief Give the column matrix permutation. + * * \returns a reference to the column matrix permutation\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$ * \sa rowsPermutation() */ @@ -251,7 +301,9 @@ class SparseLU : public SparseSolverBase>, void setPivotThreshold(const RealScalar& thresh) { m_diagpivotthresh = thresh; } #ifdef EIGEN_PARSED_BY_DOXYGEN - /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. + /** \brief Solve a system \f$ A X = B \f$ + * + * \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. * * \warning the destination matrix X in X = this->solve(B) must be colmun-major. * @@ -267,14 +319,17 @@ class SparseLU : public SparseSolverBase>, * \c NumericalIssue if the LU factorization reports a problem, zero diagonal for instance * \c InvalidInput if the input matrix is invalid * - * \sa iparm() + * You can get a readable error message with lastErrorMessage(). + * + * \sa lastErrorMessage() */ ComputationInfo info() const { eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_info; } - /** + /** \brief Give a human readable error + * * \returns A string describing the type of error */ std::string lastErrorMessage() const { return m_lastError; } @@ -302,7 +357,8 @@ class SparseLU : public SparseSolverBase>, return true; } - /** + /** \brief Give the absolute value of the determinant. + * * \returns the absolute value of the determinant of the matrix of which * *this is the QR decomposition. * @@ -330,7 +386,9 @@ class SparseLU : public SparseSolverBase>, return det; } - /** \returns the natural log of the absolute value of the determinant of the matrix + /** \brief Give the natural log of the absolute determinant. + * + * \returns the natural log of the absolute value of the determinant of the matrix * of which **this is the QR decomposition * * \note This method is useful to work around the risk of overflow/underflow that's @@ -356,7 +414,9 @@ class SparseLU : public SparseSolverBase>, return det; } - /** \returns A number representing the sign of the determinant + /** \brief Give the sign of the determinant. + * + * \returns A number representing the sign of the determinant * * \sa absDeterminant(), logAbsDeterminant() */ @@ -380,7 +440,9 @@ class SparseLU : public SparseSolverBase>, return det * m_detPermR * m_detPermC; } - /** \returns The determinant of the matrix. + /** \brief Give the determinant. + * + * \returns The determinant of the matrix. * * \sa absDeterminant(), logAbsDeterminant() */ @@ -401,7 +463,11 @@ class SparseLU : public SparseSolverBase>, return (m_detPermR * m_detPermC) > 0 ? det : -det; } + /** \brief Give the number of non zero in matrix L. + */ Index nnzL() const { return m_nnzL; } + /** \brief Give the number of non zero in matrix U. + */ Index nnzU() const { return m_nnzU; } protected: @@ -442,7 +508,8 @@ class SparseLU : public SparseSolverBase>, }; // End class SparseLU // Functions needed by the anaysis phase -/** +/** \brief Compute the column permutation. + * * Compute the column permutation to minimize the fill-in * * - Apply this permutation to the input matrix - @@ -451,6 +518,11 @@ class SparseLU : public SparseSolverBase>, * * - Postorder the elimination tree and the column permutation * + * It is possible to call compute() instead of analyzePattern() + factorize(). + * + * If the matrix is row-major this function will do an heavy copy. + * + * \sa factorize(), compute() */ template void SparseLU::analyzePattern(const MatrixType& mat) { @@ -516,23 +588,24 @@ void SparseLU::analyzePattern(const MatrixType& mat) { // Functions needed by the numerical factorization phase -/** +/** \brief Factorize the matrix to get the solver ready. + * * - Numerical factorization * - Interleaved with the symbolic factorization - * On exit, info is * - * = 0: successful factorization + * To get error of this function you should check info(), you can get more info of + * errors with lastErrorMessage(). * - * > 0: if info = i, and i is + * In the past (before 2012 (git history is not older)), this function was returning an integer. + * This exit was 0 if successful factorization. + * > 0 if info = i, and i is been completed, but the factor U is exactly singular, + * and division by zero will occur if it is used to solve a system of equation. + * > A->ncol: number of bytes allocated when memory allocation failure occured, plus A->ncol. + * If lwork = -1, it is the estimated amount of space needed, plus A->ncol. * - * <= A->ncol: U(i,i) is exactly zero. The factorization has - * been completed, but the factor U is exactly singular, - * and division by zero will occur if it is used to solve a - * system of equations. + * It seems that A was the name of the matrix in the past. * - * > A->ncol: number of bytes allocated when memory allocation - * failure occurred, plus A->ncol. If lwork = -1, it is - * the estimated amount of space needed, plus A->ncol. + * \sa analyzePattern(), compute(), SparseLU(), info(), lastErrorMessage() */ template void SparseLU::factorize(const MatrixType& matrix) { @@ -572,6 +645,8 @@ void SparseLU::factorize(const MatrixType& matrix) { Index maxpanel = m_perfv.panel_size * m; // Allocate working storage common to the factor routines Index lwork = 0; + // Return the size of actually allocated memory when allocation failed, + // and 0 on success. Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); if (info) { m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n"; @@ -656,6 +731,7 @@ void SparseLU::factorize(const MatrixType& matrix) { // Depth-first-search for the current column VectorBlock panel_lsubk(panel_lsub, k, m); VectorBlock repfnz_k(repfnz, k, m); + // Return 0 on success and > 0 number of bytes allocated when run out of space. info = Base::column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); if (info) { @@ -667,6 +743,7 @@ void SparseLU::factorize(const MatrixType& matrix) { // Numeric updates to this column VectorBlock dense_k(dense, k, m); VectorBlock segrep_k(segrep, nseg1, m - nseg1); + // Return 0 on success and > 0 number of bytes allocated when run out of space. info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); if (info) { m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() "; @@ -676,6 +753,7 @@ void SparseLU::factorize(const MatrixType& matrix) { } // Copy the U-segments to ucol(*) + // Return 0 on success and > 0 number of bytes allocated when run out of space. info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k, m_perm_r.indices(), dense_k, m_glu); if (info) { m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() "; @@ -685,6 +763,7 @@ void SparseLU::factorize(const MatrixType& matrix) { } // Form the L-segment + // Return O if success, i > 0 if U(i, i) is exactly zero. info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); if (info) { m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR";