![]() |
Eigen
3.3.3
|
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 StorageIndex> 00060 Index SparseLUImpl<Scalar,StorageIndex>::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 StorageIndex* 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(-1.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 using std::abs; 00081 rtemp = abs(lu_col_ptr[isub]); 00082 if (rtemp > pivmax) { 00083 pivmax = rtemp; 00084 pivptr = isub; 00085 } 00086 if (lsub_ptr[isub] == diagind) diag = isub; 00087 } 00088 00089 // Test for singularity 00090 if ( pivmax <= RealScalar(0.0) ) { 00091 // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero 00092 pivrow = pivmax < RealScalar(0.0) ? diagind : lsub_ptr[pivptr]; 00093 perm_r(pivrow) = StorageIndex(jcol); 00094 return (jcol+1); 00095 } 00096 00097 RealScalar thresh = diagpivotthresh * pivmax; 00098 00099 // Choose appropriate pivotal element 00100 00101 { 00102 // Test if the diagonal element can be used as a pivot (given the threshold value) 00103 if (diag >= 0 ) 00104 { 00105 // Diagonal element exists 00106 using std::abs; 00107 rtemp = abs(lu_col_ptr[diag]); 00108 if (rtemp != RealScalar(0.0) && rtemp >= thresh) pivptr = diag; 00109 } 00110 pivrow = lsub_ptr[pivptr]; 00111 } 00112 00113 // Record pivot row 00114 perm_r(pivrow) = StorageIndex(jcol); 00115 // Interchange row subscripts 00116 if (pivptr != nsupc ) 00117 { 00118 std::swap( lsub_ptr[pivptr], lsub_ptr[nsupc] ); 00119 // Interchange numerical values as well, for the two rows in the whole snode 00120 // such that L is indexed the same way as A 00121 for (icol = 0; icol <= nsupc; icol++) 00122 { 00123 itemp = pivptr + icol * lda; 00124 std::swap(lu_sup_ptr[itemp], lu_sup_ptr[nsupc + icol * lda]); 00125 } 00126 } 00127 // cdiv operations 00128 Scalar temp = Scalar(1.0) / lu_col_ptr[nsupc]; 00129 for (k = nsupc+1; k < nsupr; k++) 00130 lu_col_ptr[k] *= temp; 00131 return 0; 00132 } 00133 00134 } // end namespace internal 00135 } // end namespace Eigen 00136 00137 #endif // SPARSELU_PIVOTL_H