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

Cleaning documentation pass in ordering and ILUT

Browse Source
This commit is contained in:
Gael Guennebaud 2013-01-12 11:56:56 +01:00
parent 38fa432e07
commit 7262cf783c
3 changed files with 104 additions and 97 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

78
Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
View File

@ -15,15 +15,15 @@ namespace Eigen {
namespace internal {
/**
* Compute a quick-sort split of a vector
* On output, the vector row is permuted such that its elements satisfy
* abs(row(i)) >= abs(row(ncut)) if i<ncut
* abs(row(i)) <= abs(row(ncut)) if i>ncut
* \param row The vector of values
* \param ind The array of index for the elements in @p row
* \param ncut The number of largest elements to keep
**/
/** \internal
* Compute a quick-sort split of a vector
* On output, the vector row is permuted such that its elements satisfy
* abs(row(i)) >= abs(row(ncut)) if i<ncut
* abs(row(i)) <= abs(row(ncut)) if i>ncut
* \param row The vector of values
* \param ind The array of index for the elements in @p row
* \param ncut The number of largest elements to keep
**/
template <typename VectorV, typename VectorI>
int QuickSplit(VectorV &row, VectorI &ind, int ncut)
{
@ -60,34 +60,37 @@ int QuickSplit(VectorV &row, VectorI &ind, int ncut)
}
}// end namespace internal
/**
* \brief Incomplete LU factorization with dual-threshold strategy
* During the numerical factorization, two dropping rules are used :
* 1) any element whose magnitude is less than some tolerance is dropped.
* This tolerance is obtained by multiplying the input tolerance @p droptol
* by the average magnitude of all the original elements in the current row.
* 2) After the elimination of the row, only the @p fill largest elements in
* the L part and the @p fill largest elements in the U part are kept
* (in addition to the diagonal element ). Note that @p fill is computed from
* the input parameter @p fillfactor which is used the ratio to control the fill_in
* relatively to the initial number of nonzero elements.
*
* The two extreme cases are when @p droptol=0 (to keep all the @p fill*2 largest elements)
* and when @p fill=n/2 with @p droptol being different to zero.
*
* References : Yousef Saad, ILUT: A dual threshold incomplete LU factorization,
* Numerical Linear Algebra with Applications, 1(4), pp 387-402, 1994.
*
* NOTE : The following implementation is derived from the ILUT implementation
* in the SPARSKIT package, Copyright (C) 2005, the Regents of the University of Minnesota
* released under the terms of the GNU LGPL:
* http://www-users.cs.umn.edu/~saad/software/SPARSKIT/README
* However, Yousef Saad gave us permission to relicense his ILUT code to MPL2.
* See the Eigen mailing list archive, thread: ILUT, date: July 8, 2012:
* http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/07/msg00064.html
* alternatively, on GMANE:
* http://comments.gmane.org/gmane.comp.lib.eigen/3302
*/
/** \ingroup IterativeLinearSolvers_Module
* \class IncompleteLUT
* \brief Incomplete LU factorization with dual-threshold strategy
*
* During the numerical factorization, two dropping rules are used :
* 1) any element whose magnitude is less than some tolerance is dropped.
* This tolerance is obtained by multiplying the input tolerance @p droptol
* by the average magnitude of all the original elements in the current row.
* 2) After the elimination of the row, only the @p fill largest elements in
* the L part and the @p fill largest elements in the U part are kept
* (in addition to the diagonal element ). Note that @p fill is computed from
* the input parameter @p fillfactor which is used the ratio to control the fill_in
* relatively to the initial number of nonzero elements.
*
* The two extreme cases are when @p droptol=0 (to keep all the @p fill*2 largest elements)
* and when @p fill=n/2 with @p droptol being different to zero.
*
* References : Yousef Saad, ILUT: A dual threshold incomplete LU factorization,
* Numerical Linear Algebra with Applications, 1(4), pp 387-402, 1994.
*
* NOTE : The following implementation is derived from the ILUT implementation
* in the SPARSKIT package, Copyright (C) 2005, the Regents of the University of Minnesota
* released under the terms of the GNU LGPL:
* http://www-users.cs.umn.edu/~saad/software/SPARSKIT/README
* However, Yousef Saad gave us permission to relicense his ILUT code to MPL2.
* See the Eigen mailing list archive, thread: ILUT, date: July 8, 2012:
* http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/07/msg00064.html
* alternatively, on GMANE:
* http://comments.gmane.org/gmane.comp.lib.eigen/3302
*/
template <typename _Scalar>
class IncompleteLUT : internal::noncopyable
{
@ -462,4 +465,3 @@ struct solve_retval<IncompleteLUT<_MatrixType>, Rhs>
} // end namespace Eigen
#endif // EIGEN_INCOMPLETE_LUT_H

1
Eigen/src/OrderingMethods/Amd.h
View File

@ -86,6 +86,7 @@ Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Ind
/** \internal
* \ingroup OrderingMethods_Module
* Approximate minimum degree ordering algorithm.
* \returns the permutation P reducing the fill-in of the input matrix \a C
* The input matrix \a C must be a selfadjoint compressed column major SparseMatrix object. Both the upper and lower parts have to be stored, but the diagonal entries are optional.

122
Eigen/src/OrderingMethods/Ordering.h
View File

@ -33,30 +33,34 @@ namespace Eigen {
namespace internal {
/**
* Get the symmetric pattern A^T+A from the input matrix A.
* FIXME: The values should not be considered here
*/
template<typename MatrixType>
void ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat)
{
MatrixType C;
C = mat.transpose(); // NOTE: Could be costly
for (int i = 0; i < C.rows(); i++)
{
for (typename MatrixType::InnerIterator it(C, i); it; ++it)
it.valueRef() = 0.0;
}
symmat = C + mat;
}
/** \internal
* \ingroup OrderingMethods_Module
* \returns the symmetric pattern A^T+A from the input matrix A.
* FIXME: The values should not be considered here
*/
template<typename MatrixType>
void ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat)
{
MatrixType C;
C = mat.transpose(); // NOTE: Could be costly
for (int i = 0; i < C.rows(); i++)
{
for (typename MatrixType::InnerIterator it(C, i); it; ++it)
it.valueRef() = 0.0;
}
symmat = C + mat;
}
}
/**
* Get the approximate minimum degree ordering
* If the matrix is not structurally symmetric, an ordering of A^T+A is computed
* \tparam Index The type of indices of the matrix
*/
/** \ingroup OrderingMethods_Module
* \class AMDOrdering
*
* Functor computing the \em approximate \em minimum \em degree ordering
* If the matrix is not structurally symmetric, an ordering of A^T+A is computed
* \tparam Index The type of indices of the matrix
* \sa COLAMDOrdering
*/
template <typename Index>
class AMDOrdering
{
@ -90,12 +94,14 @@ class AMDOrdering
}
};
/**
* Get the natural ordering
*
*NOTE Returns an empty permutation matrix
* \tparam Index The type of indices of the matrix
*/
/** \ingroup OrderingMethods_Module
* \class NaturalOrdering
*
* Functor computing the natural ordering (identity)
*
* \note Returns an empty permutation matrix
* \tparam Index The type of indices of the matrix
*/
template <typename Index>
class NaturalOrdering
{
@ -111,48 +117,46 @@ class NaturalOrdering
};
/**
* Get the column approximate minimum degree ordering
* The matrix should be in column-major format
*/
template<typename Index>
class COLAMDOrdering;
#include "Eigen_Colamd.h"
/** \ingroup OrderingMethods_Module
* \class COLAMDOrdering
*
* Functor computing the \em column \em approximate \em minimum \em degree ordering
* The matrix should be in column-major format
*/
template<typename Index>
class COLAMDOrdering
{
public:
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
typedef Matrix<Index, Dynamic, 1> IndexVector;
typedef Matrix<Index, Dynamic, 1> IndexVector;
/** Compute the permutation vector form a sparse matrix */
template <typename MatrixType>
void operator() (const MatrixType& mat, PermutationType& perm)
{
int m = mat.rows();
int n = mat.cols();
int nnz = mat.nonZeros();
// Get the recommended value of Alen to be used by colamd
int Alen = internal::colamd_recommended(nnz, m, n);
// Set the default parameters
double knobs [COLAMD_KNOBS];
int stats [COLAMD_STATS];
internal::colamd_set_defaults(knobs);
int info;
IndexVector p(n+1), A(Alen);
for(int i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i];
for(int i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i];
// Call Colamd routine to compute the ordering
info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats);
eigen_assert( info && "COLAMD failed " );
perm.resize(n);
for (int i = 0; i < n; i++) perm.indices()(p(i)) = i;
int m = mat.rows();
int n = mat.cols();
int nnz = mat.nonZeros();
// Get the recommended value of Alen to be used by colamd
int Alen = internal::colamd_recommended(nnz, m, n);
// Set the default parameters
double knobs [COLAMD_KNOBS];
int stats [COLAMD_STATS];
internal::colamd_set_defaults(knobs);
int info;
IndexVector p(n+1), A(Alen);
for(int i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i];
for(int i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i];
// Call Colamd routine to compute the ordering
info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats);
eigen_assert( info && "COLAMD failed " );
perm.resize(n);
for (int i = 0; i < n; i++) perm.indices()(p(i)) = i;
}
};
} // end namespace Eigen
#endif
#endif