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