![]() |
Eigen
3.3.3
|
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 StorageIndex; 00041 panel_dfs_traits(Index jcol, StorageIndex* marker) 00042 : m_jcol(jcol), m_marker(marker) 00043 {} 00044 bool update_segrep(Index krep, StorageIndex 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 StorageIndex* m_marker; 00057 }; 00058 00059 00060 template <typename Scalar, typename StorageIndex> 00061 template <typename Traits> 00062 void SparseLUImpl<Scalar,StorageIndex>::dfs_kernel(const StorageIndex 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 StorageIndex kmark = marker(krow); 00071 00072 // For each unmarked krow of jj 00073 marker(krow) = jj; 00074 StorageIndex 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++) = StorageIndex(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 StorageIndex krep = glu.xsup(glu.supno(kperm)+1) - 1; 00087 // First nonzero element in the current column: 00088 StorageIndex 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 StorageIndex oldrep = emptyIdxLU; 00100 parent(krep) = oldrep; 00101 repfnz_col(krep) = kperm; 00102 StorageIndex xdfs = glu.xlsub(krep); 00103 Index maxdfs = xprune(krep); 00104 00105 StorageIndex kpar; 00106 do 00107 { 00108 // For each unmarked kchild of krep 00109 while (xdfs < maxdfs) 00110 { 00111 StorageIndex kchild = glu.lsub(xdfs); 00112 xdfs++; 00113 StorageIndex chmark = marker(kchild); 00114 00115 if (chmark != jj ) 00116 { 00117 marker(kchild) = jj; 00118 StorageIndex 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 StorageIndex 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 StorageIndex> 00219 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) 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 (StorageIndex jj = StorageIndex(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 StorageIndex 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