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 * NOTE: This file is the modified version of [s,d,c,z]copy_to_ucol.c file in SuperLU
12
13 * -- SuperLU routine (version 2.0) --
14 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
15 * and Lawrence Berkeley National Lab.
16 * November 15, 1997
17 *
18 * Copyright (c) 1994 by Xerox Corporation. All rights reserved.
19 *
20 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
21 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
22 *
23 * Permission is hereby granted to use or copy this program for any
24 * purpose, provided the above notices are retained on all copies.
25 * Permission to modify the code and to distribute modified code is
26 * granted, provided the above notices are retained, and a notice that
27 * the code was modified is included with the above copyright notice.
28 */
29 #ifndef SPARSELU_COPY_TO_UCOL_H
30 #define SPARSELU_COPY_TO_UCOL_H
31
32 namespace Eigen {
33 namespace internal {
34
35 /**
36 * \brief Performs numeric block updates (sup-col) in topological order
37 *
38 * \param jcol current column to update
39 * \param nseg Number of segments in the U part
40 * \param segrep segment representative ...
41 * \param repfnz First nonzero column in each row ...
42 * \param perm_r Row permutation
43 * \param dense Store the full representation of the column
44 * \param glu Global LU data.
45 * \return 0 - successful return
46 * > 0 - number of bytes allocated when run out of space
47 *
48 */
49 template <typename Scalar, typename Index>
copy_to_ucol(const Index jcol,const Index nseg,IndexVector & segrep,BlockIndexVector repfnz,IndexVector & perm_r,BlockScalarVector dense,GlobalLU_t & glu)50 Index SparseLUImpl<Scalar,Index>::copy_to_ucol(const Index jcol, const Index nseg, IndexVector& segrep, BlockIndexVector repfnz ,IndexVector& perm_r, BlockScalarVector dense, GlobalLU_t& glu)
51 {
52 Index ksub, krep, ksupno;
53
54 Index jsupno = glu.supno(jcol);
55
56 // For each nonzero supernode segment of U[*,j] in topological order
57 Index k = nseg - 1, i;
58 Index nextu = glu.xusub(jcol);
59 Index kfnz, isub, segsize;
60 Index new_next,irow;
61 Index fsupc, mem;
62 for (ksub = 0; ksub < nseg; ksub++)
63 {
64 krep = segrep(k); k--;
65 ksupno = glu.supno(krep);
66 if (jsupno != ksupno ) // should go into ucol();
67 {
68 kfnz = repfnz(krep);
69 if (kfnz != emptyIdxLU)
70 { // Nonzero U-segment
71 fsupc = glu.xsup(ksupno);
72 isub = glu.xlsub(fsupc) + kfnz - fsupc;
73 segsize = krep - kfnz + 1;
74 new_next = nextu + segsize;
75 while (new_next > glu.nzumax)
76 {
77 mem = memXpand<ScalarVector>(glu.ucol, glu.nzumax, nextu, UCOL, glu.num_expansions);
78 if (mem) return mem;
79 mem = memXpand<IndexVector>(glu.usub, glu.nzumax, nextu, USUB, glu.num_expansions);
80 if (mem) return mem;
81
82 }
83
84 for (i = 0; i < segsize; i++)
85 {
86 irow = glu.lsub(isub);
87 glu.usub(nextu) = perm_r(irow); // Unlike the L part, the U part is stored in its final order
88 glu.ucol(nextu) = dense(irow);
89 dense(irow) = Scalar(0.0);
90 nextu++;
91 isub++;
92 }
93
94 } // end nonzero U-segment
95
96 } // end if jsupno
97
98 } // end for each segment
99 glu.xusub(jcol + 1) = nextu; // close U(*,jcol)
100 return 0;
101 }
102
103 } // namespace internal
104 } // end namespace Eigen
105
106 #endif // SPARSELU_COPY_TO_UCOL_H
107