Eigen  3.3.3
SparseLU_column_dfs.h
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]column_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_COLUMN_DFS_H
00031 #define SPARSELU_COLUMN_DFS_H
00032 
00033 template <typename Scalar, typename StorageIndex> class SparseLUImpl;
00034 namespace Eigen {
00035 
00036 namespace internal {
00037 
00038 template<typename IndexVector, typename ScalarVector>
00039 struct column_dfs_traits : no_assignment_operator
00040 {
00041   typedef typename ScalarVector::Scalar Scalar;
00042   typedef typename IndexVector::Scalar StorageIndex;
00043   column_dfs_traits(Index jcol, Index& jsuper, typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& glu, SparseLUImpl<Scalar, StorageIndex>& luImpl)
00044    : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu), m_luImpl(luImpl)
00045  {}
00046   bool update_segrep(Index /*krep*/, Index /*jj*/)
00047   {
00048     return true;
00049   }
00050   void mem_expand(IndexVector& lsub, Index& nextl, Index chmark)
00051   {
00052     if (nextl >= m_glu.nzlmax)
00053       m_luImpl.memXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions); 
00054     if (chmark != (m_jcol-1)) m_jsuper_ref = emptyIdxLU;
00055   }
00056   enum { ExpandMem = true };
00057   
00058   Index m_jcol;
00059   Index& m_jsuper_ref;
00060   typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& m_glu;
00061   SparseLUImpl<Scalar, StorageIndex>& m_luImpl;
00062 };
00063 
00064 
00092 template <typename Scalar, typename StorageIndex>
00093 Index SparseLUImpl<Scalar,StorageIndex>::column_dfs(const Index m, const Index jcol, IndexVector& perm_r, Index maxsuper, Index& nseg,
00094                                                     BlockIndexVector lsub_col, IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune,
00095                                                     IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
00096 {
00097   
00098   Index jsuper = glu.supno(jcol); 
00099   Index nextl = glu.xlsub(jcol); 
00100   VectorBlock<IndexVector> marker2(marker, 2*m, m); 
00101   
00102   
00103   column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu, *this);
00104   
00105   // For each nonzero in A(*,jcol) do dfs 
00106   for (Index k = 0; ((k < m) ? lsub_col[k] != emptyIdxLU : false) ; k++)
00107   {
00108     Index krow = lsub_col(k); 
00109     lsub_col(k) = emptyIdxLU; 
00110     Index kmark = marker2(krow); 
00111     
00112     // krow was visited before, go to the next nonz; 
00113     if (kmark == jcol) continue;
00114     
00115     dfs_kernel(StorageIndex(jcol), perm_r, nseg, glu.lsub, segrep, repfnz, xprune, marker2, parent,
00116                    xplore, glu, nextl, krow, traits);
00117   } // for each nonzero ... 
00118   
00119   Index fsupc;
00120   StorageIndex nsuper = glu.supno(jcol);
00121   StorageIndex jcolp1 = StorageIndex(jcol) + 1;
00122   Index jcolm1 = jcol - 1;
00123   
00124   // check to see if j belongs in the same supernode as j-1
00125   if ( jcol == 0 )
00126   { // Do nothing for column 0 
00127     nsuper = glu.supno(0) = 0 ;
00128   }
00129   else 
00130   {
00131     fsupc = glu.xsup(nsuper); 
00132     StorageIndex jptr = glu.xlsub(jcol); // Not yet compressed
00133     StorageIndex jm1ptr = glu.xlsub(jcolm1); 
00134     
00135     // Use supernodes of type T2 : see SuperLU paper
00136     if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = emptyIdxLU;
00137     
00138     // Make sure the number of columns in a supernode doesn't
00139     // exceed threshold
00140     if ( (jcol - fsupc) >= maxsuper) jsuper = emptyIdxLU; 
00141     
00142     /* If jcol starts a new supernode, reclaim storage space in
00143      * glu.lsub from previous supernode. Note we only store 
00144      * the subscript set of the first and last columns of 
00145      * a supernode. (first for num values, last for pruning)
00146      */
00147     if (jsuper == emptyIdxLU)
00148     { // starts a new supernode 
00149       if ( (fsupc < jcolm1-1) ) 
00150       { // >= 3 columns in nsuper
00151         StorageIndex ito = glu.xlsub(fsupc+1);
00152         glu.xlsub(jcolm1) = ito; 
00153         StorageIndex istop = ito + jptr - jm1ptr; 
00154         xprune(jcolm1) = istop; // intialize xprune(jcol-1)
00155         glu.xlsub(jcol) = istop; 
00156         
00157         for (StorageIndex ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
00158           glu.lsub(ito) = glu.lsub(ifrom); 
00159         nextl = ito;  // = istop + length(jcol)
00160       }
00161       nsuper++; 
00162       glu.supno(jcol) = nsuper; 
00163     } // if a new supernode 
00164   } // end else:  jcol > 0
00165   
00166   // Tidy up the pointers before exit
00167   glu.xsup(nsuper+1) = jcolp1; 
00168   glu.supno(jcolp1) = nsuper; 
00169   xprune(jcol) = StorageIndex(nextl);  // Intialize upper bound for pruning
00170   glu.xlsub(jcolp1) = StorageIndex(nextl); 
00171   
00172   return 0; 
00173 }
00174 
00175 } // end namespace internal
00176 
00177 } // end namespace Eigen
00178 
00179 #endif
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends