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 #ifndef EIGEN_SPARSELU_UTILS_H 12 #define EIGEN_SPARSELU_UTILS_H 13 14 namespace Eigen { 15 namespace internal { 16 17 /** 18 * \brief Count Nonzero elements in the factors 19 */ 20 template <typename Scalar, typename StorageIndex> 21 void SparseLUImpl<Scalar,StorageIndex>::countnz(const Index n, Index& nnzL, Index& nnzU, GlobalLU_t& glu) 22 { 23 nnzL = 0; 24 nnzU = (glu.xusub)(n); 25 Index nsuper = (glu.supno)(n); 26 Index jlen; 27 Index i, j, fsupc; 28 if (n <= 0 ) return; 29 // For each supernode 30 for (i = 0; i <= nsuper; i++) 31 { 32 fsupc = glu.xsup(i); 33 jlen = glu.xlsub(fsupc+1) - glu.xlsub(fsupc); 34 35 for (j = fsupc; j < glu.xsup(i+1); j++) 36 { 37 nnzL += jlen; 38 nnzU += j - fsupc + 1; 39 jlen--; 40 } 41 } 42 } 43 44 /** 45 * \brief Fix up the data storage lsub for L-subscripts. 46 * 47 * It removes the subscripts sets for structural pruning, 48 * and applies permutation to the remaining subscripts 49 * 50 */ 51 template <typename Scalar, typename StorageIndex> 52 void SparseLUImpl<Scalar,StorageIndex>::fixupL(const Index n, const IndexVector& perm_r, GlobalLU_t& glu) 53 { 54 Index fsupc, i, j, k, jstart; 55 56 StorageIndex nextl = 0; 57 Index nsuper = (glu.supno)(n); 58 59 // For each supernode 60 for (i = 0; i <= nsuper; i++) 61 { 62 fsupc = glu.xsup(i); 63 jstart = glu.xlsub(fsupc); 64 glu.xlsub(fsupc) = nextl; 65 for (j = jstart; j < glu.xlsub(fsupc + 1); j++) 66 { 67 glu.lsub(nextl) = perm_r(glu.lsub(j)); // Now indexed into P*A 68 nextl++; 69 } 70 for (k = fsupc+1; k < glu.xsup(i+1); k++) 71 glu.xlsub(k) = nextl; // other columns in supernode i 72 } 73 74 glu.xlsub(n) = nextl; 75 } 76 77 } // end namespace internal 78 79 } // end namespace Eigen 80 #endif // EIGEN_SPARSELU_UTILS_H 81
