cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

SparseLU_Utils.h (2049B)

      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@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