![]() |
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]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