src/SparseLU/SparseLU_pivotL.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
0012  * NOTE: This file is the modified version of xpivotL.c file in SuperLU
0013
0014  * -- SuperLU routine (version 3.0) --
0015  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
0016  * and Lawrence Berkeley National Lab.
0017  * October 15, 2003
0018  *
0019  * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
0020  *
0021  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
0022  * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
0023  *
0024  * Permission is hereby granted to use or copy this program for any
0025  * purpose, provided the above notices are retained on all copies.
0026  * Permission to modify the code and to distribute modified code is
0027  * granted, provided the above notices are retained, and a notice that
0028  * the code was modified is included with the above copyright notice.
0029  */
0030 #ifndef SPARSELU_PIVOTL_H
0031 #define SPARSELU_PIVOTL_H
0032
0033 namespace RivetEigen {
0034 namespace internal {
0035
0036 /**
0037  * \brief Performs the numerical pivotin on the current column of L, and the CDIV operation.
0038  *
0039  * Pivot policy :
0040  * (1) Compute thresh = u * max_(i>=j) abs(A_ij);
0041  * (2) IF user specifies pivot row k and abs(A_kj) >= thresh THEN
0042  *           pivot row = k;
0043  *       ELSE IF abs(A_jj) >= thresh THEN
0044  *           pivot row = j;
0045  *       ELSE
0046  *           pivot row = m;
0047  *
0048  *   Note: If you absolutely want to use a given pivot order, then set u=0.0.
0049  *
0050  * \param jcol The current column of L
0051  * \param diagpivotthresh diagonal pivoting threshold
0052  * \param[in,out] perm_r Row permutation (threshold pivoting)
0053  * \param[in] iperm_c column permutation - used to finf diagonal of Pc*A*Pc'
0054  * \param[out] pivrow  The pivot row
0055  * \param glu Global LU data
0056  * \return 0 if success, i > 0 if U(i,i) is exactly zero
0057  *
0058  */
0059 template <typename Scalar, typename StorageIndex>
0060 Index SparseLUImpl<Scalar,StorageIndex>::pivotL(const Index jcol, const RealScalar& diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, Index& pivrow, GlobalLU_t& glu)
0061 {
0062
0063   Index fsupc = (glu.xsup)((glu.supno)(jcol)); // First column in the supernode containing the column jcol
0064   Index nsupc = jcol - fsupc; // Number of columns in the supernode portion, excluding jcol; nsupc >=0
0065   Index lptr = glu.xlsub(fsupc); // pointer to the starting location of the row subscripts for this supernode portion
0066   Index nsupr = glu.xlsub(fsupc+1) - lptr; // Number of rows in the supernode
0067   Index lda = glu.xlusup(fsupc+1) - glu.xlusup(fsupc); // leading dimension
0068   Scalar* lu_sup_ptr = &(glu.lusup.data()[glu.xlusup(fsupc)]); // Start of the current supernode
0069   Scalar* lu_col_ptr = &(glu.lusup.data()[glu.xlusup(jcol)]); // Start of jcol in the supernode
0070   StorageIndex* lsub_ptr = &(glu.lsub.data()[lptr]); // Start of row indices of the supernode
0071
0072   // Determine the largest abs numerical value for partial pivoting
0073   Index diagind = iperm_c(jcol); // diagonal index
0074   RealScalar pivmax(-1.0);
0075   Index pivptr = nsupc;
0076   Index diag = emptyIdxLU;
0077   RealScalar rtemp;
0078   Index isub, icol, itemp, k;
0079   for (isub = nsupc; isub < nsupr; ++isub) {
0080     using std::abs;
0081     rtemp = abs(lu_col_ptr[isub]);
0082     if (rtemp > pivmax) {
0083       pivmax = rtemp;
0084       pivptr = isub;
0085     }
0086     if (lsub_ptr[isub] == diagind) diag = isub;
0087   }
0088
0089   // Test for singularity
0090   if ( pivmax <= RealScalar(0.0) ) {
0091     // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero
0092     pivrow = pivmax < RealScalar(0.0) ? diagind : lsub_ptr[pivptr];
0093     perm_r(pivrow) = StorageIndex(jcol);
0094     return (jcol+1);
0095   }
0096
0097   RealScalar thresh = diagpivotthresh * pivmax;
0098
0099   // Choose appropriate pivotal element
0100
0101   {
0102     // Test if the diagonal element can be used as a pivot (given the threshold value)
0103     if (diag >= 0 )
0104     {
0105       // Diagonal element exists
0106       using std::abs;
0107       rtemp = abs(lu_col_ptr[diag]);
0108       if (rtemp != RealScalar(0.0) && rtemp >= thresh) pivptr = diag;
0109     }
0110     pivrow = lsub_ptr[pivptr];
0111   }
0112
0113   // Record pivot row
0114   perm_r(pivrow) = StorageIndex(jcol);
0115   // Interchange row subscripts
0116   if (pivptr != nsupc )
0117   {
0118     std::swap( lsub_ptr[pivptr], lsub_ptr[nsupc] );
0119     // Interchange numerical values as well, for the two rows in the whole snode
0120     // such that L is indexed the same way as A
0121     for (icol = 0; icol <= nsupc; icol++)
0122     {
0123       itemp = pivptr + icol * lda;
0124       std::swap(lu_sup_ptr[itemp], lu_sup_ptr[nsupc + icol * lda]);
0125     }
0126   }
0127   // cdiv operations
0128   Scalar temp = Scalar(1.0) / lu_col_ptr[nsupc];
0129   for (k = nsupc+1; k < nsupr; k++)
0130     lu_col_ptr[k] *= temp;
0131   return 0;
0132 }
0133
0134 } // end namespace internal
0135 } // end namespace RivetEigen
0136
0137 #endif // SPARSELU_PIVOTL_H