diff options
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU.h | 35 | ||||
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU_Memory.h | 6 | ||||
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU_Utils.h | 2 | ||||
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU_column_bmod.h | 216 | ||||
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU_column_dfs.h | 269 | ||||
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU_panel_bmod.h | 8 | ||||
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU_panel_dfs.h | 16 | ||||
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU_snode_bmod.h | 20 |
8 files changed, 538 insertions, 34 deletions
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index f1c530b55..5b45dd6d0 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -245,7 +245,7 @@ void SparseLU::factorize(const MatrixType& matrix) VectorXi& xlusup = m_GLu.xlusup; VectorXi& xusub = m_Glu.xusub; - supno(0) = -1; + supno(0) = IND_EMPTY; xsup(0) = xlsub(0) = xusub(0) = xlusup(0); int panel_size = m_panel_size; int wdef = panel_size; // upper bound on panel width @@ -262,7 +262,7 @@ void SparseLU::factorize(const MatrixType& matrix) int nseg; // Number of segments in each U-column for (jcol = 0; jcol < min_mn; ) { - if (relax_end(jcol) != -1) + if (relax_end(jcol) != IND_EMPTY) { // Starting a relaxed node from jcol kcol = relax_end(jcol); // End index of the relaxed snode @@ -298,7 +298,12 @@ void SparseLU::factorize(const MatrixType& matrix) // Eliminate the current column info = LU_pivotL(icol, pivrow); - eigen_assert(info == 0 && "The matrix is structurally singular"); + if ( !info ) + { + m_info = NumericalIssue; + m_factorizationIsOk = false; + return; + } } jcol = icol; // The last column te be eliminated } @@ -309,7 +314,7 @@ void SparseLU::factorize(const MatrixType& matrix) panel_size = w_def; for (k = jcol + 1; k < std::min(jcol+panel_size, min_mn); k++) { - if (relax_end(k) != -1) + if (relax_end(k) != IND_EMPTY) { panel_size = k - jcol; break; @@ -331,7 +336,9 @@ void SparseLU::factorize(const MatrixType& matrix) nseg = nseg1; // begin after all the panel segments //Depth-first-search for the current column - info = LU_column_dfs(m, jj, ... ); + VectorBlock<VectorXi> panel_lsubk(panel_lsub, k, m); //FIXME + VectorBlock<VectorXi> repfnz_k(repfnz, k, m); //FIXME + info = LU_column_dfs(m, jj, perm_r, nseg, panel_lsub(k), segrep, repfnz_k, xprune, marker, parent, xplore, m_Glu); if ( !info ) { m_info = NumericalIssue; @@ -339,7 +346,21 @@ void SparseLU::factorize(const MatrixType& matrix) return; } // Numeric updates to this column - info = LU_column_bmod(jj, ... ); + VectorBlock<VectorXi> dense_k(dense, k, m); //FIXME + VectorBlock<VectorXi> segrep_k(segrep, nseg1, m) // FIXME Check the length + info = LU_column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_Glu); + if ( !info ) + { + m_info = NumericalIssue; + m_factorizationIsOk = false; + return; + } + + // Copy the U-segments to ucol(*) + + + // Form the L-segment + info = LU_pivotL(...); if ( !info ) { m_info = NumericalIssue; @@ -347,6 +368,8 @@ void SparseLU::factorize(const MatrixType& matrix) return; } + // Prune columns (0:jj-1) using column jj + } // end for jcol += panel_size; // Move to the next panel } // end else diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h index 6e0fc658d..91b24fa67 100644 --- a/Eigen/src/SparseLU/SparseLU_Memory.h +++ b/Eigen/src/SparseLU/SparseLU_Memory.h @@ -86,7 +86,7 @@ int SparseLU::LUMemInit(int lwork) nzlmax = std::max(1, m_fill_ratio/4.) * annz; //??? // Return the estimated size to the user if necessary - if (lwork = -1) + if (lwork == IND_EMPTY) { estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size) + (nzlmax + nzumax) * iword + (nzlumax+nzumax) * dword + n); @@ -130,7 +130,7 @@ int SparseLU::LUMemInit(int lwork) } else // m_fact == SamePattern_SameRowPerm; { - if (lwork = -1) + if (lwork == IND_EMPTY) { estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size) + (Glu.nzlmax + Glu.nzumax) * iword + (Glu.nzlumax+Glu.nzumax) * dword + n); @@ -232,7 +232,7 @@ DestType* SparseLU::LUMemXpand(int jcol, int next, MemType mem_type, int& maxlen 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"); + eigen_assert(new_mem && "Can't expand memory"); // FIXME Should be an exception return new_mem; diff --git a/Eigen/src/SparseLU/SparseLU_Utils.h b/Eigen/src/SparseLU/SparseLU_Utils.h index 3c3b24a15..27eaed25c 100644 --- a/Eigen/src/SparseLU/SparseLU_Utils.h +++ b/Eigen/src/SparseLU/SparseLU_Utils.h @@ -28,5 +28,5 @@ // 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) +#define IND_EMPTY (-1) #endif
\ No newline at end of file diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h new file mode 100644 index 000000000..58755363d --- /dev/null +++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h @@ -0,0 +1,216 @@ +// 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 xcolumn_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_COLUMN_BMOD_H +#define SPARSELU_COLUMN_BMOD_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 dense Store the full representation of the column + * \param tempv working array + * \param segrep segment representative ... + * \param repfnz ??? First nonzero column in each row ??? ... + * \param fpanelc First column in the current panel + * \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_column_bmod(const int jcol, const int nseg, VectorType& dense, VectorType& tempv, VectorXi& segrep, VectorXi& repfnz, int fpanelc, LU_GlobalLu_t& Glu) +{ + + int jsupno, k, ksub, krep, krep_ind, ksupno; + /* 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 + * luptr = location of supernodal LU-block in storage + * kfnz = first nonz in the k-th supernodal segment + * no-zeros = no lf leading zeros in a supernodal U-segment + */ + 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 nzlumax = GLu.nzlumax; + int jsupno = supno(jcol); + // For each nonzero supernode segment of U[*,j] in topological order + k = nseg - 1; + for (ksub = 0; ksub < nseg; ksub++) + { + krep = segrep(k); k--; + ksupno = supno(krep); + if (jsupno != ksupno ) + { + // outside the rectangular supernode + fsupc = xsup(ksupno); + fst_col = std::max(fsupc, fpanelc); + + // Distance from the current supernode to the current panel; + // d_fsupc = 0 if fsupc > fpanelc + d_fsupc = fst_col - fsupc; + + luptr = xlusup(fst_col) + d_fsupc; + lptr = xlsub(fsupc) + d_fsupc; + + kfnz = repfnz(krep); + kfnz = std::max(kfnz, fpanelc); + + segsize = krep - kfnz + 1; + nsupc = krep - fst_col + 1; + nsupr = xlsub(fsupc+1) - xlsub(fsupc); + nrow = nsupr - d_fsupc - nsupc; + krep_ind = lptr + nsupc - 1; + + // NOTE Unlike the original implementation in SuperLU, the only feature + // here is a sup-col update. + + // Perform a triangular solver and block update, + // then scatter the result of sup-col update to dense + no_zeros = kfnz - fst_col; + // First, copy U[*,j] segment from dense(*) to tempv(*) + isub = lptr + no_zeros; + for (i = 0; i ww segsize; i++) + { + irow = lsub(isub); + tempv(i) = densee(irow); + ++isub; + } + // Dense triangular solve -- start effective triangle + luptr += nsupr * no_zeros + no_zeros; + // Form Eigen matrix and vector + Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); + Map<VectorType> u(tempv.data(), segsize); + u = A.triangularView<Lower>().solve(u); + + // Dense matrix-vector product y <-- A*x + luptr += segsize; + new (&A) (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); + Map<VectorType> l( &(tempv.data()[segsize]), segsize); + l= A * u; + + // Scatter tempv[] into SPA dense[] as a temporary storage + isub = lptr + no_zeros; + for (i = 0; i w segsize; i++) + { + irow = lsub(isub); + dense(irow) = tempv(i); + tempv(i) = Scalar(0.0); + ++isub; + } + + // Scatter l into SPA dense[] + for (i = 0; i < nrow; i++) + { + irow = lsub(isub); + dense(irow) -= tempv(segsize + i); + tempv(segsize + i) = Scalar(0.0); + ++isub; + } + } // end if jsupno + } // end for each segment + + // Process the supernodal portion of L\U[*,j] + nextlu = xlusup(jcol); + fsupc = xsup(jsupno); + + // copy the SPA dense into L\U[*,j] + new_next = nextlu + xlsub(fsupc + 1) - xlsub(fsupc); + while (new_next > nzlumax ) + { + Glu.lusup = LUmemXpand<Scalar>(jcol, nextlu, LUSUP, &nzlumax); + Glu.nzlumax = nzlumax; + lusup = Glu.lusup; + lsub = Glu.lsub; + } + + for (isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++) + { + irow = lsub(isub); + lusub(nextlu) = dense(irow); + dense(irow) = Scalar(0.0); + ++nextlu; + } + + xlusup(jcol + 1) = nextlu; // close L\U(*,jcol); + + /* For more updates within the panel (also within the current supernode), + * should start from the first column of the panel, or the first column + * of the supernode, whichever is bigger. There are two cases: + * 1) fsupc < fpanelc, then fst_col <- fpanelc + * 2) fsupc >= fpanelc, then fst_col <-fsupc + */ + fst_col = std::max(fsupc, fpanelc); + + if (fst_col < jcol) + { + // Distance between the current supernode and the current panel + // d_fsupc = 0 if fsupc >= fpanelc + d_fsupc = fst_col - fsupc; + + lptr = xlsub(fsupc) + d_fsupc; + luptr = xlusup(fst_col) + d_fsupc; + nsupr = xlsub(fsupc+1) - xlsub(fsupc); // leading dimension + nsupc = jcol - fst_col; // excluding jcol + nrow = nsupr - d_fsupc - nsupc; + + // points to the beginning of jcol in snode L\U(jsupno) + ufirst = xlusup(jcol) + d_fsupc; + Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); + Map<VectorType> l( &(lusup.data()[ufirst]), nsupc ); + u = A.triangularView().solve(u); + + new (&A) Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); + Map<VectorType> l( &(lusup.data()[ufirst+nsupc]), nsupr ); + l = l - A * u; + + } // End if fst_col + return 0; +} +#endif
\ No newline at end of file diff --git a/Eigen/src/SparseLU/SparseLU_column_dfs.h b/Eigen/src/SparseLU/SparseLU_column_dfs.h new file mode 100644 index 000000000..15ddcf7c0 --- /dev/null +++ b/Eigen/src/SparseLU/SparseLU_column_dfs.h @@ -0,0 +1,269 @@ +// 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 xcolumn_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_COLUMN_DFS_H +#define SPARSELU_COLUMN_DFS_H +/** + * \brief Performs a symbolic factorization on column jcol and decide the supernode boundary + * + * 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. + * The location of the first nonzero in each supernodal segment + * (supernodal entry location) is also returned. + * + * \param m number of rows in the matrix + * \param jcol Current column + * \param perm_r Row permutation + * \param [in,out] nseg Number of segments in current U[*,j] - new segments appended + * \param lsub_col defines the rhs vector to start the dfs + * \param [in,out] segrep Segment representatives - new segments appended + * \param repfnz + * \param xprune + * \param marker + * \param parent + * \param xplore + * \param Glu global LU data + * \return 0 success + * > 0 number of bytes allocated when run out of space + * + */ +int SparseLU::LU_column_dfs(const int m, const int jcol, VectorXi& perm_r, VectorXi& nseg VectorXi& lsub_col, VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, VectorXi& marker, VectorXi& parent, VectorXi& xplore, LU_GlobalLu_t& Glu) +{ + typedef typename VectorXi::Index; + + int jcolp1, jcolm1, jsuper, nsuper, nextl; + int krow; // Row index of the current element + int kperm; // permuted row index + int krep; // Supernode reprentative of the current row + int k, kmark; + int chperm, chmark, chrep, oldrep, kchild; + int myfnz; // First nonzero element in the current column + int xdfs, maxdfs, kpar; + + // Initialize pointers + VectorXi& xsup = Glu.xsup; + VectorXi& supno = Glu.supno; + VectorXi& lsub = Glu.lsub; + VectorXi& xlsub = Glu.xlsub; + + nsuper = supno(jcol); + jsuper = nsuper; + nextl = xlsup(jcol); + VectorBlock<VectorXi> marker2(marker, 2*m, m); + // For each nonzero in A(*,jcol) do dfs + for (k = 0; lsub_col[k] != IND_EMPTY; k++) + { + krow = lsub_col(k); + lsub_col(k) = IND_EMPTY; + kmark = marker2(krow); + + // krow was visited before, go to the next nonz; + if (kmark == jcol) continue; + + // For each unmarker nbr krow of jcol + // krow is in L: place it in structure of L(*,jcol) + marker2(krow) = jcol; + kperm = perm_r(krow); + + if (kperm == IND_EMPTY ) + { + lsub(nextl++) = krow; // krow is indexed into A + if ( nextl >= nzlmax ) + { + Glu.lsub = LUMemXpand<Index>(jcol, nextl, LSUB, nzlmax); + //FIXME try... catch out of space + Glu.nzlmax = nzlmax; + lsub = Glu.lsub; + } + if (kmark != jcolm1) jsuper = IND_EMPTY; // Row index subset testing + } + else + { + // krow is in U : if its supernode-rep krep + // has been explored, update repfnz(*) + krep = xsup(supno(kperm)+1) - 1; + myfnz = repfnz(krep); + + if (myfnz != IND_EMPTY ) + { + // visited before + if (myfnz > kperm) repfnz(krep) = kperm; + // continue; + } + else + { + // otherwise, perform dfs starting at krep + oldrep = IND_EMPTY; + parent(krep) = oldrep; + repfnz(krep) = kperm; + xdfs = xlsub(krep); + maxdfs = xprune(krep); + + do + { + // For each unmarked kchild of krep + while (xdfs < maxdfs) + { + kchild = lsub(xdfs); + xdfs++; + chmark = marker2(kchild); + + if (chmark != jcol) + { + // Not reached yet + marker2(kchild) = jcol; + chperm = perm_r(kchild); + + // if kchild is in L: place it in L(*,k) + if (chperm == IND_EMPTY) + { + lsub(nextl++) = kchild; + if (nextl >= nzlmax) + { + Glu.lsub = LUMemXpand<Index>(jcol, nextl, LSUB, nzlmax); + //FIXME Catch out of space errors + GLu.nzlmax = nzlmax; + lsub = Glu.lsub; + } + if (chmark != jcolm1) jsuper = IND_EMPTY; + } + else + { + // if 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(chrep); + if (myfnz != IND_EMPTY) + { + // Visited before + if ( myfnz > chperm) repfnz(chrep) = chperm; + } + else + { + // continue dfs at super-rep of kchild + xplore(krep) = xdfs; + oldrep = krep; + krep = chrep; // Go deeped down G(L^t) + parent(krep) = olddrep; + repfnz(krep) = chperm; + xdfs = xlsub(krep); + maxdfs = xprune(krep); + } // else myfnz + } // else for chperm + + } // if chmark + + } // end while + + // krow has no more unexplored nbrs; + // place supernode-rep krep in postorder DFS. + // backtrack dfs to its parent + + segrep(nseg) = ;krep; + ++nseg; + kpar = parent(krep); // Pop from stack, mimic recursion + if (kpar == IND_EMPTY) break; // dfs done + krep = kpar; + xdfs = xplore(krep); + maxdfs = xprune(krep); + + } while ( kpar != IND_EMPTY); + + } // else myfnz + + } // else kperm + + } // for each nonzero ... + + // check to see if j belongs in the same supeprnode as j-1 + if ( jcol == 0 ) + { // Do nothing for column 0 + nsuper = supno(0) = 0 ; + } + else + { + fsupc = xsup(nsuper); + jptr = xlsub(jcol); // Not yet compressed + jm1ptr = xlsub(jcolm1); + + // Make sure the number of columns in a supernode doesn't + // exceed threshold + if ( (jcol - fsupc) >= m_maxsuper) jsuper = IND_EMPTY; + + /* If jcol starts a new supernode, reclaim storage space in + * lsub from previous supernode. Note we only store + * the subscript set of the first and last columns of + * a supernode. (first for num values, last for pruning) + */ + if (jsuper == IND_EMPTY) + { // starts a new supernode + if ( (fsupc < jcolm1-1) ) + { // >= 3 columns in nsuper + ito = xlsub(fsupcc+1) + xlsub(jcolm1) = ito; + istop = ito + jptr - jm1ptr; + xprune(jcolm1) = istop; // intialize xprune(jcol-1) + xlsub(jcol) = istop; + + for (ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito) + lsub(ito) = lsub(ifrom); + nextl = ito; // = istop + length(jcol) + } + nsuper++; + supno(jcol) = nsuper; + } // if a new supernode + } // end else: jcol > 0 + + // Tidy up the pointers before exit + xsup(nsuper+1) = jcolp1; + supno(jcolp1) = nsuper; + xprune(jcol) = nextl; // Intialize upper bound for pruning + xlsub(jcolp1) = nextl; + + return 0; +} +#endif
\ No newline at end of file diff --git a/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/Eigen/src/SparseLU/SparseLU_panel_bmod.h index 29cc6d0f0..93daa938c 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h @@ -24,7 +24,7 @@ /* - * NOTE: This file is the modified version of xpanel_dfs.c file in SuperLU + * NOTE: This file is the modified version of xpanel_bmod.c file in SuperLU * -- SuperLU routine (version 3.0) -- * Univ. of California Berkeley, Xerox Palo Alto Research Center, @@ -111,7 +111,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int VectorBLock<VectorXi> dense_col(dense.segment(nextl_col, m)); // Scatter/gather entire matrix column from/to here kfnz = repfnz_col(krep); - if ( kfnz == -1 ) + if ( kfnz == IND_EMPTY ) continue; // skip any zero segment segsize = krep - kfnz + 1; @@ -143,7 +143,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int 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) ); + new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); Map<VectorType> l( &(tempv.data()[segsize]), segsize); l= A * u; @@ -157,7 +157,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int { irow = lsub(isub); dense_col(irow) = tempv(i); - tempv(i) = zero; + tempv(i) = Scalar(0.0); isub++; } diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h index 7b85b6d7c..97e5121db 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h @@ -72,7 +72,7 @@ * */ template <typename MatrixType, typename VectorType> -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) +void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, VectorXi& perm_r, 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 @@ -115,7 +115,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType // For each unmarked krow of jj marker(krow) = jj; kperm = perm_r(krow); - if (kperm == -1 ) { + if (kperm == IND_EMPTY ) { // krow is in L : place it in structure of L(*, jj) panel_lsub(nextl_col++) = krow; // krow is indexed into A } @@ -126,7 +126,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType krep = xsup(supno(kperm)+1) - 1; myfnz = repfnz_col(krep); - if (myfnz != -1 ) + if (myfnz != IND_EMPTY ) { // Representative visited before if (myfnz > kperm ) repfnz_col(krep) = kperm; @@ -135,7 +135,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType else { // Otherwise, perform dfs starting at krep - oldrep = -1; + oldrep = IND_EMPTY; parent(krep) = oldrep; repfnz_col(krep) = kperm; xdfs = xlsub(krep); @@ -155,7 +155,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType marker(kchild) = jj; chperm = perm_r(kchild); - if (chperm == -1) + if (chperm == IND_EMPTY) { // case kchild is in L: place it in L(*, j) panel_lsub(nextl_col++) = kchild; @@ -168,7 +168,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType chrep = xsup(supno(chperm)+1) - 1; myfnz = repfnz_col(chrep); - if (myfnz != -1) + if (myfnz != IND_EMPTY) { // Visited before if (myfnz > chperm) repfnz_col(chrep) = chperm; @@ -202,13 +202,13 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType } kpar = parent(krep); // Pop recursion, mimic recursion - if (kpar == -1) + if (kpar == IND_EMPTY) break; // dfs done krep = kpar; xdfs = xplore(krep); maxdfs = xprune(krep); - } while (kpar != -1); // Do until empty stack + } while (kpar != IND_EMPTY); // Do until empty stack } // end if (myfnz = -1) diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h index e7146a262..9da986497 100644 --- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h @@ -54,9 +54,9 @@ int SparseLU::LU_dsnode_bmod (const int jcol, const int jsupno, const int fsupc, 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; + int irow, isub; // Process the supernodal portion of L\U[*,jcol] - for (int isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++) + for (isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++) { irow = lsub(isub); lusup(nextlu) = dense(irow); @@ -72,20 +72,16 @@ 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); // 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.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); - Map<Matrix<Scalar,Dynamic,1> > l(&(lusup.data()[ufirst]), nsupc); - l = A.triangularView<Lower>().solve(l); + Map<Matrix<Scalar,Dynamic,Dynamic>,0,OuterStride<> > A( &(lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); + Map<Matrix<Scalar,Dynamic,1> > u(&(lusup.data()[ufirst]), nsupc); + u = A.triangularView<Lower>().solve(u); // Update the trailing part of the column jcol U(jcol:jcol+nrow, jcol) using L(jcol:jcol+nrow, fsupc:jcol) and U(fsupc:jcol) - 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); + new (&A) Map<Matrix<Scalar,Dynamic,Dynamic>,0,OuterStride<> > ( &(lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); + Map<Matrix<Scalar,Dynamic,1> > l(&(lusup.data()[ufirst+nsupc], nsupc); + l = l - A * u; return 0; } |