1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 /* 11 12 * NOTE: This file is the modified version of [s,d,c,z]column_dfs.c file in SuperLU 13 14 * -- SuperLU routine (version 2.0) -- 15 * Univ. of California Berkeley, Xerox Palo Alto Research Center, 16 * and Lawrence Berkeley National Lab. 17 * November 15, 1997 18 * 19 * Copyright (c) 1994 by Xerox Corporation. All rights reserved. 20 * 21 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 22 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 23 * 24 * Permission is hereby granted to use or copy this program for any 25 * purpose, provided the above notices are retained on all copies. 26 * Permission to modify the code and to distribute modified code is 27 * granted, provided the above notices are retained, and a notice that 28 * the code was modified is included with the above copyright notice. 29 */ 30 #ifndef SPARSELU_COLUMN_DFS_H 31 #define SPARSELU_COLUMN_DFS_H 32 33 template <typename Scalar, typename StorageIndex> class SparseLUImpl; 34 namespace Eigen { 35 36 namespace internal { 37 38 template<typename IndexVector, typename ScalarVector> 39 struct column_dfs_traits : no_assignment_operator 40 { 41 typedef typename ScalarVector::Scalar Scalar; 42 typedef typename IndexVector::Scalar StorageIndex; 43 column_dfs_traits(Index jcol, Index& jsuper, typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& glu, SparseLUImpl<Scalar, StorageIndex>& luImpl) 44 : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu), m_luImpl(luImpl) 45 {} 46 bool update_segrep(Index /*krep*/, Index /*jj*/) 47 { 48 return true; 49 } 50 void mem_expand(IndexVector& lsub, Index& nextl, Index chmark) 51 { 52 if (nextl >= m_glu.nzlmax) 53 m_luImpl.memXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions); 54 if (chmark != (m_jcol-1)) m_jsuper_ref = emptyIdxLU; 55 } 56 enum { ExpandMem = true }; 57 58 Index m_jcol; 59 Index& m_jsuper_ref; 60 typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& m_glu; 61 SparseLUImpl<Scalar, StorageIndex>& m_luImpl; 62 }; 63 64 65 /** 66 * \brief Performs a symbolic factorization on column jcol and decide the supernode boundary 67 * 68 * A supernode representative is the last column of a supernode. 69 * The nonzeros in U[*,j] are segments that end at supernodes representatives. 70 * The routine returns a list of the supernodal representatives 71 * in topological order of the dfs that generates them. 72 * The location of the first nonzero in each supernodal segment 73 * (supernodal entry location) is also returned. 74 * 75 * \param m number of rows in the matrix 76 * \param jcol Current column 77 * \param perm_r Row permutation 78 * \param maxsuper Maximum number of column allowed in a supernode 79 * \param [in,out] nseg Number of segments in current U[*,j] - new segments appended 80 * \param lsub_col defines the rhs vector to start the dfs 81 * \param [in,out] segrep Segment representatives - new segments appended 82 * \param repfnz First nonzero location in each row 83 * \param xprune 84 * \param marker marker[i] == jj, if i was visited during dfs of current column jj; 85 * \param parent 86 * \param xplore working array 87 * \param glu global LU data 88 * \return 0 success 89 * > 0 number of bytes allocated when run out of space 90 * 91 */ 92 template <typename Scalar, typename StorageIndex> 93 Index SparseLUImpl<Scalar,StorageIndex>::column_dfs(const Index m, const Index jcol, IndexVector& perm_r, Index maxsuper, Index& nseg, 94 BlockIndexVector lsub_col, IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, 95 IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu) 96 { 97 98 Index jsuper = glu.supno(jcol); 99 Index nextl = glu.xlsub(jcol); 100 VectorBlock<IndexVector> marker2(marker, 2*m, m); 101 102 103 column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu, *this); 104 105 // For each nonzero in A(*,jcol) do dfs 106 for (Index k = 0; ((k < m) ? lsub_col[k] != emptyIdxLU : false) ; k++) 107 { 108 Index krow = lsub_col(k); 109 lsub_col(k) = emptyIdxLU; 110 Index kmark = marker2(krow); 111 112 // krow was visited before, go to the next nonz; 113 if (kmark == jcol) continue; 114 115 dfs_kernel(StorageIndex(jcol), perm_r, nseg, glu.lsub, segrep, repfnz, xprune, marker2, parent, 116 xplore, glu, nextl, krow, traits); 117 } // for each nonzero ... 118 119 Index fsupc; 120 StorageIndex nsuper = glu.supno(jcol); 121 StorageIndex jcolp1 = StorageIndex(jcol) + 1; 122 Index jcolm1 = jcol - 1; 123 124 // check to see if j belongs in the same supernode as j-1 125 if ( jcol == 0 ) 126 { // Do nothing for column 0 127 nsuper = glu.supno(0) = 0 ; 128 } 129 else 130 { 131 fsupc = glu.xsup(nsuper); 132 StorageIndex jptr = glu.xlsub(jcol); // Not yet compressed 133 StorageIndex jm1ptr = glu.xlsub(jcolm1); 134 135 // Use supernodes of type T2 : see SuperLU paper 136 if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = emptyIdxLU; 137 138 // Make sure the number of columns in a supernode doesn't 139 // exceed threshold 140 if ( (jcol - fsupc) >= maxsuper) jsuper = emptyIdxLU; 141 142 /* If jcol starts a new supernode, reclaim storage space in 143 * glu.lsub from previous supernode. Note we only store 144 * the subscript set of the first and last columns of 145 * a supernode. (first for num values, last for pruning) 146 */ 147 if (jsuper == emptyIdxLU) 148 { // starts a new supernode 149 if ( (fsupc < jcolm1-1) ) 150 { // >= 3 columns in nsuper 151 StorageIndex ito = glu.xlsub(fsupc+1); 152 glu.xlsub(jcolm1) = ito; 153 StorageIndex istop = ito + jptr - jm1ptr; 154 xprune(jcolm1) = istop; // intialize xprune(jcol-1) 155 glu.xlsub(jcol) = istop; 156 157 for (StorageIndex ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito) 158 glu.lsub(ito) = glu.lsub(ifrom); 159 nextl = ito; // = istop + length(jcol) 160 } 161 nsuper++; 162 glu.supno(jcol) = nsuper; 163 } // if a new supernode 164 } // end else: jcol > 0 165 166 // Tidy up the pointers before exit 167 glu.xsup(nsuper+1) = jcolp1; 168 glu.supno(jcolp1) = nsuper; 169 xprune(jcol) = StorageIndex(nextl); // Intialize upper bound for pruning 170 glu.xlsub(jcolp1) = StorageIndex(nextl); 171 172 return 0; 173 } 174 175 } // end namespace internal 176 177 } // end namespace Eigen 178 179 #endif 180
