mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-30 17:40:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add preliminary files for SparseLU

Browse Source
This commit is contained in:
Desire NUENTSA 2012-05-25 18:17:57 +02:00
parent b202c5ed2f
commit b6267507ea
12 changed files with 1781 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

341
Eigen/src/SparseLU/SparseLU.h Normal file
View File

@ -0,0 +1,341 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_SPARSE_LU
#define EIGEN_SPARSE_LU
#include <Ordering.h>
#include <SparseLU_Utils.h>
#include <SuperNodalMatrix.h>
#include <SparseLU_Structs.h>
#include <SparseLU_Memory.h>
#include <SparseLU_Coletree.h>
namespace Eigen {
template <typename _MatrixType>
class SparseLU
{
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index;
typedef SparseMatrix<Scalar,ColMajor,Index> NCMatrix;
typedef SuperNodalMatrix<Scalar, Index> SCMatrix;
typedef GlobalLU_t<Scalar, Index> Eigen_GlobalLU_t;
typedef Matrix<Scalar,Dynamic,1> VectorType;
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
public:
SparseLU():m_isInitialized(true),m_symmetricmode(false),m_fact(DOFACT),m_diagpivotthresh(1.0)
{
initperfvalues();
}
SparseLU(const MatrixType& matrix):SparseLU()
{
compute(matrix);
}
~SparseLU()
{
}
void analyzePattern (const MatrixType& matrix);
void factorize (const MatrixType& matrix);
void compute (const MatrixType& matrix);
/** Indicate that the pattern of the input matrix is symmetric */
void isSymmetric(bool sym)
{
m_symmetricmode = sym;
}
/** Set the threshold used for a diagonal entry to be an acceptable pivot. */
void diagPivotThresh(RealScalar thresh)
{
m_diagpivotthresh = thresh;
}
protected:
// Functions
void initperfvalues();
// Variables
mutable ComputationInfo m_info;
bool m_isInitialized;
bool m_factorizationIsOk;
bool m_analysisIsOk;
fact_t m_fact;
NCMatrix m_mat; // The input (permuted ) matrix
SCMatrix m_Lstore; // The lower triangular matrix (supernodal)
NCMatrix m_Ustore; //The upper triangular matrix
PermutationType m_perm_c; // Column permutation
PermutationType m_iperm_c; // Column permutation
PermutationType m_perm_r ; // Row permutation
PermutationType m_iperm_r ; // Inverse row permutation
VectorXi m_etree; // Column elimination tree
Scalar *m_work; //
Index *m_iwork; //
static Eigen_GlobalLU_t m_Glu; // persistent data to facilitate multiple factors
// should be defined as a class member
// SuperLU/SparseLU options
bool m_symmetricmode;
// values for performance
int m_panel_size; // a panel consists of at most <panel_size> consecutive columns
int m_relax; // To control degree of relaxing supernodes. If the number of nodes (columns)
// in a subtree of the elimination tree is less than relax, this subtree is considered
// as one supernode regardless of the row structures of those columns
int m_maxsuper; // The maximum size for a supernode in complete LU
int m_rowblk; // The minimum row dimension for 2-D blocking to be used;
int m_colblk; // The minimum column dimension for 2-D blocking to be used;
int m_fillfactor; // The estimated fills factors for L and U, compared with A
RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot
private:
// Copy constructor
SparseLU (SparseLU& ) {}
}; // End class SparseLU
/* Set the default values for performance */
void SparseLU::initperfvalues()
{
m_panel_size = 12;
m_relax = 1;
m_maxsuper = 100;
m_rowblk = 200;
m_colblk = 60;
m_fillfactor = 20;
}
/**
* Compute the column permutation to minimize the fill-in (file amd.c )
* - Apply this permutation to the input matrix -
* - Compute the column elimination tree on the permuted matrix (file Eigen_Coletree.h)
* - Postorder the elimination tree and the column permutation (file Eigen_Coletree.h)
* -
*/
template <typename MatrixType>
void SparseLU::analyzePattern(const MatrixType& mat)
{
// Compute the column permutation
AMDordering amd(mat);
m_perm_c = amd.get_perm_c();
// Apply the permutation to the column of the input matrix
m_mat = mat * m_perm_c; //how is the permutation represented ???
// Compute the column elimination tree of the permuted matrix
if (m_etree.size() == 0) m_etree.resize(m_mat.cols());
internal::sp_coletree(m_mat, m_etree);
// In symmetric mode, do not do postorder here
if (m_symmetricmode == false) {
VectorXi post, iwork;
// Post order etree
post = internal::TreePostorder(m_mat.cols(), m_etree);
// Renumber etree in postorder
iwork.resize(n+1);
for (i = 0; i < n; ++i) iwork(post(i)) = post(m_etree(i));
m_etree = iwork;
// Postmultiply A*Pc by post,
// i.e reorder the matrix according to the postorder of the etree
// FIXME Check if this is available : constructor from a vector
PermutationType post_perm(post);
m_mat = m_mat * post_perm;
// Product of m_perm_c and post
for (i = 0; i < n; ++i) iwork(i) = m_perm_c(post_perm.indices()(i));
m_perm_c = iwork;
} // end postordering
}
/**
* - Numerical factorization
* - Interleaved with the symbolic factorization
* \tparam MatrixType The type of the matrix, it should be a column-major sparse matrix
* \return info where
* : successful exit
* = 0: successful exit
* > 0: if info = i, and i is
* <= 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.
* > 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.
*/
template <typename MatrixType>
int SparseLU::factorize(const MatrixType& matrix)
{
// Allocate storage common to the factor routines
int lwork = 0;
int info = LUMemInit(lwork);
eigen_assert ( (info == 0) && "Unable to allocate memory for the factors");
int m = m_mat.rows();
int n = m_mat.cols();
int maxpanel = m_panel_size * m;
// Set up pointers for integer working arrays
Map<VectorXi> segrep(m_iwork, m); //
Map<VectorXi> parent(&segrep(0) + m, m); //
Map<VectorXi> xplore(&parent(0) + m, m); //
Map<VectorXi> repfnz(&xplore(0) + m, maxpanel); //
Map<VectorXi> panel_lsub(&repfnz(0) + maxpanel, maxpanel);//
Map<VectorXi> xprune(&panel_lsub(0) + maxpanel, n); //
Map<VectorXi> marker(&xprune(0)+n, m * LU_NO_MARKER); //
repfnz.setConstant(-1);
panel_lsub.setConstant(-1);
// Set up pointers for scalar working arrays
VectorType dense(maxpanel);
dense.setZero();
VectorType tempv(LU_NUM_TEMPV(m,m_panel_size,m_maxsuper,m_rowblk);
tempv.setZero();
// Setup Permutation vectors
PermutationType iperm_r; // inverse of perm_r
if (m_fact = SamePattern_SameRowPerm)
iperm_r = m_perm_r.inverse();
// Compute the inverse of perm_c
PermutationType iperm_c;
iperm_c = m_perm_c.inverse();
// Identify initial relaxed snodes
VectorXi relax_end(n);
if ( m_symmetricmode = true )
LU_heap_relax_snode(n, m_etree, m_relax, marker, relax_end);
else
LU_relax_snode(n, m_etree, m_relax, marker, relax_end);
m_perm_r.setConstant(-1);
marker.setConstant(-1);
VectorXi& xsup = m_Glu.xsup;
VectorXi& supno = m_GLu.supno;
VectorXi& xlsub = m_Glu.xlsub;
VectorXi& xlusup = m_GLu.xlusup;
VectorXi& xusub = m_Glu.xusub;
supno(0) = -1;
xsup(0) = xlsub(0) = xusub(0) = xlusup(0);
int panel_size = m_panel_size;
int wdef = panel_size; // upper bound on panel width
// 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
// (b) panel_size contiguous columns, <panel_size> defined by the user
register int jcol,kcol;
int min_mn = std::min(m,n);
VectorXi panel_histo(n);
bool ok = true;
Index nextu, nextlu, jsupno, fsupc, new_next;
int pivrow; // Pivotal row number in the original row matrix
int nseg1; // Number of segments in U-column above panel row jcol
int nseg; // Number of segments in each U-column
for (jcol = 0; jcol < min_mn; )
{
if (relax_end(jcol) != -1)
{ // Starting a relaxed node from jcol
kcol = relax_end(jcol); // End index of the relaxed snode
// 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.innerIndexPtr(), m_mat.outerIndexPtr(), xprune, marker);
if ( !info )
{
ok = false;
break;
}
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);
nzlumax = m_Glu.nzlumax;
while (new_next > nzlumax )
{
m_Glu.lusup = LUMemXpand<Scalar>(jcol, nextlu, LUSUP, nzlumax);
m_GLu.nzlumax = nzlumax;
}
// Now, left-looking factorize each column within the snode
for (icol = jcol; icol<=kcol; icol++){
xusub(icol+1) = nextu;
// Scatter into SPA dense(*)
for (typename MatrixType::InnerIterator it(m_mat, icol); it; ++it)
dense(it.row()) = it.val();
// Numeric update within the snode
LU_snode_bmod(icol, jsupno, fsupc, dense, tempv);
// Eliminate the current column
info = LU_pivotL(icol, pivrow);
eigen_assert(info == 0 && "The matrix is structurally singular");
}
jcol = icol; // The last column te be eliminated
}
else
{ // Work on one panel of panel_size columns
// Adjust panel size so that a panel won't overlap with the next relaxed snode.
panel_size = w_def;
for (k = jcol + 1; k < std::min(jcol+panel_size, min_mn); k++)
{
if (relax_end(k) != -1)
{
panel_size = k - jcol;
break;
}
}
if (k == min_mn)
panel_size = min_mn - jcol;
// Symbolic outer factorization on a panel of columns
LU_panel_dfs(m, panel_size, jcol, m_mat, m_perm_r, 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);
// Sparse LU within the panel, and below the panel diagonal
for ( jj = jcol, j< jcol + panel_size; jj++)
{
k = (jj - jcol) * m; // Column index for w-wide arrays
} // end for
jcol += panel_size; // Move to the next panel
} // end else
} // end for -- end elimination
m_info = ok ? Success : NumericalIssue;
m_factorizationIsOk = ok;
}
} // End namespace Eigen
#endif

188
Eigen/src/SparseLU/SparseLU_Coletree.h Normal file
View File

@ -0,0 +1,188 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
/*
* NOTE: This file is the modified version of sp_coletree.c file in SuperLU
* -- SuperLU routine (version 3.1) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* August 1, 2008
*
* Copyright (c) 1994 by Xerox Corporation. All rights reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* Permission to modify the code and to distribute modified code is
* granted, provided the above notices are retained, and a notice that
* the code was modified is included with the above copyright notice.
*/
#ifndef SPARSELU_COLETREE_H
#define SPARSELU_COLETREE_H
/** Compute the column elimination tree of a sparse matrix
* NOTE : The matrix is supposed to be in column-major format.
*
*/
template<typename MatrixType>
int LU_sp_coletree(const MatrixType& mat, VectorXi& parent)
{
int nc = mat.cols(); // Number of columns
int nr = mat.rows(); // Number of rows
VectorXi root(nc); // root of subtree of etree
root.setZero();
VectorXi pp(nc); // disjoint sets
pp.setZero(); // Initialize disjoint sets
VectorXi firstcol(nr); // First nonzero column in each row
firstcol.setZero();
//Compute firstcol[row]
int row,col;
firstcol.setConstant(nc); //for (row = 0; row < nr; firstcol(row++) = nc);
for (col = 0; col < nc; col++)
{
for (typename MatrixType::InnerIterator it(mat, col); it; ++it)
{ // Is it necessary to brows the whole matrix, the lower part should do the job ??
row = it.row();
firstcol(row) = std::min(firstcol(row), col);
}
}
/* Compute etree by Liu's algorithm for symmetric matrices,
except use (firstcol[r],c) in place of an edge (r,c) of A.
Thus each row clique in A'*A is replaced by a star
centered at its first vertex, which has the same fill. */
int rset, cset, rroot;
for (col = 0; col < nc; col++)
{
pp(col) = cset = col; // Initially, each element is in its own set
root(cset) = col;
parent(col) = nc;
for (typename MatrixType::InnerIterator it(mat, col); it; ++it)
{ // A sequence of interleaved find and union is performed
row = firstcol(it.row());
if (row >= col) continue;
rset = internal::etree_find(row, pp); // Find the name of the set containing row
rroot = root(rset);
if (rroot != col)
{
parent(rroot) = col;
pp(cset) = cset = rset; // Get the union of cset and rset
root(cset) = col;
}
}
}
return 0;
}
/** Find the root of the tree/set containing the vertex i : Use Path halving */
int etree_find (int i, VectorXi& pp)
{
int p = pp(i); // Parent
int gp = pp(p); // Grand parent
while (gp != p)
{
pp(i) = gp; // Parent pointer on find path is changed to former grand parent
i = gp;
p = pp(i);
gp = pp(p);
}
return p;
}
/**
* Post order a tree
*/
VectorXi TreePostorder(int n, VectorXi& parent)
{
VectorXi first_kid, next_kid; // Linked list of children
VectorXi post; // postordered etree
int postnum;
// Allocate storage for working arrays and results
first_kid.resize(n+1);
next_kid.setZero(n+1);
post.setZero(n+1);
// Set up structure describing children
int v, dad;
first_kid.setConstant(-1);
for (v = n-1, v >= 0; v--)
{
dad = parent(v);
next_kid(v) = first_kid(dad);
first_kid(dad) = v;
}
// Depth-first search from dummy root vertex #n
postnum = 0;
internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
return post;
}
/**
* Depth-first search from vertex n. No recursion.
* This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
*/
void nr_etdfs (int n, int *parent, int* first_kid, int *next_kid, int *post, int postnum)
{
int current = n, first, next;
while (postnum != n)
{
// No kid for the current node
first = first_kid(current);
// no first kid for the current node
if (first == -1)
{
// Numbering this node because it has no kid
post(current) = postnum++;
// looking for the next kid
next = next_kid(current);
while (next == -1)
{
// No more kids : back to the parent node
current = parent(current);
// numbering the parent node
post(current) = postnum++;
// Get the next kid
next = next_kid(current);
}
// stopping criterion
if (postnum==n+1) return;
// Updating current node
current = next;
}
else
{
current = first;
}
}
}
#endif

74
Eigen/src/SparseLU/SparseLU_Matrix.h Normal file
View File

@ -0,0 +1,74 @@
// 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>
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_SPARSELU_MATRIX_H
#define EIGEN_SPARSELU_MATRIX_H
/** \ingroup SparseLU_Module
* \brief a class to manipulate the supernodal matrices in the SparseLU factorization
*
* This class extends the class SparseMatrix and should contain the data to easily store
* and manipulate the supernodes during the factorization and solution phase of Sparse LU.
* Only the lower triangular matrix has supernodes.
*
* NOTE : This class corresponds to the SCformat structure in SuperLU
*
*/
template <typename _Scalar, typename _Index>
class SuperNodalMatrix
{
public:
SCMatrix()
{
}
~SCMatrix()
{
}
operator SparseMatrix();
protected:
Index nnz; // Number of nonzero values
Index nsupper; // Index of the last supernode
Scalar *nzval; //array of nonzero values packed by (supernode ??) column
Index *nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j
Index *rowind; // Array of compressed row indices of rectangular supernodes
Index rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j
Index *col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs
Index *sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode
// Index *nzval_colptr corresponds to m_outerIndex in SparseMatrix
private :
SuperNodalMatrix(SparseMatrix& ) {}
};
SuperNodalMatrix::operator SparseMatrix()
{
}
#endif

242
Eigen/src/SparseLU/SparseLU_Memory.h Normal file
View File

@ -0,0 +1,242 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
/*
* NOTE: This file is the modified version of [s,d,c,z]memory.c files in SuperLU
* -- SuperLU routine (version 3.1) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* August 1, 2008
*
* Copyright (c) 1994 by Xerox Corporation. All rights reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* Permission to modify the code and to distribute modified code is
* granted, provided the above notices are retained, and a notice that
* the code was modified is included with the above copyright notice.
*/
#ifndef EIGEN_SPARSELU_MEMORY
#define EIGEN_SPARSELU_MEMORY
#define LU_Reduce(alpha) ((alpha + 1) / 2) // i.e (alpha-1)/2 + 1
#define LU_GluIntArray(n) (5* (n) + 5)
#define LU_TempSpace(m, w) ( (2*w + 4 + LU_NO_MARKER) * m * sizeof(Index) \
+ (w + 1) * m * sizeof(Scalar)
namespace internal {
/* Allocate various working space needed in the numerical factorization phase.
* m_work : space fot the output data structures (lwork is the size)
* m_Glu: persistent data to facilitate multiple factors : is it really necessary ??
* NOTE Unlike SuperLU, this routine does not allow the user to provide the size to allocate
* nor it return an estimated amount of space required.
*
* Useful variables :
* - m_fillratio : Ratio of fill expected
* - lwork = -1 : return an estimated size of the required memory
* = 0 : Estimate and allocate the memory
*/
template <typename Scalar,typename Index>
int SparseLU::LUMemInit(int lwork)
{
int iword = sizeof(Index);
int dword = sizeof(Scalar);
int n = m_Glu.n = m_mat.cols();
int m = m_mat.rows();
m_Glu.num_expansions = 0; // No memory expansions so far ??
int estimated_size;
if (!m_Glu.expanders)
m_Glu.expanders = new ExpHeader(NO_MEMTYPE);
if (m_fact_t != SamePattern_SameRowPerm) // Create space for a new factorization
{
// Guess the size for L\U factors
int annz = m_mat.nonZeros();
int nzlmax, nzumax, nzlumax;
nzumax = nzlumax = m_fillratio * annz; // ???
nzlmax = std::max(1, m_fill_ratio/4.) * annz; //???
// Return the estimated size to the user if necessary
if (lwork = -1)
{
estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size)
+ (nzlmax + nzumax) * iword + (nzlumax+nzumax) * dword + n);
return estimated_size;
}
// Setup the required space
// NOTE: In SuperLU, there is an option to let the user provide its own space.
// Allocate Integer pointers for L\U factors.resize(n+1);
m_Glu.supno.resize(n+1);
m_Glu.xlsub.resize(n+1);
m_Glu.xlusup.resize(n+1);
m_Glu.xusub.resize(n+1);
// Reserve memory for L/U factors
m_Glu.lusup = internal::expand<Scalar>(nzlumax, LUSUP, 0, 0, m_Glu);
m_Glu.ucol = internal::expand<Scalar>(nzumax, UCOL, 0, 0, m_Glu);
m_Glu.lsub = internal::expand<Index>(nzlmax, LSUB, 0, 0, m_Glu);
m_Glu.usub = internal::expand<Index>(nzumax, USUB, 0, 1, m_Glu);
// Check if the memory is correctly allocated,
while ( !m_Glu.lusup || !m_Glu.ucol || !m_Glu.lsub || !m_Glu.usub)
{
//otherwise reduce the estimated size and retry
delete [] m_Glu.lusup;
delete [] m_Glu.ucol;
delete [] m_Glu.lsub;
delete [] m_Glu.usub;
nzlumax /= 2;
nzumax /= 2;
nzlmax /= 2;
eigen_assert (nzlumax > annz && "Not enough memory to perform factorization");
m_Glu.lusup = internal::expand<Scalar>(nzlumax, LUSUP, 0, 0, m_Glu);
m_Glu.ucol = internal::expand<Scalar>(nzumax, UCOL, 0, 0, m_Glu);
m_Glu.lsub = internal::expand<Index>(nzlmax, LSUB, 0, 0, m_Glu);
m_Glu.usub = internal::expand<Index>(nzumax, USUB, 0, 1, m_Glu);
}
}
else // m_fact == SamePattern_SameRowPerm;
{
if (lwork = -1)
{
estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size)
+ (Glu.nzlmax + Glu.nzumax) * iword + (Glu.nzlumax+Glu.nzumax) * dword + n);
return estimated_size;
}
// Use existing space from previous factorization
// Unlike in SuperLU, this should not be necessary here since m_Glu is persistent as a member of the class
m_Glu.xsup = m_Lstore.sup_to_col;
m_Glu.supno = m_Lstore.col_to_sup;
m_Glu.xlsub = m_Lstore.rowind_colptr;
m_Glu.xlusup = m_Lstore.nzval_colptr;
xusub = m_Ustore.outerIndexPtr();
m_Glu.expanders[LSUB].size = m_Glu.nzlmax; // Maximum value from previous factorization
m_Glu.expanders[LUSUP].size = m_Glu.nzlumax;
m_Glu.expanders[USUB].size = GLu.nzumax;
m_Glu.expanders[UCOL].size = m_Glu.nzumax;
m_Glu.lsub = GLu.expanders[LSUB].mem = m_Lstore.rowind;
m_Glu.lusup = GLu.expanders[LUSUP].mem = m_Lstore.nzval;
GLu.usub = m_Glu.expanders[USUB].mem = m_Ustore.InnerIndexPtr();
m_Glu.ucol = m_Glu.expanders[UCOL].mem = m_Ustore.valuePtr();
}
// LUWorkInit : Now, allocate known working storage
int isize = (2 * m_panel_size + 3 + LU_NO_MARKER) * m + n;
int dsize = m * m_panel_size + LU_NUM_TEMPV(m, m_panel_size, m_maxsuper, m_rowblk);
m_iwork = new Index(isize);
eigen_assert( (m_iwork != 0) && "Malloc fails for iwork");
m_work = new Scalar(dsize);
eigen_assert( (m_work != 0) && "Malloc fails for dwork");
++m_Glu.num_expansions;
return 0;
} // end LuMemInit
/**
* Expand the existing storage to accomodate more fill-ins
*/
template <typename DestType >
DestType* SparseLU::expand(int& prev_len, // Length from previous call
MemType type, // Which part of the memory to expand
int len_to_copy, // Size of the memory to be copied to new store
int keep_prev) // = 1: use prev_len; Do not expand this vector
// = 0: compute new_len to expand)
{
float alpha = 1.5; // Ratio of the memory increase
int new_len; // New size of the allocated memory
if(m_Glu.num_expansions == 0 || keep_prev)
new_len = prev_len;
else
new_len = alpha * prev_len;
// Allocate new space
DestType *new_mem, *old_mem;
new_mem = new DestType(new_len);
if ( m_Glu.num_expansions != 0 ) // The memory has been expanded before
{
int tries = 0;
if (keep_prev)
{
if (!new_mem) return 0;
}
else
{
while ( !new_mem)
{
// Reduce the size and allocate again
if ( ++tries > 10) return 0;
alpha = LU_Reduce(alpha);
new_len = alpha * prev_len;
new_mem = new DestType(new_len);
}
} // keep_prev
//Copy the previous values to the newly allocated space
ExpHeader<DestType>* expanders = m_Glu.expanders;
std::memcpy(new_mem, expanders[type].mem, len_to_copy);
delete [] expanders[type].mem;
}
expanders[type].mem = new_mem;
expanders[type].size = new_len;
prev_len = new_len;
if(m_Glu.num_expansions) ++m_Glu.num_expansions;
return expanders[type].mem;
}
/**
* \brief Expand the existing storage
*
* NOTE: The calling sequence of this function is different from that of SuperLU
*
* \return a pointer to the newly allocated space
*/
template <typename DestType>
DestType* SparseLU::LUMemXpand(int jcol, int next, MemType mem_type, int& maxlen)
{
DestType *newmem;
if (memtype == USUB)
new_mem = expand<DestType>(maxlen, mem_type, next, 1);
else
new_mem = expand<DestType>(maxlen, mem_type, next, 0);
eigen_assert(new_mem && "Can't expand memory");
return new_mem;
}
}// Namespace Internal
#endif

122
Eigen/src/SparseLU/SparseLU_Structs.h Normal file
View File

@ -0,0 +1,122 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
/*
* NOTE: Part of this file is the modified version of files slu_[s,d,c,z]defs.h
* -- SuperLU routine (version 4.1) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* November, 2010
*
* Global data structures used in LU factorization -
*
* nsuper: #supernodes = nsuper + 1, numbered [0, nsuper].
* (xsup,supno): supno[i] is the supernode no to which i belongs;
* xsup(s) points to the beginning of the s-th supernode.
* e.g. supno 0 1 2 2 3 3 3 4 4 4 4 4 (n=12)
* xsup 0 1 2 4 7 12
* Note: dfs will be performed on supernode rep. relative to the new
* row pivoting ordering
*
* (xlsub,lsub): lsub[*] contains the compressed subscript of
* rectangular supernodes; xlsub[j] points to the starting
* location of the j-th column in lsub[*]. Note that xlsub
* is indexed by column.
* Storage: original row subscripts
*
* During the course of sparse LU factorization, we also use
* (xlsub,lsub) for the purpose of symmetric pruning. For each
* supernode {s,s+1,...,t=s+r} with first column s and last
* column t, the subscript set
* lsub[j], j=xlsub[s], .., xlsub[s+1]-1
* is the structure of column s (i.e. structure of this supernode).
* It is used for the storage of numerical values.
* Furthermore,
* lsub[j], j=xlsub[t], .., xlsub[t+1]-1
* is the structure of the last column t of this supernode.
* It is for the purpose of symmetric pruning. Therefore, the
* structural subscripts can be rearranged without making physical
* interchanges among the numerical values.
*
* However, if the supernode has only one column, then we
* only keep one set of subscripts. For any subscript interchange
* performed, similar interchange must be done on the numerical
* values.
*
* The last column structures (for pruning) will be removed
* after the numercial LU factorization phase.
*
* (xlusup,lusup): lusup[*] contains the numerical values of the
* rectangular supernodes; xlusup[j] points to the starting
* location of the j-th column in storage vector lusup[*]
* Note: xlusup is indexed by column.
* Each rectangular supernode is stored by column-major
* scheme, consistent with Fortran 2-dim array storage.
*
* (xusub,ucol,usub): ucol[*] stores the numerical values of
* U-columns outside the rectangular supernodes. The row
* subscript of nonzero ucol[k] is stored in usub[k].
* 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
namespace Eigen {
#define NO_MEMTYPE 4 /* 0: lusup
1: ucol
2: lsub
3: usub */
typedef enum {NATURAL, MMD_ATA, MMD_AT_PLUS_A, COLAMD, MY_PREMC} colperm_t;
typedef enum {DOFACT, SamePattern, SamePattern_SameRowPerm, Factored} fact_t;
typedef enum {LUSUP, UCOL, LSUB, USUB, LLVL, ULVL} MemType;
/** Headers for dynamically managed memory
\tparam BaseType can be int, real scalar or complex scalar*/
template <typename BaseType>
struct ExpHeader {
int size; // Length of the memory that has been used */
BaseType *mem;
} ExpHeader;
template <typename VectorType, typename Index>
struct {
VectorXi xsup; // supernode and column mapping
VectorXi supno; // Supernode number corresponding to this column
VectorXi lsub; // Compressed L subscripts of rectangular supernodes
VectorXi xlsub; // xlsub(j) points to the starting location of the j-th column in lsub
VectorXi xlusup;
VectorXi xusub;
VectorType lusup; // L supernodes
VectorType ucol; // U columns
Index nzlmax; // Current max size of lsub
Index nzumax; // Current max size of ucol
Index nzlumax; // Current max size of lusup
Index n; // Number of columns in the matrix
int num_expansions;
ExpHeader *expanders; // Array of pointers to 4 types of memory
} GlobalLU_t;
}// End namespace Eigen
#endif

32
Eigen/src/SparseLU/SparseLU_Utils.h Normal file
View File

@ -0,0 +1,32 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifdef EIGEN_SPARSELU_UTILS_H
#define EIGEN_SPARSELU_UTILS_H
// Number of marker arrays used in the symbolic factorization each of size n
#define LU_NO_MARKER 3
#define LU_NUM_TEMPV(m,w,t,b) (std::max(m, (t+b)*w) )
#define LU_EMPTY (-1)
#endif

133
Eigen/src/SparseLU/SparseLU_heap_relax_snode.h Normal file
View File

@ -0,0 +1,133 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
/* This file is a modified version of heap_relax_snode.c file in SuperLU
* -- SuperLU routine (version 3.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* October 15, 2003
*
* Copyright (c) 1994 by Xerox Corporation. All rights reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* Permission to modify the code and to distribute modified code is
* granted, provided the above notices are retained, and a notice that
* the code was modified is included with the above copyright notice.
*/
#ifndef EIGEN_HEAP_RELAX_SNODE_H
#define EIGEN_HEAP_RELAX_SNODE_H
#include <Eigen_coletree.h>
/**
* \brief Identify the initial relaxed supernodes
*
* This routine applied to a symmetric elimination tree.
* It assumes that the matrix has been reordered according to the postorder of the etree
* \param et elimination tree
* \param relax_columns Maximum number of columns allowed in a relaxed snode
* \param descendants Number of descendants of each node in the etree
* \param relax_end last column in a supernode
*/
void internal::LU_heap_relax_snode (const int n, VectorXi& et, const int relax_columns, VectorXi& descendants, VectorXi& relax_end)
{
// The etree may not be postordered, but its heap ordered
// Post order etree
VectorXi post = internal::TreePostorder(n, et);
VectorXi inv_post(n+1);
register int i;
for (i = 0; i < n+1; ++i) inv_post(post(i)) = i;
// Renumber etree in postorder
VectorXi iwork(n);
VectorXi et_save(n+1);
for (i = 0; i < n; ++i)
{
iwork(post(i)) = post(et(i));
}
et_save = et; // Save the original etree
et = iwork;
// compute the number of descendants of each node in the etree
relax_end.setConstant(-1);
register int j, parent;
descendants.setZero();
for (j = 0; j < n; j++)
{
parent = et(j);
if (parent != n) // not the dummy root
descendants(parent) += descendants(j) + 1;
}
// Identify the relaxed supernodes by postorder traversal of the etree
register int snode_start; // beginning of a snode
register int k;
int nsuper_et_post = 0; // Number of relaxed snodes in postordered etree
int nsuper_et = 0; // Number of relaxed snodes in the original etree
for (j = 0; j < n; )
{
parent = et(j);
snode_start = j;
while ( parent != n && descendants(parent) < relax_columns )
{
j = parent;
parent = et(j);
}
// Found a supernode in postordered etree, j is the last column
++nsuper_et_post;
k = n;
for (i = snode_start; i <= j; ++i)
k = std::min(k, inv_post(i));
l = inv_post(j);
if ( (l - k) == (j - snode_start) ) // Same number of columns in the snode
{
// This is also a supernode in the original etree
relax_end(k) = l; // Record last column
++nsuper_et;
}
else
{
for (i = snode_start; i <= j; ++i)
{
l = inv_post(i);
if (descendants(i) == 0)
{
relax_end(l) = l;
++nsuper_et;
}
}
}
j++;
// Search for a new leaf
while (descendants(j) != 0 && j < n) j++;
} // End postorder traversal of the etree
// Recover the original etree
et = et_save;
}
#endif

221
Eigen/src/SparseLU/SparseLU_panel_dfs.h Normal file
View File

@ -0,0 +1,221 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
/*
* NOTE: This file is the modified version of xpanel_dfs.c file in SuperLU
* -- SuperLU routine (version 2.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* November 15, 1997
*
* Copyright (c) 1994 by Xerox Corporation. All rights reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* Permission to modify the code and to distribute modified code is
* granted, provided the above notices are retained, and a notice that
* the code was modified is included with the above copyright notice.
*/
#ifndef SPARSELU_PANEL_DFS_H
#define SPARSELU_PANEL_DFS_H
/**
* \brief Performs a symbolic factorization on a panel of columns [jcol, jcol+w)
*
* A supernode representative is the last column of a supernode.
* The nonzeros in U[*,j] are segments that end at supernodes representatives
*
* The routine returns a list of the supernodal representatives
* in topological order of the dfs that generates them. This list is
* a superset of the topological order of each individual column within
* the panel.
* The location of the first nonzero in each supernodal segment
* (supernodal entry location) is also returned. Each column has
* a separate list for this purpose.
*
* Two markers arrays are used for dfs :
* marker[i] == jj, if i was visited during dfs of current column jj;
* marker1[i] >= jcol, if i was visited by earlier columns in this panel;
*
* \param m number of rows in the matrix
* \param w Panel size
* \param jcol Starting column of the panel
* \param A Input matrix in column-major storage
* \param perm_r Row permutation
* \param nseg
*
*/
template <typename MatrixType, typename VectorType>
int SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, VectorXi& perm_r, VectorXi& nseg, int& nseg, VectorType& dense, VectorXi& panel_lsub, VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, VectorXi& marker, VectorXi& parent, VectorXi& xplore, LU_GlobalLu_t& Glu)
{
int jj; // Index through each column in the panel
int nextl_col; // Next available position in panel_lsub[*,jj]
int krow; // Row index of the current element
int kperm; // permuted row index
int krep; // Supernode reprentative of the current row
int kmark;
int chperm, chmark, chrep, oldrep, kchild;
int myfnz; // First nonzero element in the current column
int xdfs, maxdfs, kpar;
// Initialize pointers
// VectorXi& marker1 = marker.block(m, m);
VectorBlock<VectorXi> marker1(marker, m, m);
nseg = 0;
VectorXi& xsup = Glu.xsup;
VectorXi& supno = Glu.supno;
VectorXi& lsub = Glu.lsub;
VectorXi& xlsub = Glu.xlsub;
// For each column in the panel
for (jj = jcol; jj < jcol + w; jj++)
{
nextl_col = (jj - jcol) * m;
//FIXME
VectorBlock<VectorXi> repfnz_col(repfnz.segment(nextl_col, m)); // First nonzero location in each row
VectorBlock<VectorXi> dense_col(dense.segment(nextl_col, m)); // Accumulate a column vector here
// For each nnz in A[*, jj] do depth first search
for (MatrixType::InnerIterator it(A, jj); it; ++it)
{
krow = it.row();
dense_col(krow) = it.val();
kmark = marker(krow);
if (kmark == jj)
continue; // krow visited before, go to the next nonzero
// For each unmarked krow of jj
marker(krow) = jj;
kperm = perm_r(krow);
if (kperm == -1 ) {
// krow is in L : place it in structure of L(*, jj)
panel_lsub(nextl_col++) = krow; // krow is indexed into A
}
else
{
// krow is in U : if its supernode-representative krep
// has been explored, update repfnz(*)
krep = xsup(supno(kperm)+1) - 1;
myfnz = repfnz_col(krep);
if (myfnz != -1 )
{
// Representative visited before
if (myfnz > kperm ) repfnz_col(krep) = kperm;
}
else
{
// Otherwise, perform dfs starting at krep
oldrep = -1;
parent(krep) = oldrep;
repfnz_col(krep) = kperm;
xdfs = xlsub(krep);
maxdfs = xprune(krep);
do
{
// For each unmarked kchild of krep
while (xdfs < maxdfs)
{
kchild = lsub(xdfs);
xdfs++;
chmark = marker(kchild);
if (chmark != jj )
{
marker(kchild) = jj;
chperm = perm_r(kchild);
if (chperm == -1)
{
// case kchild is in L: place it in L(*, j)
panel_lsub(nextl_col++) = kchild;
}
else
{
// case kchild is in U :
// chrep = its supernode-rep. If its rep has been explored,
// update its repfnz(*)
chrep = xsup(supno(chperm)+1) - 1;
myfnz = repfnz_col(chrep);
if (myfnz != -1)
{ // Visited before
if (myfnz > chperm)
repfnz_col(chrep) = chperm;
}
else
{ // Cont. dfs at snode-rep of kchild
xplore(krep) = xdfs;
oldrep = krep;
krep = chrep; // Go deeper down G(L)
parent(krep) = oldrep;
repfnz_col(krep) = chperm;
xdfs = xlsub(krep);
maxdfs = xprune(krep);
} // end if myfnz != -1
} // end if chperm == -1
} // end if chmark !=jj
} // end while xdfs < maxdfs
// krow has no more unexplored nbrs :
// Place snode-rep krep in postorder DFS, if this
// segment is seen for the first time. (Note that
// "repfnz(krep)" may change later.)
// Baktrack dfs to its parent
if (marker1(krep) < jcol )
{
segrep(nseg) = krep;
++nseg;
marker1(krep) = jj;
}
kpar = parent(krep); // Pop recursion, mimic recursion
if (kpar == -1)
break; // dfs done
krep = kpar;
xdfs = xplore(krep);
maxdfs = xprune(krep);
} while (kpar != -1); // Do until empty stack
} // end if (myfnz = -1)
} // end if (kperm == -1)
}// end for nonzeros in column jj
} // end for column jj
}
#endif

132
Eigen/src/SparseLU/SparseLU_pivotL.h Normal file
View File

@ -0,0 +1,132 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
/*
* NOTE: This file is the modified version of dpivotL.c file in SuperLU
* -- SuperLU routine (version 3.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* October 15, 2003
*
* Copyright (c) 1994 by Xerox Corporation. All rights reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* Permission to modify the code and to distribute modified code is
* granted, provided the above notices are retained, and a notice that
* the code was modified is included with the above copyright notice.
*/
#ifndef SPARSELU_PIVOTL_H
#define SPARSELU_PIVOTL_H
/**
* \brief Performs the numerical pivotin on the current column of L, and the CDIV operation.
*
* Here is the pivot policy :
* (1)
*
* \param jcol The current column of L
* \param pivrow [out] The pivot row
*
*
*/
int SparseLU::LU_pivotL(const int jcol, Index& pivrow)
{
// Initialize pointers
VectorXi& lsub = m_Glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??)
VectorXi& xlsub = m_Glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*)
Scalar* lusup = m_Glu.lusup.data(); // Numerical values of the rectangular supernodes
VectorXi& xlusup = m_Glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*)
Index fsupc = (m_Glu.xsup)((m_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
Index lptr = xlsub(fsupc); // pointer to the starting location of the row subscripts for this supernode portion
Index nsupr = xlsub(fsupc+1) - lptr; // Number of rows in the supernode
Scalar* lu_sup_ptr = &(lusup[xlusup(fsupc)]); // Start of the current supernode
Scalar* lu_col_ptr = &(lusup[xlusup(jcol)]); // Start of jcol in the supernode
Index* lsub_ptr = &(lsub.data()[lptr]); // Start of row indices of the supernode
// Determine the largest abs numerical value for partial pivoting
Index diagind = m_iperm_c(jcol); // diagonal index
Scalar pivmax = 0.0;
Index pivptr = nsupc;
Index diag = -1;
Index old_pivptr = nsupc;
Scalar rtemp;
for (isub = nsupc; isub < nsupr; ++isub) {
rtemp = std::abs(lu_col_ptr[isub]);
if (rtemp > pivmax) {
pivmax = rtemp;
pivptr = isub;
}
if (lsub_ptr[isub] == diagind) diag = isub;
}
// Test for singularity
if ( pivmax == 0.0 ) {
pivrow = lsub_ptr[pivptr];
m_perm_r(pivrow) = jcol;
return (jcol+1);
}
Scalar thresh = m_diagpivotthresh * pivmax;
// Choose appropriate pivotal element
{
// Test if the diagonal element can be used as a pivot (given the threshold value)
if (diag >= 0 )
{
// Diagonal element exists
rtemp = std::abs(lu_col_ptr[diag]);
if (rtemp != Scalar(0.0) && rtemp >= thresh) pivptr = diag;
}
pivrow = lsub_ptr[pivptr];
}
// Record pivot row
perm_r(pivrow) = jcol;
// Interchange row subscripts
if (pivptr != nsupc )
{
std::swap( lsub_ptr(pivptr), lsub_ptr(nsupc) );
// Interchange numerical values as well, for the two rows in the whole snode
// such that L is indexed the same way as A
for (icol = 0; icol <= nsupc; icol++)
{
itemp = pivptr + icol * nsupr;
std::swap(lu_sup_ptr[itemp], lu_sup_ptr[nsupc + icol * nsupr]);
}
}
// cdiv operations
Scalar temp = Scalar(1.0) / lu_col_ptr[nsupc];
for (k = nsupc+1; k < nsupr; k++)
lu_col_ptr[k] *= temp;
return 0;
}
#endif

89
Eigen/src/SparseLU/SparseLU_relax_snode.h Normal file
View File

@ -0,0 +1,89 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
/* This file is a modified version of heap_relax_snode.c file in SuperLU
* -- SuperLU routine (version 3.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* October 15, 2003
*
* Copyright (c) 1994 by Xerox Corporation. All rights reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* Permission to modify the code and to distribute modified code is
* granted, provided the above notices are retained, and a notice that
* the code was modified is included with the above copyright notice.
*/
#ifndef EIGEN_HEAP_RELAX_SNODE_H
#define EIGEN_HEAP_RELAX_SNODE_H
#include <Eigen_coletree.h>
/**
* \brief Identify the initial relaxed supernodes
*
* This routine applied to a column elimination tree.
* It assumes that the matrix has been reordered according to the postorder of the etree
* \param et elimination tree
* \param relax_columns Maximum number of columns allowed in a relaxed snode
* \param descendants Number of descendants of each node in the etree
* \param relax_end last column in a supernode
*/
void internal::LU_relax_snode (const int n, VectorXi& et, const int relax_columns, VectorXi& descendants, VectorXi& relax_end)
{
// compute the number of descendants of each node in the etree
register int j, parent;
relax_end.setConstant(-1);
descendants.setZero();
for (j = 0; j < n; j++)
{
parent = et(j);
if (parent != n) // not the dummy root
descendants(parent) += descendants(j) + 1;
}
// Identify the relaxed supernodes by postorder traversal of the etree
register int snode_start; // beginning of a snode
for (j = 0; j < n; )
{
parent = et(j);
snode_start = j;
while ( parent != n && descendants(parent) < relax_columns )
{
j = parent;
parent = et(j);
}
// Found a supernode in postordered etree, j is the last column
relax_end(snode_start) = j; // Record last column
j++;
// Search for a new leaf
while (descendants(j) != 0 && j < n) j++;
} // End postorder traversal of the etree
}
#endif

88
Eigen/src/SparseLU/SparseLU_snode_bmod.h Normal file
View File

@ -0,0 +1,88 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
/*
* NOTE: This file is the modified version of dsnode_bmod.c file in SuperLU
* -- SuperLU routine (version 3.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* October 15, 2003
*
* Copyright (c) 1994 by Xerox Corporation. All rights reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* Permission to modify the code and to distribute modified code is
* granted, provided the above notices are retained, and a notice that
* the code was modified is included with the above copyright notice.
*/
#ifndef SPARSELU_SNODE_BMOD_H
#define SPARSELU_SNODE_BMOD_H
template <typename VectorType>
int SparseLU::LU_dsnode_bmod (const int jcol, const int jsupno, const int fsupc,
VectorType& dense, VectorType& tempv)
{
VectorXi& lsub = m_Glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??)
VectorXi& xlsub = m_Glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*)
Scalar* lusup = m_Glu.lusup.data(); // Numerical values of the rectangular supernodes
VectorXi& xlusup = m_Glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*)
int nextlu = xlusup(jcol); // Starting location of the next column to add
int irow;
// Process the supernodal portion of L\U[*,jcol]
for (int isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++)
{
irow = lsub(isub);
lusup(nextlu) = dense(irow);
dense(irow) = 0;
++nextlu;
}
xlusup(jcol + 1) = nextlu; // Initialize xlusup for next column ( jcol+1 )
if (fsupc < jcol ){
int luptr = xlusup(fsupc); // points to the first column of the supernode
int nsupr = xlsub(fsupc + 1) -xlsub(fsupc); //Number of rows in the supernode
int nsupc = jcol - fsupc; // Number of columns in the supernodal portion of L\U[*,jcol]
int ufirst = xlusup(jcol); // points to the beginning of column jcol in supernode L\U(jsupno)
int nrow = nsupr - nsupc; // Number of rows in the off-diagonal blocks
int incx = 1, incy = 1;
Scalar alpha = Scalar(-1.0);
Scalar beta = Scalar(1.0);
// Solve the triangular system for U(fsupc:jcol, jcol) with L(fspuc..., fsupc:jcol)
//BLASFUNC(trsv)("L", "N", "U", &nsupc, &(lusup[luptr]), &nsupr, &(lusup[ufirst]), &incx);
Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup[luptr]), nsupc, nsupc, OuterStride<>(nsupr) );
Map<Matrix<Scalar,Dynamic,1> > l(&(lusup[ufirst]), nsupc);
l = A.triangularView<Lower>().solve(l);
// Update the trailing part of the column jcol U(jcol:jcol+nrow, jcol) using L(jcol:jcol+nrow, fsupc:jcol) and U(fsupc:jcol)
BLASFUNC(gemv)("N", &nrow, &nsupc, &alpha, &lusup[luptr+nsupc], &nsupr, &lusup[ufirst], &incx, &beta, &lusup[ufirst+nsupc], &incy);
return 0;
}
#endif

119
Eigen/src/SparseLU/SparseLU_snode_dfs.h Normal file
View File

@ -0,0 +1,119 @@
// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
/*
* NOTE: This file is the modified version of dsnode_dfs.c file in SuperLU
* -- SuperLU routine (version 2.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* November 15, 1997
*
* Copyright (c) 1994 by Xerox Corporation. All rights reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* Permission to modify the code and to distribute modified code is
* granted, provided the above notices are retained, and a notice that
* the code was modified is included with the above copyright notice.
*/
#ifdef EIGEN_SNODE_DFS_H
#define EIGEN_SNODE_DFS_H
/**
* \brief Determine the union of the row structures of those columns within the relaxed snode.
* NOTE: The relaxed snodes are leaves of the supernodal etree, therefore,
* the portion outside the rectangular supernode must be zero.
*
* \param jcol start of the supernode
* \param kcol end of the supernode
* \param asub Row indices
* \param colptr Pointer to the beginning of each column
* \param xprune (out) The pruned tree ??
* \param marker (in/out) working vector
*/
template <typename Index>
int SparseLU::LU_snode_dfs(const int jcol, const int kcol, const VectorXi* asub, const VectorXi* colptr,
VectorXi& xprune, VectorXi& marker, LU_GlobalLu_t *m_Glu)
{
VectorXi& xsup = m_Glu->xsup;
VectorXi& supno = m_Glu->supno; // Supernode number corresponding to this column
VectorXi& lsub = m_Glu->lsub;
VectorXi& xlsub = m_Glu->xlsub;
int nsuper = ++supno(jcol); // Next available supernode number
register int nextl = xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
register int i,k;
int krow,kmark;
for (i = jcol; i <=kcol; i++)
{
// For each nonzero in A(*,i)
for (k = colptr(i); k < colptr(i+1); k++)
{
krow = asub(k);
kmark = marker(krow);
if ( kmark != kcol )
{
// First time to visit krow
marker(krow) = kcol;
lsub(nextl++) = krow;
if( nextl >= nzlmax )
{
m_Glu->lsub = LUMemXpand<Index>(jcol, nextl, LSUB, nzlmax);
m_Glu->nzlmax = nzlmax;
lsub = m_Glu->lsub;
}
}
}
supno(i) = nsuper;
}
// If supernode > 1, then make a copy of the subscripts for pruning
if (jcol < kcol)
{
int new_next = nextl + (nextl - xlsub(jcol));
while (new_next > nzlmax)
{
m_Glu->lsub = LUMemXpand<Index>(jcol, nextl, LSUB, &nzlmax);
m_Glu->nzlmax= nzlmax;
lsub = m_Glu->lsub;
}
int ifrom, ito = nextl;
for (ifrom = xlsub(jcol); ifrom < nextl;)
lsub(ito++) = lsub(ifrom++);
for (i = jcol+1; i <=kcol; i++)xlsub(i) = nextl;
nextl = ito;
}
xsup(nsuper+1) = kcol + 1; // Start of next available supernode
supno(kcol+1) = nsuper;
xprune(kcol) = nextl;
xlsub(kcol+1) = nextl;
return 0;
}
#endif