SparseLU_panel_dfs.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 /* 
00011  
00012  * NOTE: This file is the modified version of [s,d,c,z]panel_dfs.c file in SuperLU 
00013  
00014  * -- SuperLU routine (version 2.0) --
00015  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
00016  * and Lawrence Berkeley National Lab.
00017  * November 15, 1997
00018  *
00019  * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
00020  *
00021  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
00022  * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
00023  *
00024  * Permission is hereby granted to use or copy this program for any
00025  * purpose, provided the above notices are retained on all copies.
00026  * Permission to modify the code and to distribute modified code is
00027  * granted, provided the above notices are retained, and a notice that
00028  * the code was modified is included with the above copyright notice.
00029  */
00030 #ifndef SPARSELU_PANEL_DFS_H
00031 #define SPARSELU_PANEL_DFS_H
00032 
00033 namespace Eigen {
00034 
00035 namespace internal {
00036   
00037 template<typename IndexVector>
00038 struct panel_dfs_traits
00039 {
00040   typedef typename IndexVector::Scalar Index;
00041   panel_dfs_traits(Index jcol, Index* marker)
00042     : m_jcol(jcol), m_marker(marker)
00043   {}
00044   bool update_segrep(Index krep, Index jj)
00045   {
00046     if(m_marker[krep]<m_jcol)
00047     {
00048       m_marker[krep] = jj; 
00049       return true;
00050     }
00051     return false;
00052   }
00053   void mem_expand(IndexVector& /*glu.lsub*/, Index /*nextl*/, Index /*chmark*/) {}
00054   enum { ExpandMem = false };
00055   Index m_jcol;
00056   Index* m_marker;
00057 };
00058 
00059 
00060 template <typename Scalar, typename Index>
00061 template <typename Traits>
00062 void SparseLUImpl<Scalar,Index>::dfs_kernel(const Index jj, IndexVector& perm_r,
00063                    Index& nseg, IndexVector& panel_lsub, IndexVector& segrep,
00064                    Ref<IndexVector> repfnz_col, IndexVector& xprune, Ref<IndexVector> marker, IndexVector& parent,
00065                    IndexVector& xplore, GlobalLU_t& glu,
00066                    Index& nextl_col, Index krow, Traits& traits
00067                   )
00068 {
00069   
00070   Index kmark = marker(krow);
00071       
00072   // For each unmarked krow of jj
00073   marker(krow) = jj; 
00074   Index kperm = perm_r(krow); 
00075   if (kperm == emptyIdxLU ) {
00076     // krow is in L : place it in structure of L(*, jj)
00077     panel_lsub(nextl_col++) = krow;  // krow is indexed into A
00078     
00079     traits.mem_expand(panel_lsub, nextl_col, kmark);
00080   }
00081   else 
00082   {
00083     // krow is in U : if its supernode-representative krep
00084     // has been explored, update repfnz(*)
00085     // krep = supernode representative of the current row
00086     Index krep = glu.xsup(glu.supno(kperm)+1) - 1; 
00087     // First nonzero element in the current column:
00088     Index myfnz = repfnz_col(krep); 
00089     
00090     if (myfnz != emptyIdxLU )
00091     {
00092       // Representative visited before
00093       if (myfnz > kperm ) repfnz_col(krep) = kperm; 
00094       
00095     }
00096     else 
00097     {
00098       // Otherwise, perform dfs starting at krep
00099       Index oldrep = emptyIdxLU; 
00100       parent(krep) = oldrep; 
00101       repfnz_col(krep) = kperm; 
00102       Index xdfs =  glu.xlsub(krep); 
00103       Index maxdfs = xprune(krep); 
00104       
00105       Index kpar;
00106       do 
00107       {
00108         // For each unmarked kchild of krep
00109         while (xdfs < maxdfs) 
00110         {
00111           Index kchild = glu.lsub(xdfs); 
00112           xdfs++; 
00113           Index chmark = marker(kchild); 
00114           
00115           if (chmark != jj ) 
00116           {
00117             marker(kchild) = jj; 
00118             Index chperm = perm_r(kchild); 
00119             
00120             if (chperm == emptyIdxLU) 
00121             {
00122               // case kchild is in L: place it in L(*, j)
00123               panel_lsub(nextl_col++) = kchild;
00124               traits.mem_expand(panel_lsub, nextl_col, chmark);
00125             }
00126             else
00127             {
00128               // case kchild is in U :
00129               // chrep = its supernode-rep. If its rep has been explored, 
00130               // update its repfnz(*)
00131               Index chrep = glu.xsup(glu.supno(chperm)+1) - 1; 
00132               myfnz = repfnz_col(chrep); 
00133               
00134               if (myfnz != emptyIdxLU) 
00135               { // Visited before 
00136                 if (myfnz > chperm) 
00137                   repfnz_col(chrep) = chperm; 
00138               }
00139               else 
00140               { // Cont. dfs at snode-rep of kchild
00141                 xplore(krep) = xdfs; 
00142                 oldrep = krep; 
00143                 krep = chrep; // Go deeper down G(L)
00144                 parent(krep) = oldrep; 
00145                 repfnz_col(krep) = chperm; 
00146                 xdfs = glu.xlsub(krep); 
00147                 maxdfs = xprune(krep); 
00148                 
00149               } // end if myfnz != -1
00150             } // end if chperm == -1 
00151                 
00152           } // end if chmark !=jj
00153         } // end while xdfs < maxdfs
00154         
00155         // krow has no more unexplored nbrs :
00156         //    Place snode-rep krep in postorder DFS, if this 
00157         //    segment is seen for the first time. (Note that 
00158         //    "repfnz(krep)" may change later.)
00159         //    Baktrack dfs to its parent
00160         if(traits.update_segrep(krep,jj))
00161         //if (marker1(krep) < jcol )
00162         {
00163           segrep(nseg) = krep; 
00164           ++nseg; 
00165           //marker1(krep) = jj; 
00166         }
00167         
00168         kpar = parent(krep); // Pop recursion, mimic recursion 
00169         if (kpar == emptyIdxLU) 
00170           break; // dfs done 
00171         krep = kpar; 
00172         xdfs = xplore(krep); 
00173         maxdfs = xprune(krep); 
00174 
00175       } while (kpar != emptyIdxLU); // Do until empty stack 
00176       
00177     } // end if (myfnz = -1)
00178 
00179   } // end if (kperm == -1)   
00180 }
00181 
00218 template <typename Scalar, typename Index>
00219 void SparseLUImpl<Scalar,Index>::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)
00220 {
00221   Index nextl_col; // Next available position in panel_lsub[*,jj] 
00222   
00223   // Initialize pointers 
00224   VectorBlock<IndexVector> marker1(marker, m, m); 
00225   nseg = 0; 
00226   
00227   panel_dfs_traits<IndexVector> traits(jcol, marker1.data());
00228   
00229   // For each column in the panel 
00230   for (Index jj = jcol; jj < jcol + w; jj++) 
00231   {
00232     nextl_col = (jj - jcol) * m; 
00233     
00234     VectorBlock<IndexVector> repfnz_col(repfnz, nextl_col, m); // First nonzero location in each row
00235     VectorBlock<ScalarVector> dense_col(dense,nextl_col, m); // Accumulate a column vector here
00236     
00237     
00238     // For each nnz in A[*, jj] do depth first search
00239     for (typename MatrixType::InnerIterator it(A, jj); it; ++it)
00240     {
00241       Index krow = it.row(); 
00242       dense_col(krow) = it.value();
00243       
00244       Index kmark = marker(krow); 
00245       if (kmark == jj) 
00246         continue; // krow visited before, go to the next nonzero
00247       
00248       dfs_kernel(jj, perm_r, nseg, panel_lsub, segrep, repfnz_col, xprune, marker, parent,
00249                    xplore, glu, nextl_col, krow, traits);
00250     }// end for nonzeros in column jj
00251     
00252   } // end for column jj
00253 }
00254 
00255 } // end namespace internal
00256 } // end namespace Eigen
00257 
00258 #endif // SPARSELU_PANEL_DFS_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:39:23