src/SparseLU/SparseLU_panel_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]panel_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_PANEL_DFS_H
0031 #define SPARSELU_PANEL_DFS_H
0032
0033 namespace RivetEigen {
0034
0035 namespace internal {
0036
0037 template<typename IndexVector>
0038 struct panel_dfs_traits
0039 {
0040   typedef typename IndexVector::Scalar StorageIndex;
0041   panel_dfs_traits(Index jcol, StorageIndex* marker)
0042     : m_jcol(jcol), m_marker(marker)
0043   {}
0044   bool update_segrep(Index krep, StorageIndex jj)
0045   {
0046     if(m_marker[krep]<m_jcol)
0047     {
0048       m_marker[krep] = jj;
0049       return true;
0050     }
0051     return false;
0052   }
0053   void mem_expand(IndexVector& /*glu.lsub*/, Index /*nextl*/, Index /*chmark*/) {}
0054   enum { ExpandMem = false };
0055   Index m_jcol;
0056   StorageIndex* m_marker;
0057 };
0058
0059
0060 template <typename Scalar, typename StorageIndex>
0061 template <typename Traits>
0062 void SparseLUImpl<Scalar,StorageIndex>::dfs_kernel(const StorageIndex jj, IndexVector& perm_r,
0063                    Index& nseg, IndexVector& panel_lsub, IndexVector& segrep,
0064                    Ref<IndexVector> repfnz_col, IndexVector& xprune, Ref<IndexVector> marker, IndexVector& parent,
0065                    IndexVector& xplore, GlobalLU_t& glu,
0066                    Index& nextl_col, Index krow, Traits& traits
0067                   )
0068 {
0069
0070   StorageIndex kmark = marker(krow);
0071
0072   // For each unmarked krow of jj
0073   marker(krow) = jj;
0074   StorageIndex kperm = perm_r(krow);
0075   if (kperm == emptyIdxLU ) {
0076     // krow is in L : place it in structure of L(*, jj)
0077     panel_lsub(nextl_col++) = StorageIndex(krow);  // krow is indexed into A
0078
0079     traits.mem_expand(panel_lsub, nextl_col, kmark);
0080   }
0081   else
0082   {
0083     // krow is in U : if its supernode-representative krep
0084     // has been explored, update repfnz(*)
0085     // krep = supernode representative of the current row
0086     StorageIndex krep = glu.xsup(glu.supno(kperm)+1) - 1;
0087     // First nonzero element in the current column:
0088     StorageIndex myfnz = repfnz_col(krep);
0089
0090     if (myfnz != emptyIdxLU )
0091     {
0092       // Representative visited before
0093       if (myfnz > kperm ) repfnz_col(krep) = kperm;
0094
0095     }
0096     else
0097     {
0098       // Otherwise, perform dfs starting at krep
0099       StorageIndex oldrep = emptyIdxLU;
0100       parent(krep) = oldrep;
0101       repfnz_col(krep) = kperm;
0102       StorageIndex xdfs =  glu.xlsub(krep);
0103       Index maxdfs = xprune(krep);
0104
0105       StorageIndex kpar;
0106       do
0107       {
0108         // For each unmarked kchild of krep
0109         while (xdfs < maxdfs)
0110         {
0111           StorageIndex kchild = glu.lsub(xdfs);
0112           xdfs++;
0113           StorageIndex chmark = marker(kchild);
0114
0115           if (chmark != jj )
0116           {
0117             marker(kchild) = jj;
0118             StorageIndex chperm = perm_r(kchild);
0119
0120             if (chperm == emptyIdxLU)
0121             {
0122               // case kchild is in L: place it in L(*, j)
0123               panel_lsub(nextl_col++) = kchild;
0124               traits.mem_expand(panel_lsub, nextl_col, chmark);
0125             }
0126             else
0127             {
0128               // case kchild is in U :
0129               // chrep = its supernode-rep. If its rep has been explored,
0130               // update its repfnz(*)
0131               StorageIndex chrep = glu.xsup(glu.supno(chperm)+1) - 1;
0132               myfnz = repfnz_col(chrep);
0133
0134               if (myfnz != emptyIdxLU)
0135               { // Visited before
0136                 if (myfnz > chperm)
0137                   repfnz_col(chrep) = chperm;
0138               }
0139               else
0140               { // Cont. dfs at snode-rep of kchild
0141                 xplore(krep) = xdfs;
0142                 oldrep = krep;
0143                 krep = chrep; // Go deeper down G(L)
0144                 parent(krep) = oldrep;
0145                 repfnz_col(krep) = chperm;
0146                 xdfs = glu.xlsub(krep);
0147                 maxdfs = xprune(krep);
0148
0149               } // end if myfnz != -1
0150             } // end if chperm == -1
0151
0152           } // end if chmark !=jj
0153         } // end while xdfs < maxdfs
0154
0155         // krow has no more unexplored nbrs :
0156         //    Place snode-rep krep in postorder DFS, if this
0157         //    segment is seen for the first time. (Note that
0158         //    "repfnz(krep)" may change later.)
0159         //    Baktrack dfs to its parent
0160         if(traits.update_segrep(krep,jj))
0161         //if (marker1(krep) < jcol )
0162         {
0163           segrep(nseg) = krep;
0164           ++nseg;
0165           //marker1(krep) = jj;
0166         }
0167
0168         kpar = parent(krep); // Pop recursion, mimic recursion
0169         if (kpar == emptyIdxLU)
0170           break; // dfs done
0171         krep = kpar;
0172         xdfs = xplore(krep);
0173         maxdfs = xprune(krep);
0174
0175       } while (kpar != emptyIdxLU); // Do until empty stack
0176
0177     } // end if (myfnz = -1)
0178
0179   } // end if (kperm == -1)
0180 }
0181
0182 /**
0183  * \brief Performs a symbolic factorization on a panel of columns [jcol, jcol+w)
0184  *
0185  * A supernode representative is the last column of a supernode.
0186  * The nonzeros in U[*,j] are segments that end at supernodes representatives
0187  *
0188  * The routine returns a list of the supernodal representatives
0189  * in topological order of the dfs that generates them. This list is
0190  * a superset of the topological order of each individual column within
0191  * the panel.
0192  * The location of the first nonzero in each supernodal segment
0193  * (supernodal entry location) is also returned. Each column has
0194  * a separate list for this purpose.
0195  *
0196  * Two markers arrays are used for dfs :
0197  *    marker[i] == jj, if i was visited during dfs of current column jj;
0198  *    marker1[i] >= jcol, if i was visited by earlier columns in this panel;
0199  *
0200  * \param[in] m number of rows in the matrix
0201  * \param[in] w Panel size
0202  * \param[in] jcol Starting  column of the panel
0203  * \param[in] A Input matrix in column-major storage
0204  * \param[in] perm_r Row permutation
0205  * \param[out] nseg Number of U segments
0206  * \param[out] dense Accumulate the column vectors of the panel
0207  * \param[out] panel_lsub Subscripts of the row in the panel
0208  * \param[out] segrep Segment representative i.e first nonzero row of each segment
0209  * \param[out] repfnz First nonzero location in each row
0210  * \param[out] xprune The pruned elimination tree
0211  * \param[out] marker work vector
0212  * \param  parent The elimination tree
0213  * \param xplore work vector
0214  * \param glu The global data structure
0215  *
0216  */
0217
0218 template <typename Scalar, typename StorageIndex>
0219 void SparseLUImpl<Scalar,StorageIndex>::panel_dfs(const Index m, const Index w, const Index jcol, MatrixType& A, IndexVector& perm_r, Index& nseg, ScalarVector& dense, IndexVector& panel_lsub, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
0220 {
0221   Index nextl_col; // Next available position in panel_lsub[*,jj]
0222
0223   // Initialize pointers
0224   VectorBlock<IndexVector> marker1(marker, m, m);
0225   nseg = 0;
0226
0227   panel_dfs_traits<IndexVector> traits(jcol, marker1.data());
0228
0229   // For each column in the panel
0230   for (StorageIndex jj = StorageIndex(jcol); jj < jcol + w; jj++)
0231   {
0232     nextl_col = (jj - jcol) * m;
0233
0234     VectorBlock<IndexVector> repfnz_col(repfnz, nextl_col, m); // First nonzero location in each row
0235     VectorBlock<ScalarVector> dense_col(dense,nextl_col, m); // Accumulate a column vector here
0236
0237
0238     // For each nnz in A[*, jj] do depth first search
0239     for (typename MatrixType::InnerIterator it(A, jj); it; ++it)
0240     {
0241       Index krow = it.row();
0242       dense_col(krow) = it.value();
0243
0244       StorageIndex kmark = marker(krow);
0245       if (kmark == jj)
0246         continue; // krow visited before, go to the next nonzero
0247
0248       dfs_kernel(jj, perm_r, nseg, panel_lsub, segrep, repfnz_col, xprune, marker, parent,
0249                    xplore, glu, nextl_col, krow, traits);
0250     }// end for nonzeros in column jj
0251
0252   } // end for column jj
0253 }
0254
0255 } // end namespace internal
0256 } // end namespace RivetEigen
0257
0258 #endif // SPARSELU_PANEL_DFS_H