SparseLU_pivotL.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 /* 
00011  
00012  * NOTE: This file is the modified version of xpivotL.c file in SuperLU 
00013  
00014  * -- SuperLU routine (version 3.0) --
00015  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
00016  * and Lawrence Berkeley National Lab.
00017  * October 15, 2003
00018  *
00019  * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
00020  *
00021  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
00022  * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
00023  *
00024  * Permission is hereby granted to use or copy this program for any
00025  * purpose, provided the above notices are retained on all copies.
00026  * Permission to modify the code and to distribute modified code is
00027  * granted, provided the above notices are retained, and a notice that
00028  * the code was modified is included with the above copyright notice.
00029  */
00030 #ifndef SPARSELU_PIVOTL_H
00031 #define SPARSELU_PIVOTL_H
00032 
00033 namespace Eigen {
00034 namespace internal {
00035   
00059 template <typename Scalar, typename Index>
00060 Index SparseLUImpl<Scalar,Index>::pivotL(const Index jcol, const RealScalar& diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, Index& pivrow, GlobalLU_t& glu)
00061 {
00062   
00063   Index fsupc = (glu.xsup)((glu.supno)(jcol)); // First column in the supernode containing the column jcol
00064   Index nsupc = jcol - fsupc; // Number of columns in the supernode portion, excluding jcol; nsupc >=0
00065   Index lptr = glu.xlsub(fsupc); // pointer to the starting location of the row subscripts for this supernode portion
00066   Index nsupr = glu.xlsub(fsupc+1) - lptr; // Number of rows in the supernode
00067   Index lda = glu.xlusup(fsupc+1) - glu.xlusup(fsupc); // leading dimension
00068   Scalar* lu_sup_ptr = &(glu.lusup.data()[glu.xlusup(fsupc)]); // Start of the current supernode
00069   Scalar* lu_col_ptr = &(glu.lusup.data()[glu.xlusup(jcol)]); // Start of jcol in the supernode
00070   Index* lsub_ptr = &(glu.lsub.data()[lptr]); // Start of row indices of the supernode
00071   
00072   // Determine the largest abs numerical value for partial pivoting 
00073   Index diagind = iperm_c(jcol); // diagonal index 
00074   RealScalar pivmax = 0.0; 
00075   Index pivptr = nsupc; 
00076   Index diag = emptyIdxLU; 
00077   RealScalar rtemp;
00078   Index isub, icol, itemp, k; 
00079   for (isub = nsupc; isub < nsupr; ++isub) {
00080     rtemp = std::abs(lu_col_ptr[isub]);
00081     if (rtemp > pivmax) {
00082       pivmax = rtemp; 
00083       pivptr = isub;
00084     } 
00085     if (lsub_ptr[isub] == diagind) diag = isub;
00086   }
00087   
00088   // Test for singularity
00089   if ( pivmax == 0.0 ) {
00090     pivrow = lsub_ptr[pivptr];
00091     perm_r(pivrow) = jcol;
00092     return (jcol+1);
00093   }
00094   
00095   RealScalar thresh = diagpivotthresh * pivmax; 
00096   
00097   // Choose appropriate pivotal element 
00098   
00099   {
00100     // Test if the diagonal element can be used as a pivot (given the threshold value)
00101     if (diag >= 0 ) 
00102     {
00103       // Diagonal element exists
00104       rtemp = std::abs(lu_col_ptr[diag]);
00105       if (rtemp != 0.0 && rtemp >= thresh) pivptr = diag;
00106     }
00107     pivrow = lsub_ptr[pivptr];
00108   }
00109   
00110   // Record pivot row
00111   perm_r(pivrow) = jcol; 
00112   // Interchange row subscripts
00113   if (pivptr != nsupc )
00114   {
00115     std::swap( lsub_ptr[pivptr], lsub_ptr[nsupc] );
00116     // Interchange numerical values as well, for the two rows in the whole snode
00117     // such that L is indexed the same way as A
00118     for (icol = 0; icol <= nsupc; icol++)
00119     {
00120       itemp = pivptr + icol * lda; 
00121       std::swap(lu_sup_ptr[itemp], lu_sup_ptr[nsupc + icol * lda]);
00122     }
00123   }
00124   // cdiv operations
00125   Scalar temp = Scalar(1.0) / lu_col_ptr[nsupc];
00126   for (k = nsupc+1; k < nsupr; k++)
00127     lu_col_ptr[k] *= temp; 
00128   return 0;
00129 }
00130 
00131 } // end namespace internal
00132 } // end namespace Eigen
00133 
00134 #endif // SPARSELU_PIVOTL_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 12:00:34