SparseLU_column_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]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 Index> 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 Index;
00043   column_dfs_traits(Index jcol, Index& jsuper, typename SparseLUImpl<Scalar, Index>::GlobalLU_t& glu, SparseLUImpl<Scalar, Index>& 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, Index>::GlobalLU_t& m_glu;
00061   SparseLUImpl<Scalar, Index>& m_luImpl;
00062 };
00063 
00064 
00092 template <typename Scalar, typename Index>
00093 Index SparseLUImpl<Scalar,Index>::column_dfs(const Index m, const Index jcol, IndexVector& perm_r, Index maxsuper, Index& nseg,  BlockIndexVector lsub_col, IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
00094 {
00095   
00096   Index jsuper = glu.supno(jcol); 
00097   Index nextl = glu.xlsub(jcol); 
00098   VectorBlock<IndexVector> marker2(marker, 2*m, m); 
00099   
00100   
00101   column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu, *this);
00102   
00103   // For each nonzero in A(*,jcol) do dfs 
00104   for (Index k = 0; ((k < m) ? lsub_col[k] != emptyIdxLU : false) ; k++)
00105   {
00106     Index krow = lsub_col(k); 
00107     lsub_col(k) = emptyIdxLU; 
00108     Index kmark = marker2(krow); 
00109     
00110     // krow was visited before, go to the next nonz; 
00111     if (kmark == jcol) continue;
00112     
00113     dfs_kernel(jcol, perm_r, nseg, glu.lsub, segrep, repfnz, xprune, marker2, parent,
00114                    xplore, glu, nextl, krow, traits);
00115   } // for each nonzero ... 
00116   
00117   Index fsupc, jptr, jm1ptr, ito, ifrom, istop;
00118   Index nsuper = glu.supno(jcol);
00119   Index jcolp1 = jcol + 1;
00120   Index jcolm1 = jcol - 1;
00121   
00122   // check to see if j belongs in the same supernode as j-1
00123   if ( jcol == 0 )
00124   { // Do nothing for column 0 
00125     nsuper = glu.supno(0) = 0 ;
00126   }
00127   else 
00128   {
00129     fsupc = glu.xsup(nsuper); 
00130     jptr = glu.xlsub(jcol); // Not yet compressed
00131     jm1ptr = glu.xlsub(jcolm1); 
00132     
00133     // Use supernodes of type T2 : see SuperLU paper
00134     if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = emptyIdxLU;
00135     
00136     // Make sure the number of columns in a supernode doesn't
00137     // exceed threshold
00138     if ( (jcol - fsupc) >= maxsuper) jsuper = emptyIdxLU; 
00139     
00140     /* If jcol starts a new supernode, reclaim storage space in
00141      * glu.lsub from previous supernode. Note we only store 
00142      * the subscript set of the first and last columns of 
00143      * a supernode. (first for num values, last for pruning)
00144      */
00145     if (jsuper == emptyIdxLU)
00146     { // starts a new supernode 
00147       if ( (fsupc < jcolm1-1) ) 
00148       { // >= 3 columns in nsuper
00149         ito = glu.xlsub(fsupc+1);
00150         glu.xlsub(jcolm1) = ito; 
00151         istop = ito + jptr - jm1ptr; 
00152         xprune(jcolm1) = istop; // intialize xprune(jcol-1)
00153         glu.xlsub(jcol) = istop; 
00154         
00155         for (ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
00156           glu.lsub(ito) = glu.lsub(ifrom); 
00157         nextl = ito;  // = istop + length(jcol)
00158       }
00159       nsuper++; 
00160       glu.supno(jcol) = nsuper; 
00161     } // if a new supernode 
00162   } // end else:  jcol > 0
00163   
00164   // Tidy up the pointers before exit
00165   glu.xsup(nsuper+1) = jcolp1; 
00166   glu.supno(jcolp1) = nsuper; 
00167   xprune(jcol) = nextl;  // Intialize upper bound for pruning
00168   glu.xlsub(jcolp1) = nextl; 
00169   
00170   return 0; 
00171 }
00172 
00173 } // end namespace internal
00174 
00175 } // end namespace Eigen
00176 
00177 #endif


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