1*bf2c3715SXin Li // This file is part of Eigen, a lightweight C++ template library
2*bf2c3715SXin Li // for linear algebra.
3*bf2c3715SXin Li //
4*bf2c3715SXin Li // Copyright (C) 2012 Désiré Nuentsa-Wakam <[email protected]>
5*bf2c3715SXin Li //
6*bf2c3715SXin Li // This Source Code Form is subject to the terms of the Mozilla
7*bf2c3715SXin Li // Public License v. 2.0. If a copy of the MPL was not distributed
8*bf2c3715SXin Li // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9*bf2c3715SXin Li
10*bf2c3715SXin Li
11*bf2c3715SXin Li #ifndef EIGEN_SPARSELU_UTILS_H
12*bf2c3715SXin Li #define EIGEN_SPARSELU_UTILS_H
13*bf2c3715SXin Li
14*bf2c3715SXin Li namespace Eigen {
15*bf2c3715SXin Li namespace internal {
16*bf2c3715SXin Li
17*bf2c3715SXin Li /**
18*bf2c3715SXin Li * \brief Count Nonzero elements in the factors
19*bf2c3715SXin Li */
20*bf2c3715SXin Li template <typename Scalar, typename StorageIndex>
countnz(const Index n,Index & nnzL,Index & nnzU,GlobalLU_t & glu)21*bf2c3715SXin Li void SparseLUImpl<Scalar,StorageIndex>::countnz(const Index n, Index& nnzL, Index& nnzU, GlobalLU_t& glu)
22*bf2c3715SXin Li {
23*bf2c3715SXin Li nnzL = 0;
24*bf2c3715SXin Li nnzU = (glu.xusub)(n);
25*bf2c3715SXin Li Index nsuper = (glu.supno)(n);
26*bf2c3715SXin Li Index jlen;
27*bf2c3715SXin Li Index i, j, fsupc;
28*bf2c3715SXin Li if (n <= 0 ) return;
29*bf2c3715SXin Li // For each supernode
30*bf2c3715SXin Li for (i = 0; i <= nsuper; i++)
31*bf2c3715SXin Li {
32*bf2c3715SXin Li fsupc = glu.xsup(i);
33*bf2c3715SXin Li jlen = glu.xlsub(fsupc+1) - glu.xlsub(fsupc);
34*bf2c3715SXin Li
35*bf2c3715SXin Li for (j = fsupc; j < glu.xsup(i+1); j++)
36*bf2c3715SXin Li {
37*bf2c3715SXin Li nnzL += jlen;
38*bf2c3715SXin Li nnzU += j - fsupc + 1;
39*bf2c3715SXin Li jlen--;
40*bf2c3715SXin Li }
41*bf2c3715SXin Li }
42*bf2c3715SXin Li }
43*bf2c3715SXin Li
44*bf2c3715SXin Li /**
45*bf2c3715SXin Li * \brief Fix up the data storage lsub for L-subscripts.
46*bf2c3715SXin Li *
47*bf2c3715SXin Li * It removes the subscripts sets for structural pruning,
48*bf2c3715SXin Li * and applies permutation to the remaining subscripts
49*bf2c3715SXin Li *
50*bf2c3715SXin Li */
51*bf2c3715SXin Li template <typename Scalar, typename StorageIndex>
fixupL(const Index n,const IndexVector & perm_r,GlobalLU_t & glu)52*bf2c3715SXin Li void SparseLUImpl<Scalar,StorageIndex>::fixupL(const Index n, const IndexVector& perm_r, GlobalLU_t& glu)
53*bf2c3715SXin Li {
54*bf2c3715SXin Li Index fsupc, i, j, k, jstart;
55*bf2c3715SXin Li
56*bf2c3715SXin Li StorageIndex nextl = 0;
57*bf2c3715SXin Li Index nsuper = (glu.supno)(n);
58*bf2c3715SXin Li
59*bf2c3715SXin Li // For each supernode
60*bf2c3715SXin Li for (i = 0; i <= nsuper; i++)
61*bf2c3715SXin Li {
62*bf2c3715SXin Li fsupc = glu.xsup(i);
63*bf2c3715SXin Li jstart = glu.xlsub(fsupc);
64*bf2c3715SXin Li glu.xlsub(fsupc) = nextl;
65*bf2c3715SXin Li for (j = jstart; j < glu.xlsub(fsupc + 1); j++)
66*bf2c3715SXin Li {
67*bf2c3715SXin Li glu.lsub(nextl) = perm_r(glu.lsub(j)); // Now indexed into P*A
68*bf2c3715SXin Li nextl++;
69*bf2c3715SXin Li }
70*bf2c3715SXin Li for (k = fsupc+1; k < glu.xsup(i+1); k++)
71*bf2c3715SXin Li glu.xlsub(k) = nextl; // other columns in supernode i
72*bf2c3715SXin Li }
73*bf2c3715SXin Li
74*bf2c3715SXin Li glu.xlsub(n) = nextl;
75*bf2c3715SXin Li }
76*bf2c3715SXin Li
77*bf2c3715SXin Li } // end namespace internal
78*bf2c3715SXin Li
79*bf2c3715SXin Li } // end namespace Eigen
80*bf2c3715SXin Li #endif // EIGEN_SPARSELU_UTILS_H
81