src/SparseLU/SparseLU_copy_to_ucol.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 /*
0010
0011  * NOTE: This file is the modified version of [s,d,c,z]copy_to_ucol.c file in SuperLU
0012
0013  * -- SuperLU routine (version 2.0) --
0014  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
0015  * and Lawrence Berkeley National Lab.
0016  * November 15, 1997
0017  *
0018  * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
0019  *
0020  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
0021  * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
0022  *
0023  * Permission is hereby granted to use or copy this program for any
0024  * purpose, provided the above notices are retained on all copies.
0025  * Permission to modify the code and to distribute modified code is
0026  * granted, provided the above notices are retained, and a notice that
0027  * the code was modified is included with the above copyright notice.
0028  */
0029 #ifndef SPARSELU_COPY_TO_UCOL_H
0030 #define SPARSELU_COPY_TO_UCOL_H
0031
0032 namespace RivetEigen {
0033 namespace internal {
0034
0035 /**
0036  * \brief Performs numeric block updates (sup-col) in topological order
0037  *
0038  * \param jcol current column to update
0039  * \param nseg Number of segments in the U part
0040  * \param segrep segment representative ...
0041  * \param repfnz First nonzero column in each row  ...
0042  * \param perm_r Row permutation
0043  * \param dense Store the full representation of the column
0044  * \param glu Global LU data.
0045  * \return 0 - successful return
0046  *         > 0 - number of bytes allocated when run out of space
0047  *
0048  */
0049 template <typename Scalar, typename StorageIndex>
0050 Index SparseLUImpl<Scalar,StorageIndex>::copy_to_ucol(const Index jcol, const Index nseg, IndexVector& segrep,
0051                                                       BlockIndexVector repfnz ,IndexVector& perm_r, BlockScalarVector dense, GlobalLU_t& glu)
0052 {
0053   Index ksub, krep, ksupno;
0054
0055   Index jsupno = glu.supno(jcol);
0056
0057   // For each nonzero supernode segment of U[*,j] in topological order
0058   Index k = nseg - 1, i;
0059   StorageIndex nextu = glu.xusub(jcol);
0060   Index kfnz, isub, segsize;
0061   Index new_next,irow;
0062   Index fsupc, mem;
0063   for (ksub = 0; ksub < nseg; ksub++)
0064   {
0065     krep = segrep(k); k--;
0066     ksupno = glu.supno(krep);
0067     if (jsupno != ksupno ) // should go into ucol();
0068     {
0069       kfnz = repfnz(krep);
0070       if (kfnz != emptyIdxLU)
0071       { // Nonzero U-segment
0072         fsupc = glu.xsup(ksupno);
0073         isub = glu.xlsub(fsupc) + kfnz - fsupc;
0074         segsize = krep - kfnz + 1;
0075         new_next = nextu + segsize;
0076         while (new_next > glu.nzumax)
0077         {
0078           mem = memXpand<ScalarVector>(glu.ucol, glu.nzumax, nextu, UCOL, glu.num_expansions);
0079           if (mem) return mem;
0080           mem = memXpand<IndexVector>(glu.usub, glu.nzumax, nextu, USUB, glu.num_expansions);
0081           if (mem) return mem;
0082
0083         }
0084
0085         for (i = 0; i < segsize; i++)
0086         {
0087           irow = glu.lsub(isub);
0088           glu.usub(nextu) = perm_r(irow); // Unlike the L part, the U part is stored in its final order
0089           glu.ucol(nextu) = dense(irow);
0090           dense(irow) = Scalar(0.0);
0091           nextu++;
0092           isub++;
0093         }
0094
0095       } // end nonzero U-segment
0096
0097     } // end if jsupno
0098
0099   } // end for each segment
0100   glu.xusub(jcol + 1) = nextu; // close U(*,jcol)
0101   return 0;
0102 }
0103
0104 } // namespace internal
0105 } // end namespace RivetEigen
0106
0107 #endif // SPARSELU_COPY_TO_UCOL_H