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

Eliminate and prune columns in a panel

Browse Source
This commit is contained in:
Desire NUENTSA 2012-05-31 17:10:29 +02:00
parent 8608d08d65
commit b26d6b02de
5 changed files with 412 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
Eigen/src/SparseLU/SparseLU.h
View File

@ -113,7 +113,7 @@ class SparseLU
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
int nnzL, nnzU; // Nonzeros in L and U factors
private:
// Copy constructor
SparseLU (SparseLU& ) {}
@ -260,6 +260,7 @@ void SparseLU::factorize(const MatrixType& matrix)
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
int irep,ir;
for (jcol = 0; jcol < min_mn; )
{
if (relax_end(jcol) != IND_EMPTY)
@ -297,7 +298,7 @@ void SparseLU::factorize(const MatrixType& matrix)
LU_snode_bmod(icol, jsupno, fsupc, dense, tempv);
// Eliminate the current column
info = LU_pivotL(icol, pivrow);
info = LU_pivotL(icol, m_diagpivotthresh, m_perm_r, m_iperm_c, pivrow, m_Glu);
if ( !info )
{
m_info = NumericalIssue;
@ -357,10 +358,17 @@ void SparseLU::factorize(const MatrixType& matrix)
}
// Copy the U-segments to ucol(*)
//FIXME Check that repfnz_k, dense_k... have stored references to modified columns
info = LU_copy_to_col(jj, nseg, segrep, repfnz_k, perm_r, dense_k, m_Glu);
if ( !info )
{
m_info = NumericalIssue;
m_factorizationIsOk = false;
return;
}
// Form the L-segment
info = LU_pivotL(...);
info = LU_pivotL(jj, m_diagpivotthresh, m_perm_r, iperm_c, pivrow, m_Glu);
if ( !info )
{
m_info = NumericalIssue;
@ -369,11 +377,44 @@ void SparseLU::factorize(const MatrixType& matrix)
}
// Prune columns (0:jj-1) using column jj
LU_pruneL(jj, m_perm_r, pivrow, nseg, segrep, repfnz_k, xprune, m_Glu);
} // end for
// Reset repfnz for this column
for (i = 0; i < nseg; i++)
{
irep = segrep(i);
repfnz(irep) = IND_EMPTY;
}
} // end SparseLU within the panel
jcol += panel_size; // Move to the next panel
} // end else
} // end for -- end elimination
// Adjust row permutation in the case of rectangular matrices
if (m > n )
{
k = 0;
for (i = 0; i < m; ++i)
{
if ( perm_r(i) == IND_EMPTY )
{
perm_r(i) = n + k;
++k;
}
}
}
// Count the number of nonzeros in factors
LU_countnz(min_mn, xprune, m_nnzL, m_nnzU, m_Glu);
// Apply permutation to the L subscripts
LU_fixupL(min_mn, m_perm_r, m_Glu);
// Free work space and compress storage iwork and work
// ?? Should it be done automatically by C++
//...
// Create supernode matrix L and the column major matrix U
// ...
m_info = Success;
m_factorizationIsOk = ok;
}

65
Eigen/src/SparseLU/SparseLU_Utils.h
View File

@ -29,4 +29,67 @@
#define LU_NO_MARKER 3
#define LU_NUM_TEMPV(m,w,t,b) (std::max(m, (t+b)*w) )
#define IND_EMPTY (-1)
#endif
void SparseLU::LU_countnz(const int n, VectorXi& xprune, int& nnzL, int& nnzU, GlobalLU_t& Glu)
{
VectorXi& xsup = Glu.xsup;
VectorXi& xlsub = Glu.xlsub;
nnzL = 0;
nnzU = (Glu.xusub)(n);
int nnzL0 = 0;
int nsuper = (Glu.supno)(n);
int jlen, irep;
if (n <= 0 ) return;
// For each supernode
for (i = 0; i <= nsuper; i++)
{
fsupc = xsup(i);
jlen = xlsub(fsupc+1) - xlsub(fsupc);
for (j = fsupc; j < xsup(i+1); j++)
{
nnzL += jlen;
nnzLU += j - fsupc + 1;
jlen--;
}
irep = xsup(i+1) - 1;
nnzL0 += xprune(irep) - xlsub(irep);
}
}
/**
* \brief Fix up the data storage lsub for L-subscripts.
*
* It removes the subscripts sets for structural pruning,
* and applies permutation to the remaining subscripts
*
*/
void SparseLU::LU_fixupL(const int n, const VectorXi& perm_r, GlobalLU_t& Glu)
{
int nsuper, fsupc, i, j, k, jstart;
VectorXi& xsup = GLu.xsup;
VectorXi& lsub = Glu.lsub;
VectorXi& xlsub = Glu.xlsub;
int nextl = 0;
int nsuper = (Glu.supno)(n);
// For each supernode
for (i = 0; i <= nsuper; i++)
{
fsupc = xsup(i);
jstart = xlsub(fsupc);
xlsub(fsupc) = nextl;
for (j = jstart; j < xlsub(fsupc + 1); j++)
{
lsub(nextl) = perm_r(lsub(j)); // Now indexed into P*A
nextl++
}
for (k = fsupc+1; k < xsup(i+1); k++)
xlsub(k) = nextl; // other columns in supernode i
}
xlsub(n) = nextl;
}
#endif

123
Eigen/src/SparseLU/SparseLU_copy_to_ucol.h Normal file
View File

@ -0,0 +1,123 @@
// 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 xcopy_to_ucol.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_COPY_TO_UCOL_H
#define SPARSELU_COPY_TO_UCOL_H
/**
* \brief Performs numeric block updates (sup-col) in topological order
*
* \param jcol current column to update
* \param nseg Number of segments in the U part
* \param segrep segment representative ...
* \param repfnz ??? First nonzero column in each row ??? ...
* \param perm_r Row permutation
* \param dense Store the full representation of the column
* \param Glu Global LU data.
* \return 0 - successful return
* > 0 - number of bytes allocated when run out of space
*
*/
template <typename VectorType>
int SparseLU::LU_copy_to_ucol(const int jcol, const int nseg, VectorXi& segrep, VectorXi& repfnz, VectorXi& perm_r, VectorType& dense, LU_GlobalLu_t& Glu)
{
int ksupno, k, ksub, krep, ksupno;
VectorXi& xsup = Glu.xsup;
VectorXi& supno = Glu.supno;
VectorXi& lsub = Glu.lsub;
VectorXi& xlsub = Glu.xlsub;
VectorType& ucol = GLu.ucol;
VectorXi& usub = Glu.usub;
VectorXi& xusub = Glu.xusub;
int nzumax = GLu.nzumax;
int jsupno = supno(jcol);
// For each nonzero supernode segment of U[*,j] in topological order
k = nseg - 1;
int nextu = xusub(jcol);
int kfnz, isub, segsize;
int new_next,irow;
for (ksub = 0; ksub < nseg; ksub++)
{
krep = segrep(k); k--;
ksupno = supno(krep);
if (jsupno != ksupno ) // should go into ucol();
{
kfnz = repfnz(krep);
if (kfnz != IND_EMPTY)
{ // Nonzero U-segment
fsupc = xsup(ksupno);
isub = xlsub(fsupc) + kfnz - fsupc;
segsize = krep - kfnz + 1;
new_next = nextu + segsize;
while (new_next > nzumax)
{
Glu.ucol = LU_MemXpand<Scalar>(jcol, nextu, UCOL, nzumax); //FIXME try and catch errors
ucol = Glu.ucol;
Glu.nzumax = nzumax;
Glu.usub = LU_MemXpand<Index>(jcol, nextu, USUB, nzumax); //FIXME try and catch errors
Glu.nzumax = nzumax;
usub = Glu.usub;
lsub = Glu.lsub;
}
for (i = 0; i < segsize; i++)
{
irow = lsub(isub);
usub(nextu) = perm_r(irow); // Unlike teh L part, the U part is stored in its final order
ucol(nextu) = dense(irow);
dense(irow) = Scalar(0.0);
nextu++;
isub++;
}
} // end nonzero U-segment
} // end if jsupno
} // end for each segment
xusub(jcol + 1) = nextu; // close U(*,jcol)
return 0;
}
#endif

41
Eigen/src/SparseLU/SparseLU_pivotL.h
View File

@ -24,7 +24,7 @@
/*
* NOTE: This file is the modified version of dpivotL.c file in SuperLU
* NOTE: This file is the modified version of xpivotL.c file in SuperLU
* -- SuperLU routine (version 3.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
@ -47,23 +47,36 @@
/**
* \brief Performs the numerical pivotin on the current column of L, and the CDIV operation.
*
* Here is the pivot policy :
* (1)
* Pivot policy :
* (1) Compute thresh = u * max_(i>=j) abs(A_ij);
* (2) IF user specifies pivot row k and abs(A_kj) >= thresh THEN
* pivot row = k;
* ELSE IF abs(A_jj) >= thresh THEN
* pivot row = j;
* ELSE
* pivot row = m;
*
* Note: If you absolutely want to use a given pivot order, then set u=0.0.
*
* \param jcol The current column of L
* \param pivrow [out] The pivot row
*
* \param u diagonal pivoting threshold
* \param [in,out]perm_r Row permutation (threshold pivoting)
* \param [in] iperm_c column permutation - used to finf diagonal of Pc*A*Pc'
* \param [out]pivrow The pivot row
* \param Glu Global LU data
* \return 0 if success, i > 0 if U(i,i) is exactly zero
*
*/
int SparseLU::LU_pivotL(const int jcol, Index& pivrow)
template <typename Scalar>
int SparseLU::LU_pivotL(const int jcol, const Scalar u, VectorXi& perm_r, VectorXi& iperm_c, int& pivrow, GlobalLU_t& Glu)
{
// 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(*)
VectorXi& lsub = Glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??)
VectorXi& xlsub = Glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*)
Scalar* lusup = Glu.lusup.data(); // Numerical values of the rectangular supernodes
VectorXi& xlusup = 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 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
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
@ -72,7 +85,7 @@ int SparseLU::LU_pivotL(const int jcol, Index& pivrow)
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
Index diagind = iperm_c(jcol); // diagonal index
Scalar pivmax = 0.0;
Index pivptr = nsupc;
Index diag = -1;
@ -90,11 +103,11 @@ int SparseLU::LU_pivotL(const int jcol, Index& pivrow)
// Test for singularity
if ( pivmax == 0.0 ) {
pivrow = lsub_ptr[pivptr];
m_perm_r(pivrow) = jcol;
perm_r(pivrow) = jcol;
return (jcol+1);
}
Scalar thresh = m_diagpivotthresh * pivmax;
Scalar thresh = diagpivotthresh * pivmax;
// Choose appropriate pivotal element

152
Eigen/src/SparseLU/SparseLU_pruneL.h Normal file
View File

@ -0,0 +1,152 @@
// 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 xpruneL.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_PRUNEL_H
#define SPARSELU_PRUNEL_H
/**
* \brief Prunes the L-structure.
*
* It prunes the L-structure of supernodes whose L-structure constains the current pivot row "pivrow"
*
*
* \param jcol The current column of L
* \param [in]perm_r Row permutation
* \param [out]pivrow The pivot row
* \param nseg Number of segments ???
* \param segrep
* \param repfnz
* \param [out]xprune
* \param Glu Global LU data
*
*/
template <typename VectorType>
void SparseLU::LU_pruneL(const int jcol, const VectorXi& perm_r, const int pivrow, const int nseg, const VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, GlobalLU_t& Glu)
{
// Initialize pointers
VectorXi& xsup = Glu.xsup;
VectorXi& supno = Glu.supno;
VectorXi& lsub = Glu.lsub;
VectorXi& xlsub = Glu.xlsub;
VectorType& lusup = Glu.lusup;
VectorXi& xlusup = Glu.xlusup;
// For each supernode-rep irep in U(*,j]
int jsupno = supno(jcol);
int i,irep,irep1;
bool movnum, do_prune = false;
int kmin, kmax, ktemp, minloc, maxloc;
for (i = 0; i < nseg; i++)
{
irep = segrep(i);
irep1 = irep + 1;
do_prune = false;
// Don't prune with a zero U-segment
if (repfnz(irep) == IND_EMPTY) continue;
// If a snode overlaps with the next panel, then the U-segment
// is fragmented into two parts -- irep and irep1. We should let
// pruning occur at the rep-column in irep1s snode.
if (supno(irep) == supno(irep1) continue; // don't prune
// If it has not been pruned & it has a nonz in row L(pivrow,i)
if (supno(irep) != jsupno )
{
if ( xprune (irep) >= xlsub(irep1)
{
kmin = xlsub(irep);
kmax = xlsub(irep1) - 1;
for (krow = kmin; krow <= kmax; krow++)
{
if (lsub(krow) == pivrow)
{
do_prune = true;
break;
}
}
}
if (do_prune)
{
// do a quicksort-type partition
// movnum=true means that the num values have to be exchanged
movnum = false;
if (irep == xsup(supno(irep)) ) // Snode of size 1
movnum = true;
while (kmin <= kmax)
{
if (perm_r(lsub(kmax)) == IND_EMPTY)
kmax--;
else if ( perm_r(lsub(kmin)) != IND_EMPTY)
kmin--;
else
{
// kmin below pivrow (not yet pivoted), and kmax
// above pivrow: interchange the two suscripts
ktemp = lsub(kmin);
lsub(kmin) = lsub(kmax);
lsub(kmax) = ktemp;
// If the supernode has only one column, then we
// only keep one set of subscripts. For any subscript
// intercnahge performed, similar interchange must be
// done on the numerical values.
if (movnum)
{
minloc = xlusup(irep) + ( kmin - xlsub(irep) );
maxloc = xlusup(irep) + ( kmax - xlsub(irep) );
std::swap(lusup(minloc), lusup(maxloc));
}
kmin++;
kmax--;
}
} // end while
xprune(irep) = kmin;
} // end if do_prune
} // end pruning
} // End for each U-segment
}
#endif