diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 5b45dd6d0..833832f3b 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -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; } diff --git a/Eigen/src/SparseLU/SparseLU_Utils.h b/Eigen/src/SparseLU/SparseLU_Utils.h index 27eaed25c..88d1c8b80 100644 --- a/Eigen/src/SparseLU/SparseLU_Utils.h +++ b/Eigen/src/SparseLU/SparseLU_Utils.h @@ -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 \ No newline at end of file + +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 diff --git a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h new file mode 100644 index 000000000..3f8d8abe2 --- /dev/null +++ b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h @@ -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 +// +// 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 . + +/* + + * 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 +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(jcol, nextu, UCOL, nzumax); //FIXME try and catch errors + ucol = Glu.ucol; + Glu.nzumax = nzumax; + Glu.usub = LU_MemXpand(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 \ No newline at end of file diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h index f939ef939..3bfe14e7e 100644 --- a/Eigen/src/SparseLU/SparseLU_pivotL.h +++ b/Eigen/src/SparseLU/SparseLU_pivotL.h @@ -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 +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 diff --git a/Eigen/src/SparseLU/SparseLU_pruneL.h b/Eigen/src/SparseLU/SparseLU_pruneL.h new file mode 100644 index 000000000..687717d52 --- /dev/null +++ b/Eigen/src/SparseLU/SparseLU_pruneL.h @@ -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 +// +// 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 . + +/* + + * 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 +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 \ No newline at end of file