diff options
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/src/SparseLU/SparseLU.h | 31 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_panel_bmod.h | 178 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_panel_dfs.h | 5 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_snode_bmod.h | 26 |
4 files changed, 221 insertions, 19 deletions
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index f5a1c787e..f1c530b55 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -192,7 +192,7 @@ void SparseLU::analyzePattern(const MatrixType& mat) * the estimated amount of space needed, plus A->ncol. */ template <typename MatrixType> -int SparseLU::factorize(const MatrixType& matrix) +void SparseLU::factorize(const MatrixType& matrix) { // Allocate storage common to the factor routines @@ -256,7 +256,6 @@ int SparseLU::factorize(const MatrixType& matrix) 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 @@ -272,8 +271,9 @@ int SparseLU::factorize(const MatrixType& matrix) info = LU_snode_dfs(jcol, kcol, m_mat.innerIndexPtr(), m_mat.outerIndexPtr(), xprune, marker); if ( !info ) { - ok = false; - break; + m_info = NumericalIssue; + m_factorizationIsOk = false; + return; } nextu = xusub(jcol); //starting location of column jcol in ucol nextlu = xlusup(jcol); //Starting location of column jcol in lusup (rectangular supernodes) @@ -322,17 +322,36 @@ int SparseLU::factorize(const MatrixType& matrix) 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); + LU_panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_Glu); // 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 + + nseg = nseg1; // begin after all the panel segments + //Depth-first-search for the current column + info = LU_column_dfs(m, jj, ... ); + if ( !info ) + { + m_info = NumericalIssue; + m_factorizationIsOk = false; + return; + } + // Numeric updates to this column + info = LU_column_bmod(jj, ... ); + if ( !info ) + { + m_info = NumericalIssue; + m_factorizationIsOk = false; + return; + } + } // end for jcol += panel_size; // Move to the next panel } // end else } // end for -- end elimination - m_info = ok ? Success : NumericalIssue; + m_info = Success; m_factorizationIsOk = ok; } diff --git a/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/Eigen/src/SparseLU/SparseLU_panel_bmod.h new file mode 100644 index 000000000..29cc6d0f0 --- /dev/null +++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h @@ -0,0 +1,178 @@ +// 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 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_PANEL_BMOD_H +#define SPARSELU_PANEL_BMOD_H +/** + * \brief Performs numeric block updates (sup-panel) in topological order. + * + * Before entering this routine, the original nonzeros in the panel + * were already copied i nto the spa[m,w] ... FIXME to be checked + * + * \param m number of rows in the matrix + * \param w Panel size + * \param jcol Starting column of the panel + * \param nseg Number of segments in the U part + * \param dense Store the full representation of the panel + * \param tempv working array + * \param segrep in ... + * \param repfnz in ... + * \param Glu Global LU data. + * + * + */ +template <typename VectorType> +void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, VectorType& dense, VectorType& tempv, VectorXi& segrep, VectorXi& repfnz, LU_GlobalLu_t& Glu) +{ + VectorXi& xsup = Glu.xsup; + VectorXi& supno = Glu.supno; + VectorXi& lsub = Glu.lsub; + VectorXi& xlsub = Glu.xlsub; + VectorXi& xlusup = Glu.xlusup; + VectorType& lusup = Glu.lusup; + + int i,ksub,jj,nextl_col,irow; + int fsupc, nsupc, nsupr, nrow; + int krep, krep_ind; + int nrow; + int lptr; // points to the row subscripts of a supernode + int luptr; // ... + int segsze,no_zeros,irow ; + // For each nonz supernode segment of U[*,j] in topological order + int k = nseg - 1; + for (ksub = 0; ksub < nseg; ksub++) + { // For each updating supernode + + /* krep = representative of current k-th supernode + * fsupc = first supernodal column + * nsupc = number of columns in a supernode + * nsupr = number of rows in a supernode + */ + krep = segrep(k); k--; + fsupc = xsup(supno(krep)); + nsupc = krep - fsupc + 1; + nsupr = xlsub(fsupc+1) - xlsub(fsupc); + nrow = nsupr - nsupc; + lptr = xlsub(fsupc); + krep_ind = lptr + nsupc - 1; + + repfnz_col = repfnz; + dense_col = dense; + + // NOTE : Unlike the original implementation in SuperLU, the present implementation + // does not include a 2-D block update. + + // Sequence through each column in the panel + for (jj = jcol; jj < jcol + w; jj++) + { + nextl_col = (jj-jcol) * m; + VectorBlock<VectorXi> repfnz_col(repfnz.segment(nextl_col, m)); // First nonzero column index for each row + VectorBLock<VectorXi> dense_col(dense.segment(nextl_col, m)); // Scatter/gather entire matrix column from/to here + + kfnz = repfnz_col(krep); + if ( kfnz == -1 ) + continue; // skip any zero segment + + segsize = krep - kfnz + 1; + luptr = xlusup(fsupc); + + // NOTE : Unlike the original implementation in SuperLU, + // there is no update feature for col-col, 2col-col ... + + // Perform a trianglar solve and block update, + // then scatter the result of sup-col update to dense[] + no_zeros = kfnz - fsupc; + + // Copy U[*,j] segment from dense[*] to tempv[*] : + // The result of triangular solve is in tempv[*]; + // The result of matric-vector update is in dense_col[*] + isub = lptr + no_zeros; + for (i = 0; i < segsize; ++i) + { + irow = lsub(isub); + tempv(i) = dense_col(irow); // Gather to a compact vector + ++isub; + } + // Start effective triangle + luptr += nsupr * no_zeros + no_zeros; + // triangular solve with Eigen + Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); + Map<Matrix<Scalar,Dynamic,1> > u( tempv.data(), segsize); + u = A.triangularView<Lower>().solve(u); + + luptr += segsize; + // Dense Matrix vector product y <-- A*x; + new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); + Map<VectorType> l( &(tempv.data()[segsize]), segsize); + l= A * u; + + // Scatter tempv(*) into SPA dense(*) such that tempv(*) + // can be used for the triangular solve of the next + // column of the panel. The y will be copied into ucol(*) + // after the whole panel has been finished. + + isub = lptr + no_zeros; + for (i = 0; i < segsize; i++) + { + irow = lsub(isub); + dense_col(irow) = tempv(i); + tempv(i) = zero; + isub++; + } + + // Scatter the update from &tempv[segsize] into SPA dense(*) + // Start dense rectangular L + for (i = 0; i < nrow; i++) + { + irow = lsub(isub); + dense_col(irow) -= tempv(segsize + i); + tempv(segsize + i) = 0; + ++isub; + } + + } // End for each column in the panel + + } // End for each updating supernode +} +#endif
\ No newline at end of file diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h index 550544d05..7b85b6d7c 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h @@ -67,11 +67,12 @@ * \param jcol Starting column of the panel * \param A Input matrix in column-major storage * \param perm_r Row permutation - * \param nseg + * \param nseg Number of U segments + * ... * */ 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) +void 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 diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h index fc6ffc320..e7146a262 100644 --- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h @@ -46,12 +46,12 @@ #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) + VectorType& dense, VectorType& tempv, LU_GlobalLu_t& Glu) { - 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(*) + VectorType& lusup = Glu.lusup; // Numerical values of the rectangular supernodes + VectorXi& xlusup = 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; @@ -72,16 +72,20 @@ int SparseLU::LU_dsnode_bmod (const int jcol, const int jsupno, const int fsupc, 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); +// 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); + Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); + Map<Matrix<Scalar,Dynamic,1> > l(&(lusup.data()[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); + Map<Matrix<Scalar,Dynamic,1> > u(&(lusup.data()[ufirst+nsupc], nsupc); + u = A * l; +// BLASFUNC(gemv)("N", &nrow, &nsupc, &alpha, &lusup[luptr+nsupc], &nsupr, &lusup[ufirst], &incx, &beta, &lusup[ufirst+nsupc], &incy); return 0; } |