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 * 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 StorageIndex> 50 Index SparseLUImpl<Scalar,StorageIndex>::copy_to_ucol(const Index jcol, const Index nseg, IndexVector& segrep, 51 BlockIndexVector repfnz ,IndexVector& perm_r, BlockScalarVector dense, GlobalLU_t& glu) 52 { 53 Index ksub, krep, ksupno; 54 55 Index jsupno = glu.supno(jcol); 56 57 // For each nonzero supernode segment of U[*,j] in topological order 58 Index k = nseg - 1, i; 59 StorageIndex nextu = glu.xusub(jcol); 60 Index kfnz, isub, segsize; 61 Index new_next,irow; 62 Index fsupc, mem; 63 for (ksub = 0; ksub < nseg; ksub++) 64 { 65 krep = segrep(k); k--; 66 ksupno = glu.supno(krep); 67 if (jsupno != ksupno ) // should go into ucol(); 68 { 69 kfnz = repfnz(krep); 70 if (kfnz != emptyIdxLU) 71 { // Nonzero U-segment 72 fsupc = glu.xsup(ksupno); 73 isub = glu.xlsub(fsupc) + kfnz - fsupc; 74 segsize = krep - kfnz + 1; 75 new_next = nextu + segsize; 76 while (new_next > glu.nzumax) 77 { 78 mem = memXpand<ScalarVector>(glu.ucol, glu.nzumax, nextu, UCOL, glu.num_expansions); 79 if (mem) return mem; 80 mem = memXpand<IndexVector>(glu.usub, glu.nzumax, nextu, USUB, glu.num_expansions); 81 if (mem) return mem; 82 83 } 84 85 for (i = 0; i < segsize; i++) 86 { 87 irow = glu.lsub(isub); 88 glu.usub(nextu) = perm_r(irow); // Unlike the L part, the U part is stored in its final order 89 glu.ucol(nextu) = dense(irow); 90 dense(irow) = Scalar(0.0); 91 nextu++; 92 isub++; 93 } 94 95 } // end nonzero U-segment 96 97 } // end if jsupno 98 99 } // end for each segment 100 glu.xusub(jcol + 1) = nextu; // close U(*,jcol) 101 return 0; 102 } 103 104 } // namespace internal 105 } // end namespace Eigen 106 107 #endif // SPARSELU_COPY_TO_UCOL_H 108
