src/SparseLU/SparseLU_column_dfs.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 /*
0011
0012  * NOTE: This file is the modified version of [s,d,c,z]column_dfs.c file in SuperLU
0013
0014  * -- SuperLU routine (version 2.0) --
0015  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
0016  * and Lawrence Berkeley National Lab.
0017  * November 15, 1997
0018  *
0019  * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
0020  *
0021  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
0022  * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
0023  *
0024  * Permission is hereby granted to use or copy this program for any
0025  * purpose, provided the above notices are retained on all copies.
0026  * Permission to modify the code and to distribute modified code is
0027  * granted, provided the above notices are retained, and a notice that
0028  * the code was modified is included with the above copyright notice.
0029  */
0030 #ifndef SPARSELU_COLUMN_DFS_H
0031 #define SPARSELU_COLUMN_DFS_H
0032
0033 template <typename Scalar, typename StorageIndex> class SparseLUImpl;
0034 namespace RivetEigen {
0035
0036 namespace internal {
0037
0038 template<typename IndexVector, typename ScalarVector>
0039 struct column_dfs_traits : no_assignment_operator
0040 {
0041   typedef typename ScalarVector::Scalar Scalar;
0042   typedef typename IndexVector::Scalar StorageIndex;
0043   column_dfs_traits(Index jcol, Index& jsuper, typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& glu, SparseLUImpl<Scalar, StorageIndex>& luImpl)
0044    : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu), m_luImpl(luImpl)
0045  {}
0046   bool update_segrep(Index /*krep*/, Index /*jj*/)
0047   {
0048     return true;
0049   }
0050   void mem_expand(IndexVector& lsub, Index& nextl, Index chmark)
0051   {
0052     if (nextl >= m_glu.nzlmax)
0053       m_luImpl.memXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions);
0054     if (chmark != (m_jcol-1)) m_jsuper_ref = emptyIdxLU;
0055   }
0056   enum { ExpandMem = true };
0057
0058   Index m_jcol;
0059   Index& m_jsuper_ref;
0060   typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& m_glu;
0061   SparseLUImpl<Scalar, StorageIndex>& m_luImpl;
0062 };
0063
0064
0065 /**
0066  * \brief Performs a symbolic factorization on column jcol and decide the supernode boundary
0067  *
0068  * A supernode representative is the last column of a supernode.
0069  * The nonzeros in U[*,j] are segments that end at supernodes representatives.
0070  * The routine returns a list of the supernodal representatives
0071  * in topological order of the dfs that generates them.
0072  * The location of the first nonzero in each supernodal segment
0073  * (supernodal entry location) is also returned.
0074  *
0075  * \param m number of rows in the matrix
0076  * \param jcol Current column
0077  * \param perm_r Row permutation
0078  * \param maxsuper  Maximum number of column allowed in a supernode
0079  * \param [in,out] nseg Number of segments in current U[*,j] - new segments appended
0080  * \param lsub_col defines the rhs vector to start the dfs
0081  * \param [in,out] segrep Segment representatives - new segments appended
0082  * \param repfnz  First nonzero location in each row
0083  * \param xprune
0084  * \param marker  marker[i] == jj, if i was visited during dfs of current column jj;
0085  * \param parent
0086  * \param xplore working array
0087  * \param glu global LU data
0088  * \return 0 success
0089  *         > 0 number of bytes allocated when run out of space
0090  *
0091  */
0092 template <typename Scalar, typename StorageIndex>
0093 Index SparseLUImpl<Scalar,StorageIndex>::column_dfs(const Index m, const Index jcol, IndexVector& perm_r, Index maxsuper, Index& nseg,
0094                                                     BlockIndexVector lsub_col, IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune,
0095                                                     IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
0096 {
0097
0098   Index jsuper = glu.supno(jcol);
0099   Index nextl = glu.xlsub(jcol);
0100   VectorBlock<IndexVector> marker2(marker, 2*m, m);
0101
0102
0103   column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu, *this);
0104
0105   // For each nonzero in A(*,jcol) do dfs
0106   for (Index k = 0; ((k < m) ? lsub_col[k] != emptyIdxLU : false) ; k++)
0107   {
0108     Index krow = lsub_col(k);
0109     lsub_col(k) = emptyIdxLU;
0110     Index kmark = marker2(krow);
0111
0112     // krow was visited before, go to the next nonz;
0113     if (kmark == jcol) continue;
0114
0115     dfs_kernel(StorageIndex(jcol), perm_r, nseg, glu.lsub, segrep, repfnz, xprune, marker2, parent,
0116                    xplore, glu, nextl, krow, traits);
0117   } // for each nonzero ...
0118
0119   Index fsupc;
0120   StorageIndex nsuper = glu.supno(jcol);
0121   StorageIndex jcolp1 = StorageIndex(jcol) + 1;
0122   Index jcolm1 = jcol - 1;
0123
0124   // check to see if j belongs in the same supernode as j-1
0125   if ( jcol == 0 )
0126   { // Do nothing for column 0
0127     nsuper = glu.supno(0) = 0 ;
0128   }
0129   else
0130   {
0131     fsupc = glu.xsup(nsuper);
0132     StorageIndex jptr = glu.xlsub(jcol); // Not yet compressed
0133     StorageIndex jm1ptr = glu.xlsub(jcolm1);
0134
0135     // Use supernodes of type T2 : see SuperLU paper
0136     if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = emptyIdxLU;
0137
0138     // Make sure the number of columns in a supernode doesn't
0139     // exceed threshold
0140     if ( (jcol - fsupc) >= maxsuper) jsuper = emptyIdxLU;
0141
0142     /* If jcol starts a new supernode, reclaim storage space in
0143      * glu.lsub from previous supernode. Note we only store
0144      * the subscript set of the first and last columns of
0145      * a supernode. (first for num values, last for pruning)
0146      */
0147     if (jsuper == emptyIdxLU)
0148     { // starts a new supernode
0149       if ( (fsupc < jcolm1-1) )
0150       { // >= 3 columns in nsuper
0151         StorageIndex ito = glu.xlsub(fsupc+1);
0152         glu.xlsub(jcolm1) = ito;
0153         StorageIndex istop = ito + jptr - jm1ptr;
0154         xprune(jcolm1) = istop; // initialize xprune(jcol-1)
0155         glu.xlsub(jcol) = istop;
0156
0157         for (StorageIndex ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
0158           glu.lsub(ito) = glu.lsub(ifrom);
0159         nextl = ito;  // = istop + length(jcol)
0160       }
0161       nsuper++;
0162       glu.supno(jcol) = nsuper;
0163     } // if a new supernode
0164   } // end else:  jcol > 0
0165
0166   // Tidy up the pointers before exit
0167   glu.xsup(nsuper+1) = jcolp1;
0168   glu.supno(jcolp1) = nsuper;
0169   xprune(jcol) = StorageIndex(nextl);  // Initialize upper bound for pruning
0170   glu.xlsub(jcolp1) = StorageIndex(nextl);
0171
0172   return 0;
0173 }
0174
0175 } // end namespace internal
0176
0177 } // end namespace RivetEigen
0178
0179 #endif