[136] | 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>
|
---|
| 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
|
---|