17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/* 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * NOTE: This file is the modified version of [s,d,c,z]pruneL.c file in SuperLU 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * -- SuperLU routine (version 2.0) -- 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Univ. of California Berkeley, Xerox Palo Alto Research Center, 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and Lawrence Berkeley National Lab. 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * November 15, 1997 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Copyright (c) 1994 by Xerox Corporation. All rights reserved. 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Permission is hereby granted to use or copy this program for any 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * purpose, provided the above notices are retained on all copies. 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Permission to modify the code and to distribute modified code is 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * granted, provided the above notices are retained, and a notice that 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * the code was modified is included with the above copyright notice. 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef SPARSELU_PRUNEL_H 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define SPARSELU_PRUNEL_H 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen { 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal { 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief Prunes the L-structure. 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * It prunes the L-structure of supernodes whose L-structure contains the current pivot row "pivrow" 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param jcol The current column of L 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] perm_r Row permutation 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[out] pivrow The pivot row 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param nseg Number of segments 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param segrep 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param repfnz 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[out] xprune 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param glu Global LU data 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <typename Scalar, typename StorageIndex> 532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangvoid SparseLUImpl<Scalar,StorageIndex>::pruneL(const Index jcol, const IndexVector& perm_r, const Index pivrow, const Index nseg, 542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, GlobalLU_t& glu) 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // For each supernode-rep irep in U(*,j] 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index jsupno = glu.supno(jcol); 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index i,irep,irep1; 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool movnum, do_prune = false; 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index kmin = 0, kmax = 0, minloc, maxloc,krow; 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (i = 0; i < nseg; i++) 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez irep = segrep(i); 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez irep1 = irep + 1; 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez do_prune = false; 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Don't prune with a zero U-segment 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (repfnz(irep) == emptyIdxLU) continue; 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // If a snode overlaps with the next panel, then the U-segment 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // is fragmented into two parts -- irep and irep1. We should let 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // pruning occur at the rep-column in irep1s snode. 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (glu.supno(irep) == glu.supno(irep1) ) continue; // don't prune 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // If it has not been pruned & it has a nonz in row L(pivrow,i) 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (glu.supno(irep) != jsupno ) 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if ( xprune (irep) >= glu.xlsub(irep1) ) 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez kmin = glu.xlsub(irep); 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez kmax = glu.xlsub(irep1) - 1; 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (krow = kmin; krow <= kmax; krow++) 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (glu.lsub(krow) == pivrow) 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez do_prune = true; 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez break; 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (do_prune) 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do a quicksort-type partition 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // movnum=true means that the num values have to be exchanged 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez movnum = false; 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (irep == glu.xsup(glu.supno(irep)) ) // Snode of size 1 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez movnum = true; 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez while (kmin <= kmax) 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (perm_r(glu.lsub(kmax)) == emptyIdxLU) 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez kmax--; 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if ( perm_r(glu.lsub(kmin)) != emptyIdxLU) 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez kmin++; 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // kmin below pivrow (not yet pivoted), and kmax 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // above pivrow: interchange the two suscripts 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::swap(glu.lsub(kmin), glu.lsub(kmax)); 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // If the supernode has only one column, then we 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // only keep one set of subscripts. For any subscript 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // intercnahge performed, similar interchange must be 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // done on the numerical values. 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (movnum) 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez minloc = glu.xlusup(irep) + ( kmin - glu.xlsub(irep) ); 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez maxloc = glu.xlusup(irep) + ( kmax - glu.xlsub(irep) ); 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::swap(glu.lusup(minloc), glu.lusup(maxloc)); 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez kmin++; 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez kmax--; 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // end while 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang xprune(irep) = StorageIndex(kmin); //Pruning 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // end if do_prune 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // end pruning 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // End for each U-segment 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace Eigen 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif // SPARSELU_PRUNEL_H 137