mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-03-01 18:26:24 +08:00
Define sparseLU functions as static
This commit is contained in:
Desire NUENTSA
parent
7740127e3d
commit
a01371548d
@ -18,6 +18,8 @@ namespace Eigen {
|
||||
#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
|
||||
@ -40,6 +42,7 @@ class SparseLU
|
||||
typedef Matrix<Scalar,Dynamic,1> ScalarVector;
|
||||
typedef Matrix<Index,Dynamic,1> IndexVector;
|
||||
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
|
||||
|
||||
public:
|
||||
SparseLU():m_isInitialized(true),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0)
|
||||
{
|
||||
@ -58,6 +61,7 @@ class SparseLU
|
||||
|
||||
void analyzePattern (const MatrixType& matrix);
|
||||
void factorize (const MatrixType& matrix);
|
||||
void simplicialfactorize(const MatrixType& matrix);
|
||||
|
||||
/**
|
||||
* Compute the symbolic and numeric factorization of the input sparse matrix.
|
||||
@ -224,8 +228,7 @@ class SparseLU
|
||||
PermutationType m_perm_r ; // Row permutation
|
||||
IndexVector m_etree; // Column elimination tree
|
||||
|
||||
LU_GlobalLU_t<IndexVector, ScalarVector> m_glu; // persistent data to facilitate multiple factors
|
||||
// FIXME All fields of this struct can be defined separately as class members
|
||||
LU_GlobalLU_t<IndexVector, ScalarVector> m_glu;
|
||||
|
||||
// SuperLU/SparseLU options
|
||||
bool m_symmetricmode;
|
||||
@ -243,7 +246,6 @@ class SparseLU
|
||||
|
||||
|
||||
// Functions needed by the anaysis phase
|
||||
#include "SparseLU_Coletree.h"
|
||||
/**
|
||||
* Compute the column permutation to minimize the fill-in (file amd.c )
|
||||
*
|
||||
@ -262,9 +264,6 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
|
||||
|
||||
OrderingType ord;
|
||||
ord(mat,m_perm_c);
|
||||
//FIXME Check the right semantic behind m_perm_c
|
||||
// that is, column j of mat goes to column m_perm_c(j) of mat * m_perm_c;
|
||||
|
||||
|
||||
// Apply the permutation to the column of the input matrix
|
||||
// m_mat = mat * m_perm_c.inverse(); //FIXME It should be less expensive here to permute only the structural pattern of the matrix
|
||||
@ -282,13 +281,13 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
|
||||
// Compute the column elimination tree of the permuted matrix
|
||||
/*if (m_etree.size() == 0) */m_etree.resize(m_mat.cols());
|
||||
|
||||
LU_sp_coletree(m_mat, m_etree);
|
||||
SparseLUBase<Scalar,Index>::LU_sp_coletree(m_mat, m_etree);
|
||||
|
||||
// In symmetric mode, do not do postorder here
|
||||
if (!m_symmetricmode) {
|
||||
IndexVector post, iwork;
|
||||
// Post order etree
|
||||
LU_TreePostorder(m_mat.cols(), m_etree, post);
|
||||
SparseLUBase<Scalar,Index>::LU_TreePostorder(m_mat.cols(), m_etree, post);
|
||||
|
||||
|
||||
// Renumber etree in postorder
|
||||
@ -310,21 +309,7 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
|
||||
m_analysisIsOk = true;
|
||||
}
|
||||
|
||||
// Functions needed by the numerical factorization phase
|
||||
#include "SparseLU_Memory.h"
|
||||
#include "SparseLU_heap_relax_snode.h"
|
||||
#include "SparseLU_relax_snode.h"
|
||||
#include "SparseLU_snode_dfs.h"
|
||||
#include "SparseLU_snode_bmod.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"
|
||||
// Functions needed by the numerical factorization phase
|
||||
|
||||
|
||||
/**
|
||||
@ -370,7 +355,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
int maxpanel = m_perfv.panel_size * m;
|
||||
// Allocate working storage common to the factor routines
|
||||
int lwork = 0;
|
||||
int info = LUMemInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu);
|
||||
int info = SparseLUBase<Scalar,Index>::LUMemInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu);
|
||||
if (info)
|
||||
{
|
||||
std::cerr << "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
|
||||
@ -402,25 +387,17 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
// Identify initial relaxed snodes
|
||||
IndexVector relax_end(n);
|
||||
if ( m_symmetricmode == true )
|
||||
LU_heap_relax_snode<IndexVector>(n, m_etree, m_perfv.relax, marker, relax_end);
|
||||
SparseLUBase<Scalar,Index>::LU_heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
|
||||
else
|
||||
LU_relax_snode<IndexVector>(n, m_etree, m_perfv.relax, marker, relax_end);
|
||||
SparseLUBase<Scalar,Index>::LU_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
|
||||
|
||||
|
||||
m_perm_r.resize(m);
|
||||
m_perm_r.indices().setConstant(-1);
|
||||
marker.setConstant(-1);
|
||||
|
||||
IndexVector& xsup = m_glu.xsup;
|
||||
IndexVector& supno = m_glu.supno;
|
||||
IndexVector& xlsub = m_glu.xlsub;
|
||||
IndexVector& xlusup = m_glu.xlusup;
|
||||
IndexVector& xusub = m_glu.xusub;
|
||||
ScalarVector& lusup = m_glu.lusup;
|
||||
Index& nzlumax = m_glu.nzlumax;
|
||||
|
||||
supno(0) = IND_EMPTY; xsup.setConstant(0);
|
||||
xsup(0) = xlsub(0) = xusub(0) = xlusup(0) = Index(0);
|
||||
m_glu.supno(0) = IND_EMPTY; m_glu.xsup.setConstant(0);
|
||||
m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
|
||||
|
||||
// Work on one 'panel' at a time. A panel is one of the following :
|
||||
// (a) a relaxed supernode at the bottom of the etree, or
|
||||
@ -441,7 +418,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
|
||||
// Factorize the relaxed supernode(jcol:kcol)
|
||||
// First, determine the union of the row structure of the snode
|
||||
info = LU_snode_dfs(jcol, kcol, m_mat, xprune, marker, m_glu);
|
||||
info = SparseLUBase<Scalar,Index>::LU_snode_dfs(jcol, kcol, m_mat, xprune, marker, m_glu);
|
||||
if ( info )
|
||||
{
|
||||
std::cerr << "MEMORY ALLOCATION FAILED IN SNODE_DFS() \n";
|
||||
@ -449,15 +426,15 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
m_factorizationIsOk = false;
|
||||
return;
|
||||
}
|
||||
nextu = xusub(jcol); //starting location of column jcol in ucol
|
||||
nextlu = xlusup(jcol); //Starting location of column jcol in lusup (rectangular supernodes)
|
||||
jsupno = supno(jcol); // Supernode number which column jcol belongs to
|
||||
fsupc = xsup(jsupno); //First column number of the current supernode
|
||||
new_next = nextlu + (xlsub(fsupc+1)-xlsub(fsupc)) * (kcol - jcol + 1);
|
||||
nextu = m_glu.xusub(jcol); //starting location of column jcol in ucol
|
||||
nextlu = m_glu.xlusup(jcol); //Starting location of column jcol in lusup (rectangular supernodes)
|
||||
jsupno = m_glu.supno(jcol); // Supernode number which column jcol belongs to
|
||||
fsupc = m_glu.xsup(jsupno); //First column number of the current supernode
|
||||
new_next = nextlu + (m_glu.xlsub(fsupc+1)-m_glu.xlsub(fsupc)) * (kcol - jcol + 1);
|
||||
int mem;
|
||||
while (new_next > nzlumax )
|
||||
while (new_next > m_glu.nzlumax )
|
||||
{
|
||||
mem = LUMemXpand(lusup, nzlumax, nextlu, LUSUP, m_glu.num_expansions);
|
||||
mem = SparseLUBase<Scalar,Index>::LUMemXpand(m_glu.lusup, m_glu.nzlumax, nextlu, LUSUP, m_glu.num_expansions);
|
||||
if (mem)
|
||||
{
|
||||
std::cerr << "MEMORY ALLOCATION FAILED FOR L FACTOR \n";
|
||||
@ -468,16 +445,16 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
|
||||
// Now, left-looking factorize each column within the snode
|
||||
for (icol = jcol; icol<=kcol; icol++){
|
||||
xusub(icol+1) = nextu;
|
||||
m_glu.xusub(icol+1) = nextu;
|
||||
// Scatter into SPA dense(*)
|
||||
for (typename MatrixType::InnerIterator it(m_mat, icol); it; ++it)
|
||||
dense(it.row()) = it.value();
|
||||
|
||||
// Numeric update within the snode
|
||||
LU_snode_bmod(icol, fsupc, dense, m_glu);
|
||||
SparseLUBase<Scalar,Index>::LU_snode_bmod(icol, fsupc, dense, m_glu);
|
||||
|
||||
// Eliminate the current column
|
||||
info = LU_pivotL(icol, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
|
||||
info = SparseLUBase<Scalar,Index>::LU_pivotL(icol, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
|
||||
if ( info )
|
||||
{
|
||||
m_info = NumericalIssue;
|
||||
@ -505,10 +482,10 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
panel_size = n - jcol;
|
||||
|
||||
// Symbolic outer factorization on a panel of columns
|
||||
LU_panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu);
|
||||
SparseLUBase<Scalar,Index>::LU_panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu);
|
||||
|
||||
// Numeric sup-panel updates in topological order
|
||||
LU_panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_perfv, m_glu);
|
||||
SparseLUBase<Scalar,Index>::LU_panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_perfv, m_glu);
|
||||
|
||||
// Sparse LU within the panel, and below the panel diagonal
|
||||
for ( jj = jcol; jj< jcol + panel_size; jj++)
|
||||
@ -519,7 +496,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
//Depth-first-search for the current column
|
||||
VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m);
|
||||
VectorBlock<IndexVector> repfnz_k(repfnz, k, m);
|
||||
info = LU_column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu);
|
||||
info = SparseLUBase<Scalar,Index>::LU_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 )
|
||||
{
|
||||
std::cerr << "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() \n";
|
||||
@ -530,7 +507,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
// Numeric updates to this column
|
||||
VectorBlock<ScalarVector> dense_k(dense, k, m);
|
||||
VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1);
|
||||
info = LU_column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
|
||||
info = SparseLUBase<Scalar,Index>::LU_column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
|
||||
if ( info )
|
||||
{
|
||||
std::cerr << "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() \n";
|
||||
@ -540,7 +517,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
}
|
||||
|
||||
// Copy the U-segments to ucol(*)
|
||||
info = LU_copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu);
|
||||
info = SparseLUBase<Scalar,Index>::LU_copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu);
|
||||
if ( info )
|
||||
{
|
||||
std::cerr << "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() \n";
|
||||
@ -550,7 +527,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
}
|
||||
|
||||
// Form the L-segment
|
||||
info = LU_pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
|
||||
info = SparseLUBase<Scalar,Index>::LU_pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
|
||||
if ( info )
|
||||
{
|
||||
std::cerr<< "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT " << info <<std::endl;
|
||||
@ -560,7 +537,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
}
|
||||
|
||||
// Prune columns (0:jj-1) using column jj
|
||||
LU_pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
|
||||
SparseLUBase<Scalar,Index>::LU_pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
|
||||
|
||||
// Reset repfnz for this column
|
||||
for (i = 0; i < nseg; i++)
|
||||
@ -574,11 +551,9 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
} // end for -- end elimination
|
||||
|
||||
// Count the number of nonzeros in factors
|
||||
LU_countnz(n, m_nnzL, m_nnzU, m_glu);
|
||||
SparseLUBase<Scalar,Index>::LU_countnz(n, m_nnzL, m_nnzU, m_glu);
|
||||
// Apply permutation to the L subscripts
|
||||
LU_fixupL(n, m_perm_r.indices(), m_glu);
|
||||
|
||||
|
||||
SparseLUBase<Scalar,Index>::LU_fixupL(n, m_perm_r.indices(), m_glu);
|
||||
|
||||
// Create supernode matrix L
|
||||
m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup);
|
||||
@ -589,7 +564,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
|
||||
m_factorizationIsOk = true;
|
||||
}
|
||||
|
||||
|
||||
// #include "SparseLU_simplicialfactorize.h"
|
||||
namespace internal {
|
||||
|
||||
template<typename _MatrixType, typename Derived, typename Rhs>
|
||||
@ -607,7 +582,5 @@ struct solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
|
||||
|
||||
} // end namespace internal
|
||||
|
||||
|
||||
|
||||
} // End namespace Eigen
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -31,8 +31,8 @@
|
||||
#ifndef SPARSELU_COLETREE_H
|
||||
#define SPARSELU_COLETREE_H
|
||||
/** Find the root of the tree/set containing the vertex i : Use Path halving */
|
||||
template<typename IndexVector>
|
||||
int etree_find (int i, IndexVector& pp)
|
||||
template< typename Scalar,typename Index>
|
||||
int SparseLUBase<Scalar,Index>::etree_find (int i, IndexVector& pp)
|
||||
{
|
||||
int p = pp(i); // Parent
|
||||
int gp = pp(p); // Grand parent
|
||||
@ -50,8 +50,8 @@ int etree_find (int i, IndexVector& pp)
|
||||
* NOTE : The matrix is supposed to be in column-major format.
|
||||
*
|
||||
*/
|
||||
template<typename MatrixType, typename IndexVector>
|
||||
int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
|
||||
template <typename Scalar, typename Index>
|
||||
int SparseLUBase<Scalar,Index>::LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
|
||||
{
|
||||
int nc = mat.cols(); // Number of columns
|
||||
int nr = mat.rows(); // Number of rows
|
||||
@ -106,8 +106,8 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
|
||||
* Depth-first search from vertex n. No recursion.
|
||||
* This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
|
||||
*/
|
||||
template<typename IndexVector>
|
||||
void LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum)
|
||||
template <typename Scalar, typename Index>
|
||||
void SparseLUBase<Scalar,Index>::LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum)
|
||||
{
|
||||
int current = n, first, next;
|
||||
while (postnum != n)
|
||||
@ -152,8 +152,8 @@ void LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVecto
|
||||
* \param parent Input tree
|
||||
* \param post postordered tree
|
||||
*/
|
||||
template<typename IndexVector>
|
||||
void LU_TreePostorder(int n, IndexVector& parent, IndexVector& post)
|
||||
template <typename Scalar, typename Index>
|
||||
void SparseLUBase<Scalar,Index>::LU_TreePostorder(int n, IndexVector& parent, IndexVector& post)
|
||||
{
|
||||
IndexVector first_kid, next_kid; // Linked list of children
|
||||
int postnum;
|
||||
|
||||
@ -49,8 +49,9 @@
|
||||
* \param keep_prev 1: use length and do not expand the vector; 0: compute new_len and expand
|
||||
* \param [in,out]num_expansions Number of times the memory has been expanded
|
||||
*/
|
||||
template <typename VectorType >
|
||||
int expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_expansions)
|
||||
template <typename Scalar, typename Index>
|
||||
template <typename VectorType>
|
||||
int SparseLUBase<Scalar,Index>::expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_expansions)
|
||||
{
|
||||
|
||||
float alpha = 1.5; // Ratio of the memory increase
|
||||
@ -122,12 +123,9 @@ int expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_ex
|
||||
* \return an estimated size of the required memory if lwork = -1; otherwise, return the size of actually allocated memory when allocation failed, and 0 on success
|
||||
* NOTE Unlike SuperLU, this routine does not support successive factorization with the same pattern and the same row permutation
|
||||
*/
|
||||
template <typename IndexVector,typename ScalarVector>
|
||||
int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size, LU_GlobalLU_t<IndexVector,ScalarVector>& glu)
|
||||
template <typename Scalar, typename Index>
|
||||
int SparseLUBase<Scalar,Index>::LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size, GlobalLU_t& glu)
|
||||
{
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
typedef typename IndexVector::Scalar Index;
|
||||
|
||||
int& num_expansions = glu.num_expansions; //No memory expansions so far
|
||||
num_expansions = 0;
|
||||
glu.nzumax = glu.nzlumax = std::max(fillratio * annz, m*n); // estimated number of nonzeros in U
|
||||
@ -187,8 +185,9 @@ int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size,
|
||||
* \param glu Global data structure
|
||||
* \return 0 on success, > 0 size of the memory allocated so far
|
||||
*/
|
||||
template <typename Scalar, typename Index>
|
||||
template <typename VectorType>
|
||||
int LUMemXpand(VectorType& vec, int& maxlen, int nbElts, LU_MemType memtype, int& num_expansions)
|
||||
int SparseLUBase<Scalar,Index>::LUMemXpand(VectorType& vec, int& maxlen, int nbElts, LU_MemType memtype, int& num_expansions)
|
||||
{
|
||||
int failed_size;
|
||||
if (memtype == USUB)
|
||||
|
||||
@ -15,8 +15,8 @@
|
||||
/**
|
||||
* \brief Count Nonzero elements in the factors
|
||||
*/
|
||||
template <typename IndexVector, typename ScalarVector>
|
||||
void LU_countnz(const int n, int& nnzL, int& nnzU, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
|
||||
template <typename Scalar, typename Index>
|
||||
void SparseLUBase<Scalar,Index>::LU_countnz(const int n, int& nnzL, int& nnzU, GlobalLU_t& glu)
|
||||
{
|
||||
nnzL = 0;
|
||||
nnzU = (glu.xusub)(n);
|
||||
@ -46,8 +46,8 @@ void LU_countnz(const int n, int& nnzL, int& nnzU, LU_GlobalLU_t<IndexVector, Sc
|
||||
* and applies permutation to the remaining subscripts
|
||||
*
|
||||
*/
|
||||
template <typename IndexVector, typename ScalarVector>
|
||||
void LU_fixupL(const int n, const IndexVector& perm_r, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
|
||||
template <typename Scalar, typename Index>
|
||||
void SparseLUBase<Scalar,Index>::LU_fixupL(const int n, const IndexVector& perm_r, GlobalLU_t& glu)
|
||||
{
|
||||
int fsupc, i, j, k, jstart;
|
||||
|
||||
|
||||
@ -46,11 +46,9 @@
|
||||
* > 0 - number of bytes allocated when run out of space
|
||||
*
|
||||
*/
|
||||
template <typename IndexVector, typename ScalarVector, typename BlockIndexVector, typename BlockScalarVector>
|
||||
int LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, ScalarVector& tempv, BlockIndexVector& segrep, BlockIndexVector& repfnz, int fpanelc, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
|
||||
template <typename Scalar, typename Index>
|
||||
int SparseLUBase<Scalar,Index>::LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, ScalarVector& tempv, BlockIndexVector& segrep, BlockIndexVector& repfnz, int fpanelc, GlobalLU_t& glu)
|
||||
{
|
||||
typedef typename IndexVector::Scalar Index;
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
int jsupno, k, ksub, krep, ksupno;
|
||||
int lptr, nrow, isub, irow, nextlu, new_next, ufirst;
|
||||
int fsupc, nsupc, nsupr, luptr, kfnz, no_zeros;
|
||||
@ -95,9 +93,6 @@ int LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, Sca
|
||||
nsupr = glu.xlsub(fsupc+1) - glu.xlsub(fsupc);
|
||||
nrow = nsupr - d_fsupc - nsupc;
|
||||
|
||||
// NOTE Unlike the original implementation in SuperLU, the only feature
|
||||
// available here is a sup-col update.
|
||||
|
||||
// Perform a triangular solver and block update,
|
||||
// then scatter the result of sup-col update to dense
|
||||
no_zeros = kfnz - fst_col;
|
||||
|
||||
@ -42,15 +42,15 @@
|
||||
* \param m number of rows in the matrix
|
||||
* \param jcol Current column
|
||||
* \param perm_r Row permutation
|
||||
* \param maxsuper
|
||||
* \param maxsuper Maximum number of column allowed in a supernode
|
||||
* \param [in,out] nseg Number of segments in current U[*,j] - new segments appended
|
||||
* \param lsub_col defines the rhs vector to start the dfs
|
||||
* \param [in,out] segrep Segment representatives - new segments appended
|
||||
* \param repfnz
|
||||
* \param repfnz First nonzero location in each row
|
||||
* \param xprune
|
||||
* \param marker
|
||||
* \param marker marker[i] == jj, if i was visited during dfs of current column jj;
|
||||
* \param parent
|
||||
* \param xplore
|
||||
* \param xplore working array
|
||||
* \param glu global LU data
|
||||
* \return 0 success
|
||||
* > 0 number of bytes allocated when run out of space
|
||||
@ -60,6 +60,7 @@ 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)
|
||||
{}
|
||||
@ -70,7 +71,7 @@ struct LU_column_dfs_traits
|
||||
void mem_expand(IndexVector& lsub, int& nextl, int chmark)
|
||||
{
|
||||
if (nextl >= m_glu.nzlmax)
|
||||
LUMemXpand<IndexVector>(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions);
|
||||
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 };
|
||||
@ -80,11 +81,9 @@ struct LU_column_dfs_traits
|
||||
LU_GlobalLU_t<IndexVector, ScalarVector>& m_glu;
|
||||
};
|
||||
|
||||
template <typename IndexVector, typename ScalarVector, typename BlockIndexVector>
|
||||
int 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, LU_GlobalLU_t<IndexVector, ScalarVector>& 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)
|
||||
{
|
||||
typedef typename IndexVector::Scalar Index;
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
|
||||
int jsuper = glu.supno(jcol);
|
||||
int nextl = glu.xlsub(jcol);
|
||||
|
||||
@ -43,11 +43,9 @@
|
||||
* > 0 - number of bytes allocated when run out of space
|
||||
*
|
||||
*/
|
||||
template <typename IndexVector, typename ScalarVector, typename SegRepType, typename RepfnzType, typename DenseType>
|
||||
int LU_copy_to_ucol(const int jcol, const int nseg, SegRepType& segrep, RepfnzType& repfnz ,IndexVector& perm_r, DenseType& dense, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
|
||||
{
|
||||
typedef typename IndexVector::Scalar Index;
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
template <typename Scalar, typename Index>
|
||||
int SparseLUBase<Scalar,Index>::LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segrep, BlockIndexVector& repfnz ,IndexVector& perm_r, BlockScalarVector& dense, GlobalLU_t& glu)
|
||||
{
|
||||
Index ksub, krep, ksupno;
|
||||
|
||||
Index jsupno = glu.supno(jcol);
|
||||
|
||||
@ -38,8 +38,8 @@
|
||||
* \param descendants Number of descendants of each node in the etree
|
||||
* \param relax_end last column in a supernode
|
||||
*/
|
||||
template <typename IndexVector>
|
||||
void LU_heap_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end)
|
||||
template <typename Scalar, typename Index>
|
||||
void SparseLUBase<Scalar,Index>::LU_heap_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end)
|
||||
{
|
||||
|
||||
// The etree may not be postordered, but its heap ordered
|
||||
|
||||
@ -30,7 +30,7 @@ template <int SegSizeAtCompileTime> struct LU_kernel_bmod
|
||||
template <typename BlockScalarVector, typename ScalarVector, typename IndexVector>
|
||||
EIGEN_DONT_INLINE static void run(const int segsize, BlockScalarVector& dense, ScalarVector& tempv, ScalarVector& lusup, int& luptr, const int nsupr, const int nrow, IndexVector& lsub, const int lptr, const int no_zeros)
|
||||
{
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
// First, copy U[*,j] segment from dense(*) to tempv(*)
|
||||
// The result of triangular solve is in tempv[*];
|
||||
// The result of matric-vector update is in dense[*]
|
||||
|
||||
@ -34,7 +34,7 @@
|
||||
* \brief Performs numeric block updates (sup-panel) in topological order.
|
||||
*
|
||||
* Before entering this routine, the original nonzeros in the panel
|
||||
* were already copied i nto the spa[m,w] ... FIXME to be checked
|
||||
* were already copied i nto the spa[m,w]
|
||||
*
|
||||
* \param m number of rows in the matrix
|
||||
* \param w Panel size
|
||||
@ -48,10 +48,9 @@
|
||||
*
|
||||
*
|
||||
*/
|
||||
template <typename DenseIndexBlock, typename IndexVector, typename ScalarVector>
|
||||
void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, DenseIndexBlock& segrep, DenseIndexBlock& repfnz, LU_perfvalues& perfv, LU_GlobalLU_t<IndexVector,ScalarVector>& glu)
|
||||
template <typename Scalar, typename Index>
|
||||
void SparseLUBase<Scalar,Index>::LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, IndexVector& repfnz, LU_perfvalues& perfv, GlobalLU_t& glu)
|
||||
{
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
|
||||
int ksub,jj,nextl_col;
|
||||
int fsupc, nsupc, nsupr, nrow;
|
||||
@ -190,9 +189,6 @@ void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, Sca
|
||||
segsize = krep - kfnz + 1;
|
||||
luptr = glu.xlusup(fsupc);
|
||||
|
||||
// NOTE : Unlike the original implementation in SuperLU,
|
||||
// there is no update feature for col-col, 2col-col ...
|
||||
|
||||
// Perform a trianglar solve and block update,
|
||||
// then scatter the result of sup-col update to dense[]
|
||||
no_zeros = kfnz - fsupc;
|
||||
@ -205,4 +201,4 @@ void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, Sca
|
||||
|
||||
} // End for each updating supernode
|
||||
}
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -29,12 +29,12 @@
|
||||
*/
|
||||
#ifndef SPARSELU_PANEL_DFS_H
|
||||
#define SPARSELU_PANEL_DFS_H
|
||||
|
||||
template <typename ScalarVector, typename IndexVector, typename MarkerType, typename Traits>
|
||||
void LU_dfs_kernel(const int jj, IndexVector& perm_r,
|
||||
template <typename Scalar, typename Index>
|
||||
template <typename RepfnzType, typename MarkerType,typename Traits>
|
||||
void SparseLUBase<Scalar,Index>::LU_dfs_kernel(const int jj, IndexVector& perm_r,
|
||||
int& nseg, IndexVector& panel_lsub, IndexVector& segrep,
|
||||
VectorBlock<IndexVector>& repfnz_col, IndexVector& xprune, MarkerType& marker, IndexVector& parent,
|
||||
IndexVector& xplore, LU_GlobalLU_t<IndexVector, ScalarVector>& glu,
|
||||
RepfnzType& repfnz_col, IndexVector& xprune, MarkerType& marker, IndexVector& parent,
|
||||
IndexVector& xplore, GlobalLU_t& glu,
|
||||
int& nextl_col, int krow, Traits& traits
|
||||
)
|
||||
{
|
||||
@ -207,8 +207,8 @@ struct LU_panel_dfs_traits
|
||||
Index* m_marker;
|
||||
};
|
||||
|
||||
template <typename MatrixType, typename ScalarVector, typename IndexVector>
|
||||
void 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, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
|
||||
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)
|
||||
{
|
||||
int nextl_col; // Next available position in panel_lsub[*,jj]
|
||||
|
||||
|
||||
@ -52,12 +52,9 @@
|
||||
* \return 0 if success, i > 0 if U(i,i) is exactly zero
|
||||
*
|
||||
*/
|
||||
template <typename IndexVector, typename ScalarVector>
|
||||
int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, int& pivrow, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
|
||||
template <typename Scalar, typename Index>
|
||||
int SparseLUBase<Scalar,Index>::LU_pivotL(const int jcol, const RealScalar diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, int& pivrow, GlobalLU_t& glu)
|
||||
{
|
||||
typedef typename IndexVector::Scalar Index;
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
typedef typename ScalarVector::RealScalar RealScalar;
|
||||
|
||||
Index fsupc = (glu.xsup)((glu.supno)(jcol)); // First column in the supernode containing the column jcol
|
||||
Index nsupc = jcol - fsupc; // Number of columns in the supernode portion, excluding jcol; nsupc >=0
|
||||
|
||||
@ -46,12 +46,9 @@
|
||||
* \param glu Global LU data
|
||||
*
|
||||
*/
|
||||
template <typename IndexVector, typename ScalarVector, typename BlockIndexVector>
|
||||
void LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, const IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
|
||||
template <typename Scalar, typename Index>
|
||||
void SparseLUBase<Scalar,Index>::LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, const IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, GlobalLU_t& glu)
|
||||
{
|
||||
typedef typename IndexVector::Scalar Index;
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
|
||||
// For each supernode-rep irep in U(*,j]
|
||||
int jsupno = glu.supno(jcol);
|
||||
int i,irep,irep1;
|
||||
|
||||
@ -37,8 +37,8 @@
|
||||
* \param descendants Number of descendants of each node in the etree
|
||||
* \param relax_end last column in a supernode
|
||||
*/
|
||||
template <typename IndexVector>
|
||||
void LU_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end)
|
||||
template <typename Scalar, typename Index>
|
||||
void SparseLUBase<Scalar,Index>::LU_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end)
|
||||
{
|
||||
|
||||
// compute the number of descendants of each node in the etree
|
||||
|
||||
@ -29,11 +29,9 @@
|
||||
*/
|
||||
#ifndef SPARSELU_SNODE_BMOD_H
|
||||
#define SPARSELU_SNODE_BMOD_H
|
||||
template <typename IndexVector, typename ScalarVector>
|
||||
int LU_snode_bmod (const int jcol, const int fsupc, ScalarVector& dense, LU_GlobalLU_t<IndexVector,ScalarVector>& glu)
|
||||
{
|
||||
typedef typename ScalarVector::Scalar Scalar;
|
||||
|
||||
template <typename Scalar, typename Index>
|
||||
int SparseLUBase<Scalar,Index>::LU_snode_bmod (const int jcol, const int fsupc, ScalarVector& dense, GlobalLU_t& glu)
|
||||
{
|
||||
/* lsub : Compressed row subscripts of ( rectangular supernodes )
|
||||
* xlsub : xlsub[j] is the starting location of the j-th column in lsub(*)
|
||||
* lusup : Numerical values of the rectangular supernodes
|
||||
|
||||
@ -42,55 +42,54 @@
|
||||
* \param marker (in/out) working vector
|
||||
* \return 0 on success, > 0 size of the memory when memory allocation failed
|
||||
*/
|
||||
template <typename MatrixType, typename IndexVector, typename ScalarVector>
|
||||
int LU_snode_dfs(const int jcol, const int kcol,const MatrixType& mat, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
|
||||
template <typename Scalar, typename Index>
|
||||
int SparseLUBase<Scalar,Index>::LU_snode_dfs(const int jcol, const int kcol,const MatrixType& mat, IndexVector& xprune, IndexVector& marker, GlobalLU_t& glu)
|
||||
{
|
||||
int mem;
|
||||
Index nsuper = ++glu.supno(jcol); // Next available supernode number
|
||||
int nextl = glu.xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
|
||||
int krow,kmark;
|
||||
for (int i = jcol; i <=kcol; i++)
|
||||
{
|
||||
typedef typename IndexVector::Scalar Index;
|
||||
int mem;
|
||||
Index nsuper = ++glu.supno(jcol); // Next available supernode number
|
||||
int nextl = glu.xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
|
||||
int krow,kmark;
|
||||
for (int i = jcol; i <=kcol; i++)
|
||||
// For each nonzero in A(*,i)
|
||||
for (typename MatrixType::InnerIterator it(mat, i); it; ++it)
|
||||
{
|
||||
// For each nonzero in A(*,i)
|
||||
for (typename MatrixType::InnerIterator it(mat, i); it; ++it)
|
||||
krow = it.row();
|
||||
kmark = marker(krow);
|
||||
if ( kmark != kcol )
|
||||
{
|
||||
krow = it.row();
|
||||
kmark = marker(krow);
|
||||
if ( kmark != kcol )
|
||||
// First time to visit krow
|
||||
marker(krow) = kcol;
|
||||
glu.lsub(nextl++) = krow;
|
||||
if( nextl >= glu.nzlmax )
|
||||
{
|
||||
// First time to visit krow
|
||||
marker(krow) = kcol;
|
||||
glu.lsub(nextl++) = krow;
|
||||
if( nextl >= glu.nzlmax )
|
||||
{
|
||||
mem = LUMemXpand<IndexVector>(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions);
|
||||
if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
|
||||
}
|
||||
mem = LUMemXpand<IndexVector>(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions);
|
||||
if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
|
||||
}
|
||||
}
|
||||
glu.supno(i) = nsuper;
|
||||
}
|
||||
|
||||
// If supernode > 1, then make a copy of the subscripts for pruning
|
||||
if (jcol < kcol)
|
||||
{
|
||||
Index new_next = nextl + (nextl - glu.xlsub(jcol));
|
||||
while (new_next > glu.nzlmax)
|
||||
{
|
||||
mem = LUMemXpand<IndexVector>(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions);
|
||||
if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
|
||||
}
|
||||
Index ifrom, ito = nextl;
|
||||
for (ifrom = glu.xlsub(jcol); ifrom < nextl;)
|
||||
glu.lsub(ito++) = glu.lsub(ifrom++);
|
||||
for (int i = jcol+1; i <=kcol; i++) glu.xlsub(i) = nextl;
|
||||
nextl = ito;
|
||||
}
|
||||
glu.xsup(nsuper+1) = kcol + 1; // Start of next available supernode
|
||||
glu.supno(kcol+1) = nsuper;
|
||||
xprune(kcol) = nextl;
|
||||
glu.xlsub(kcol+1) = nextl;
|
||||
return 0;
|
||||
glu.supno(i) = nsuper;
|
||||
}
|
||||
|
||||
// If supernode > 1, then make a copy of the subscripts for pruning
|
||||
if (jcol < kcol)
|
||||
{
|
||||
Index new_next = nextl + (nextl - glu.xlsub(jcol));
|
||||
while (new_next > glu.nzlmax)
|
||||
{
|
||||
mem = LUMemXpand<IndexVector>(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions);
|
||||
if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
|
||||
}
|
||||
Index ifrom, ito = nextl;
|
||||
for (ifrom = glu.xlsub(jcol); ifrom < nextl;)
|
||||
glu.lsub(ito++) = glu.lsub(ifrom++);
|
||||
for (int i = jcol+1; i <=kcol; i++) glu.xlsub(i) = nextl;
|
||||
nextl = ito;
|
||||
}
|
||||
glu.xsup(nsuper+1) = kcol + 1; // Start of next available supernode
|
||||
glu.supno(kcol+1) = nsuper;
|
||||
xprune(kcol) = nextl;
|
||||
glu.xlsub(kcol+1) = nextl;
|
||||
return 0;
|
||||
}
|
||||
#endif
|
||||
@ -118,7 +118,7 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
|
||||
{
|
||||
eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
|
||||
|
||||
// FIXME Stability: We should probably compute the scaling factors and the shifts that are needed to ensure an efficient LLT preconditioner.
|
||||
// FIXME Stability: We should probably compute the scaling factors and the shifts that are needed to ensure a succesful LLT factorization and an efficient preconditioner.
|
||||
|
||||
// Dropping strategies : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added
|
||||
|
||||
@ -177,8 +177,8 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
|
||||
// p is the original number of elements in the column (without the diagonal)
|
||||
int p = colPtr[j+1] - colPtr[j] - 2 ;
|
||||
internal::QuickSplit(curCol, irow, p);
|
||||
if(RealScalar(diag) <= 0)
|
||||
{
|
||||
if(RealScalar(diag) <= 0)
|
||||
{ //FIXME We can use heuristics (Kershaw, 1978 or above reference ) to get a dynamic shift
|
||||
m_info = NumericalIssue;
|
||||
return;
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user