Clean inclusion, namespace definition, and documentation of SparseLU

This commit is contained in:
Gael Guennebaud 2013-01-12 11:55:16 +01:00
parent 50625834e6
commit 38fa432e07
18 changed files with 235 additions and 145 deletions

View File

@ -1,9 +1,17 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@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
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_SPARSELU_MODULE_H
#define EIGEN_SPARSELU_MODULE_H
#include "SparseCore"
/**
* \defgroup SparseLU_Module SparseLU module
*
@ -14,6 +22,22 @@
#include "src/SparseLU/SparseLU_gemm_kernel.h"
#include "src/SparseLU/SparseLU_Structs.h"
#include "src/SparseLU/SparseLU_Matrix.h"
#include "src/SparseLU/SparseLUBase.h"
#include "src/SparseLU/SparseLU_Coletree.h"
#include "src/SparseLU/SparseLU_Memory.h"
#include "src/SparseLU/SparseLU_heap_relax_snode.h"
#include "src/SparseLU/SparseLU_relax_snode.h"
#include "src/SparseLU/SparseLU_pivotL.h"
#include "src/SparseLU/SparseLU_panel_dfs.h"
#include "src/SparseLU/SparseLU_kernel_bmod.h"
#include "src/SparseLU/SparseLU_panel_bmod.h"
#include "src/SparseLU/SparseLU_column_dfs.h"
#include "src/SparseLU/SparseLU_column_bmod.h"
#include "src/SparseLU/SparseLU_copy_to_ucol.h"
#include "src/SparseLU/SparseLU_pruneL.h"
#include "src/SparseLU/SparseLU_Utils.h"
#include "src/SparseLU/SparseLU.h"
#endif // EIGEN_SPARSELU_MODULE_H

View File

@ -13,63 +13,57 @@
namespace Eigen {
// Data structure needed by all routines
#include "SparseLU_Structs.h"
#include "SparseLU_Matrix.h"
// Base structure containing all the factorization routines
#include "SparseLUBase.h"
/**
* \ingroup SparseLU_Module
* \brief Sparse supernodal LU factorization for general matrices
*
* This class implements the supernodal LU factorization for general matrices.
* It uses the main techniques from the sequential SuperLU package
* (http://crd-legacy.lbl.gov/~xiaoye/SuperLU/). It handles transparently real
* and complex arithmetics with single and double precision, depending on the
* scalar type of your input matrix.
* The code has been optimized to provide BLAS-3 operations during supernode-panel updates.
* It benefits directly from the built-in high-performant Eigen BLAS routines.
* Moreover, when the size of a supernode is very small, the BLAS calls are avoided to
* enable a better optimization from the compiler. For best performance,
* you should compile it with NDEBUG flag to avoid the numerous bounds checking on vectors.
*
* An important parameter of this class is the ordering method. It is used to reorder the columns
* (and eventually the rows) of the matrix to reduce the number of new elements that are created during
* numerical factorization. The cheapest method available is COLAMD.
* See \link OrderingMethods_Module the OrderingMethods module \endlink for the list of
* built-in and external ordering methods.
*
* Simple example with key steps
* \code
* VectorXd x(n), b(n);
* SparseMatrix<double, ColMajor> A;
* SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<int> > solver;
* // fill A and b;
* // Compute the ordering permutation vector from the structural pattern of A
* solver.analyzePattern(A);
* // Compute the numerical factorization
* solver.factorize(A);
* //Use the factors to solve the linear system
* x = solver.solve(b);
* \endcode
*
* \warning The input matrix A should be in a \b compressed and \b column-major form.
* Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix.
*
* \note Unlike the initial SuperLU implementation, there is no step to equilibrate the matrix.
* For badly scaled matrices, this step can be useful to reduce the pivoting during factorization.
* If this is the case for your matrices, you can try the basic scaling method at
* "unsupported/Eigen/src/IterativeSolvers/Scaling.h"
*
* \tparam _MatrixType The type of the sparse matrix. It must be a column-major SparseMatrix<>
* \tparam _OrderingType The ordering method to use, either AMD, COLAMD or METIS
*
*
* \sa \ref TutorialSparseDirectSolvers
* \sa \ref OrderingMethods_Module
*/
/** \ingroup SparseLU_Module
* \class SparseLU
*
* \brief Sparse supernodal LU factorization for general matrices
*
* This class implements the supernodal LU factorization for general matrices.
* It uses the main techniques from the sequential SuperLU package
* (http://crd-legacy.lbl.gov/~xiaoye/SuperLU/). It handles transparently real
* and complex arithmetics with single and double precision, depending on the
* scalar type of your input matrix.
* The code has been optimized to provide BLAS-3 operations during supernode-panel updates.
* It benefits directly from the built-in high-performant Eigen BLAS routines.
* Moreover, when the size of a supernode is very small, the BLAS calls are avoided to
* enable a better optimization from the compiler. For best performance,
* you should compile it with NDEBUG flag to avoid the numerous bounds checking on vectors.
*
* An important parameter of this class is the ordering method. It is used to reorder the columns
* (and eventually the rows) of the matrix to reduce the number of new elements that are created during
* numerical factorization. The cheapest method available is COLAMD.
* See \link OrderingMethods_Module the OrderingMethods module \endlink for the list of
* built-in and external ordering methods.
*
* Simple example with key steps
* \code
* VectorXd x(n), b(n);
* SparseMatrix<double, ColMajor> A;
* SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<int> > solver;
* // fill A and b;
* // Compute the ordering permutation vector from the structural pattern of A
* solver.analyzePattern(A);
* // Compute the numerical factorization
* solver.factorize(A);
* //Use the factors to solve the linear system
* x = solver.solve(b);
* \endcode
*
* \warning The input matrix A should be in a \b compressed and \b column-major form.
* Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix.
*
* \note Unlike the initial SuperLU implementation, there is no step to equilibrate the matrix.
* For badly scaled matrices, this step can be useful to reduce the pivoting during factorization.
* If this is the case for your matrices, you can try the basic scaling method at
* "unsupported/Eigen/src/IterativeSolvers/Scaling.h"
*
* \tparam _MatrixType The type of the sparse matrix. It must be a column-major SparseMatrix<>
* \tparam _OrderingType The ordering method to use, either AMD, COLAMD or METIS
*
*
* \sa \ref TutorialSparseDirectSolvers
* \sa \ref OrderingMethods_Module
*/
template <typename _MatrixType, typename _OrderingType>
class SparseLU
{
@ -548,7 +542,6 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
m_factorizationIsOk = true;
}
// #include "SparseLU_simplicialfactorize.h"
namespace internal {
template<typename _MatrixType, typename Derived, typename Rhs>

View File

@ -8,9 +8,13 @@
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef SPARSELUBASE_H
#define SPARSELUBASE_H
/**
* Base class for sparseLU
*/
namespace Eigen {
/** \ingroup SparseLU_Module
* \class SparseLUBase
* Base class for sparseLU
*/
template <typename Scalar, typename Index>
struct SparseLUBase
{
@ -49,18 +53,6 @@ struct SparseLUBase
};
#include "SparseLU_Coletree.h"
#include "SparseLU_Memory.h"
#include "SparseLU_heap_relax_snode.h"
#include "SparseLU_relax_snode.h"
#include "SparseLU_pivotL.h"
#include "SparseLU_panel_dfs.h"
#include "SparseLU_kernel_bmod.h"
#include "SparseLU_panel_bmod.h"
#include "SparseLU_column_dfs.h"
#include "SparseLU_column_bmod.h"
#include "SparseLU_copy_to_ucol.h"
#include "SparseLU_pruneL.h"
#include "SparseLU_Utils.h"
} // end namespace Eigen
#endif

View File

@ -31,7 +31,10 @@
#ifndef SPARSELU_COLETREE_H
#define SPARSELU_COLETREE_H
namespace Eigen {
namespace internal {
/** Find the root of the tree/set containing the vertex i : Use Path halving */
template<typename IndexVector>
int etree_find (int i, IndexVector& pp)
@ -187,5 +190,8 @@ void treePostorder(int n, IndexVector& parent, IndexVector& post)
internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
}
} // internal
#endif
} // end namespace internal
} // end namespace Eigen
#endif // SPARSELU_COLETREE_H

View File

@ -11,6 +11,8 @@
#ifndef EIGEN_SPARSELU_MATRIX_H
#define EIGEN_SPARSELU_MATRIX_H
namespace Eigen {
/** \ingroup SparseLU_Module
* \brief a class to manipulate the L supernodal factor from the SparseLU factorization
*
@ -309,5 +311,6 @@ void SuperNodalMatrix<Scalar,Index>::solveInPlace( MatrixBase<Dest>&X) const
}
}
} // end namespace Eigen
#endif
#endif // EIGEN_SPARSELU_MATRIX_H

View File

@ -31,6 +31,8 @@
#ifndef EIGEN_SPARSELU_MEMORY
#define EIGEN_SPARSELU_MEMORY
namespace Eigen {
#define LU_NO_MARKER 3
#define LU_NUM_TEMPV(m,w,t,b) ((std::max)(m, (t+b)*w) )
#define IND_EMPTY (-1)
@ -198,7 +200,9 @@ int SparseLUBase<Scalar,Index>::LUMemXpand(VectorType& vec, int& maxlen, int nbE
if (failed_size)
return failed_size;
return 0 ;
return 0 ;
}
#endif
} // end namespace Eigen
#endif // EIGEN_SPARSELU_MEMORY

View File

@ -65,10 +65,13 @@
* xusub[i] points to the starting location of column i in ucol.
* Storage: new row subscripts; that is subscripts of PA.
*/
#ifndef EIGEN_LU_STRUCTS
#define EIGEN_LU_STRUCTS
typedef enum {LUSUP, UCOL, LSUB, USUB, LLVL, ULVL} LU_MemType;
namespace Eigen {
typedef enum {LUSUP, UCOL, LSUB, USUB, LLVL, ULVL} LU_MemType;
template <typename IndexVector, typename ScalarVector>
struct LU_GlobalLU_t {
@ -100,4 +103,7 @@ struct LU_perfvalues {
int colblk; // The minimum column dimension for 2-D blocking to be used;
int fillfactor; // The estimated fills factors for L and U, compared with A
};
#endif
} // end namespace Eigen
#endif // EIGEN_LU_STRUCTS

View File

@ -11,6 +11,7 @@
#ifndef EIGEN_SPARSELU_UTILS_H
#define EIGEN_SPARSELU_UTILS_H
namespace Eigen {
/**
* \brief Count Nonzero elements in the factors
@ -37,8 +38,8 @@ void SparseLUBase<Scalar,Index>::LU_countnz(const int n, int& nnzL, int& nnzU, G
jlen--;
}
}
}
/**
* \brief Fix up the data storage lsub for L-subscripts.
*
@ -72,4 +73,6 @@ void SparseLUBase<Scalar,Index>::LU_fixupL(const int n, const IndexVector& perm_
glu.xlsub(n) = nextl;
}
#endif
} // end namespace Eigen
#endif // EIGEN_SPARSELU_UTILS_H

View File

@ -31,6 +31,8 @@
#ifndef SPARSELU_COLUMN_BMOD_H
#define SPARSELU_COLUMN_BMOD_H
namespace Eigen {
/**
* \brief Performs numeric block updates (sup-col) in topological order
*
@ -170,4 +172,7 @@ int SparseLUBase<Scalar,Index>::LU_column_bmod(const int jcol, const int nseg, B
} // End if fst_col
return 0;
}
#endif
} // end namespace Eigen
#endif // SPARSELU_COLUMN_BMOD_H

View File

@ -29,6 +29,38 @@
*/
#ifndef SPARSELU_COLUMN_DFS_H
#define SPARSELU_COLUMN_DFS_H
namespace Eigen {
namespace internal {
template<typename IndexVector, typename ScalarVector>
struct LU_column_dfs_traits
{
typedef typename IndexVector::Scalar Index;
typedef typename ScalarVector::Scalar Scalar;
LU_column_dfs_traits(Index jcol, Index& jsuper, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
: m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu)
{}
bool update_segrep(Index /*krep*/, Index /*jj*/)
{
return true;
}
void mem_expand(IndexVector& lsub, int& nextl, int chmark)
{
if (nextl >= m_glu.nzlmax)
SparseLUBase<Scalar,Index>::LUMemXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions);
if (chmark != (m_jcol-1)) m_jsuper_ref = IND_EMPTY;
}
enum { ExpandMem = true };
int m_jcol;
int& m_jsuper_ref;
LU_GlobalLU_t<IndexVector, ScalarVector>& m_glu;
};
} // end namespace internal
/**
* \brief Performs a symbolic factorization on column jcol and decide the supernode boundary
*
@ -56,31 +88,6 @@
* > 0 number of bytes allocated when run out of space
*
*/
template<typename IndexVector, typename ScalarVector>
struct LU_column_dfs_traits
{
typedef typename IndexVector::Scalar Index;
typedef typename ScalarVector::Scalar Scalar;
LU_column_dfs_traits(Index jcol, Index& jsuper, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
: m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu)
{}
bool update_segrep(Index /*krep*/, Index /*jj*/)
{
return true;
}
void mem_expand(IndexVector& lsub, int& nextl, int chmark)
{
if (nextl >= m_glu.nzlmax)
SparseLUBase<Scalar,Index>::LUMemXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions);
if (chmark != (m_jcol-1)) m_jsuper_ref = IND_EMPTY;
}
enum { ExpandMem = true };
int m_jcol;
int& m_jsuper_ref;
LU_GlobalLU_t<IndexVector, ScalarVector>& m_glu;
};
template <typename Scalar, typename Index>
int SparseLUBase<Scalar,Index>::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, int maxsuper, int& nseg, BlockIndexVector lsub_col, IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
{
@ -90,7 +97,7 @@ int SparseLUBase<Scalar,Index>::LU_column_dfs(const int m, const int jcol, Index
VectorBlock<IndexVector> marker2(marker, 2*m, m);
LU_column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu);
internal::LU_column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu);
// For each nonzero in A(*,jcol) do dfs
for (int k = 0; lsub_col[k] != IND_EMPTY; k++)
@ -161,4 +168,7 @@ int SparseLUBase<Scalar,Index>::LU_column_dfs(const int m, const int jcol, Index
return 0;
}
#endif
} // end namespace Eigen
#endif

View File

@ -29,6 +29,8 @@
#ifndef SPARSELU_COPY_TO_UCOL_H
#define SPARSELU_COPY_TO_UCOL_H
namespace Eigen {
/**
* \brief Performs numeric block updates (sup-col) in topological order
*
@ -97,4 +99,6 @@ int SparseLUBase<Scalar,Index>::LU_copy_to_ucol(const int jcol, const int nseg,
return 0;
}
#endif
} // end namespace Eigen
#endif // SPARSELU_COPY_TO_UCOL_H

View File

@ -27,7 +27,9 @@
#ifndef SPARSELU_HEAP_RELAX_SNODE_H
#define SPARSELU_HEAP_RELAX_SNODE_H
#include "SparseLU_Coletree.h"
namespace Eigen {
/**
* \brief Identify the initial relaxed supernodes
*
@ -116,4 +118,7 @@ void SparseLUBase<Scalar,Index>::LU_heap_relax_snode (const int n, IndexVector&
// Recover the original etree
et = et_save;
}
#endif
} // end namespace Eigen
#endif // SPARSELU_HEAP_RELAX_SNODE_H

View File

@ -11,6 +11,8 @@
#ifndef SPARSELU_KERNEL_BMOD_H
#define SPARSELU_KERNEL_BMOD_H
namespace Eigen {
/**
* \brief Performs numeric block updates from a given supernode to a single column
*
@ -108,4 +110,7 @@ template <> struct LU_kernel_bmod<1>
dense.coeffRef(*(irow++)) -= f * *(a++);
}
};
#endif
} // end namespace Eigen
#endif // SPARSELU_KERNEL_BMOD_H

View File

@ -31,6 +31,8 @@
#ifndef SPARSELU_PANEL_BMOD_H
#define SPARSELU_PANEL_BMOD_H
namespace Eigen {
/**
* \brief Performs numeric block updates (sup-panel) in topological order.
*
@ -211,4 +213,7 @@ void SparseLUBase<Scalar,Index>::LU_panel_bmod(const int m, const int w, const i
} // End for each updating supernode
}
#endif
} // end namespace Eigen
#endif // SPARSELU_PANEL_BMOD_H

View File

@ -29,6 +29,35 @@
*/
#ifndef SPARSELU_PANEL_DFS_H
#define SPARSELU_PANEL_DFS_H
namespace Eigen {
namespace internal {
template<typename IndexVector>
struct LU_panel_dfs_traits
{
typedef typename IndexVector::Scalar Index;
LU_panel_dfs_traits(Index jcol, Index* marker)
: m_jcol(jcol), m_marker(marker)
{}
bool update_segrep(Index krep, Index jj)
{
if(m_marker[krep]<m_jcol)
{
m_marker[krep] = jj;
return true;
}
return false;
}
void mem_expand(IndexVector& /*glu.lsub*/, int /*nextl*/, int /*chmark*/) {}
enum { ExpandMem = false };
Index m_jcol;
Index* m_marker;
};
} // end namespace internal
template <typename Scalar, typename Index>
template <typename Traits>
void SparseLUBase<Scalar,Index>::LU_dfs_kernel(const int jj, IndexVector& perm_r,
@ -185,28 +214,6 @@ void SparseLUBase<Scalar,Index>::LU_dfs_kernel(const int jj, IndexVector& perm_r
*
*/
template<typename IndexVector>
struct LU_panel_dfs_traits
{
typedef typename IndexVector::Scalar Index;
LU_panel_dfs_traits(Index jcol, Index* marker)
: m_jcol(jcol), m_marker(marker)
{}
bool update_segrep(Index krep, Index jj)
{
if(m_marker[krep]<m_jcol)
{
m_marker[krep] = jj;
return true;
}
return false;
}
void mem_expand(IndexVector& /*glu.lsub*/, int /*nextl*/, int /*chmark*/) {}
enum { ExpandMem = false };
Index m_jcol;
Index* m_marker;
};
template <typename Scalar, typename Index>
void SparseLUBase<Scalar,Index>::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, IndexVector& perm_r, int& nseg, ScalarVector& dense, IndexVector& panel_lsub, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
{
@ -216,7 +223,7 @@ void SparseLUBase<Scalar,Index>::LU_panel_dfs(const int m, const int w, const in
VectorBlock<IndexVector> marker1(marker, m, m);
nseg = 0;
LU_panel_dfs_traits<IndexVector> traits(jcol, marker1.data());
internal::LU_panel_dfs_traits<IndexVector> traits(jcol, marker1.data());
// For each column in the panel
for (int jj = jcol; jj < jcol + w; jj++)
@ -242,6 +249,8 @@ void SparseLUBase<Scalar,Index>::LU_panel_dfs(const int m, const int w, const in
}// end for nonzeros in column jj
} // end for column jj
}
#endif
} // end namespace Eigen
#endif // SPARSELU_PANEL_DFS_H

View File

@ -29,6 +29,9 @@
*/
#ifndef SPARSELU_PIVOTL_H
#define SPARSELU_PIVOTL_H
namespace Eigen {
/**
* \brief Performs the numerical pivotin on the current column of L, and the CDIV operation.
*
@ -123,4 +126,7 @@ int SparseLUBase<Scalar,Index>::LU_pivotL(const int jcol, const RealScalar diagp
lu_col_ptr[k] *= temp;
return 0;
}
#endif
} // end namespace Eigen
#endif // SPARSELU_PIVOTL_H

View File

@ -30,6 +30,8 @@
#ifndef SPARSELU_PRUNEL_H
#define SPARSELU_PRUNEL_H
namespace Eigen {
/**
* \brief Prunes the L-structure.
*
@ -126,4 +128,6 @@ void SparseLUBase<Scalar,Index>::LU_pruneL(const int jcol, const IndexVector& pe
} // End for each U-segment
}
#endif
} // end namespace Eigen
#endif // SPARSELU_PRUNEL_H

View File

@ -27,6 +27,9 @@
#ifndef SPARSELU_RELAX_SNODE_H
#define SPARSELU_RELAX_SNODE_H
namespace Eigen {
/**
* \brief Identify the initial relaxed supernodes
*
@ -70,4 +73,7 @@ void SparseLUBase<Scalar,Index>::LU_relax_snode (const int n, IndexVector& et, c
} // End postorder traversal of the etree
}
} // end namespace Eigen
#endif