GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-14 04:35:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

correct bug in the complex version

Browse Source
This commit is contained in:
Desire NUENTSA 2012-07-19 18:15:23 +02:00
parent 59642da88b
commit 925ace196c
2 changed files with 8 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Eigen/src/SparseLU/SparseLU_pivotL.h
View File

@ -72,6 +72,7 @@ int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivott
{
typedef typename IndexVector::Scalar Index;
typedef typename ScalarVector::Scalar Scalar;
typedef typename ScalarVector::RealScalar RealScalar;
// Initialize pointers
IndexVector& lsub = glu.lsub; // Compressed row subscripts of L rectangular supernodes.
IndexVector& xlsub = glu.xlsub; // pointers to the beginning of each column subscript in lsub
@ -88,10 +89,10 @@ int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivott
// Determine the largest abs numerical value for partial pivoting
Index diagind = iperm_c(jcol); // diagonal index
Scalar pivmax = 0.0;
RealScalar pivmax = 0.0;
Index pivptr = nsupc;
Index diag = IND_EMPTY;
Scalar rtemp;
RealScalar rtemp;
Index isub, icol, itemp, k;
for (isub = nsupc; isub < nsupr; ++isub) {
rtemp = std::abs(lu_col_ptr[isub]);
@ -109,7 +110,7 @@ int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivott
return (jcol+1);
}
Scalar thresh = diagpivotthresh * pivmax;
RealScalar thresh = diagpivotthresh * pivmax;
// Choose appropriate pivotal element
@ -119,7 +120,7 @@ int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivott
{
// Diagonal element exists
rtemp = std::abs(lu_col_ptr[diag]);
if (rtemp != Scalar(0.0) && rtemp >= thresh) pivptr = diag;
if (rtemp != 0.0 && rtemp >= thresh) pivptr = diag;
}
pivrow = lsub_ptr[pivptr];
}

3
bench/spbench/test_sparseLU.cpp
View File

@ -14,6 +14,7 @@ using namespace Eigen;
int main(int argc, char **args)
{
typedef complex<double> scalar;
// typedef double scalar;
SparseMatrix<scalar, ColMajor> A;
typedef SparseMatrix<scalar, ColMajor>::Index Index;
typedef Matrix<scalar, Dynamic, Dynamic> DenseMatrix;
@ -34,7 +35,7 @@ int main(int argc, char **args)
bool iscomplex=false, isvector=false;
int sym;
getMarketHeader(args[1], sym, iscomplex, isvector);
if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; }
// if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; }
if (isvector) { cout << "The provided file is not a matrix file\n"; return -1;}
if (sym != 0) { // symmetric matrices, only the lower part is stored
SparseMatrix<scalar, ColMajor> temp;