devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix SparseLU and SparseQR for custom scalar types

Browse Source 
libeigen/eigen!2345
This commit is contained in:
Tyler Veness
2026-03-27 00:13:11 -07:00
committed by Rasmus Munk Larsen
parent 002229ce47
commit 9939a4c6e3
2 changed files with 9 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
Eigen/src/SparseLU/SparseLU_panel_bmod.h
View File

@@ -118,7 +118,7 @@ void SparseLUImpl<Scalar, StorageIndex>::panel_bmod(const Index m, const Index w
Index isub = lptr + no_zeros; Index isub = lptr + no_zeros;
Index off = u_rows - segsize; Index off = u_rows - segsize;
for (Index i = 0; i < off; i++) U(i, u_col) = 0; for (Index i = 0; i < off; i++) U(i, u_col) = Scalar(0);
for (Index i = 0; i < segsize; i++) { for (Index i = 0; i < segsize; i++) {
Index irow = glu.lsub(isub); Index irow = glu.lsub(isub);
U(i + off, u_col) = dense_col(irow); U(i + off, u_col) = dense_col(irow);
@@ -163,14 +163,14 @@ void SparseLUImpl<Scalar, StorageIndex>::panel_bmod(const Index m, const Index w
for (Index i = 0; i < segsize; i++) { for (Index i = 0; i < segsize; i++) {
Index irow = glu.lsub(isub++); Index irow = glu.lsub(isub++);
dense_col(irow) = U.coeff(i + off, u_col); dense_col(irow) = U.coeff(i + off, u_col);
U.coeffRef(i + off, u_col) = 0; U.coeffRef(i + off, u_col) = Scalar(0);
} }
// Scatter l into SPA dense[] // Scatter l into SPA dense[]
for (Index i = 0; i < nrow; i++) { for (Index i = 0; i < nrow; i++) {
Index irow = glu.lsub(isub++); Index irow = glu.lsub(isub++);
dense_col(irow) -= L.coeff(i, u_col); dense_col(irow) -= L.coeff(i, u_col);
L.coeffRef(i, u_col) = 0; L.coeffRef(i, u_col) = Scalar(0);
} }
u_col++; u_col++;
} }

12
Eigen/src/SparseQR/SparseQR.h
View File

@@ -401,10 +401,10 @@ void SparseQR<MatrixType, OrderingType>::factorize(const MatrixType& mat) {
*/ */
RealScalar pivotThreshold; RealScalar pivotThreshold;
if (m_useDefaultThreshold) { if (m_useDefaultThreshold) {
RealScalar max2Norm = 0.0; RealScalar max2Norm = RealScalar(0.0);
for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm()); for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
if (max2Norm == RealScalar(0)) max2Norm = RealScalar(1); if (max2Norm == RealScalar(0)) max2Norm = RealScalar(1);
pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon(); pivotThreshold = RealScalar(20 * (m + n)) * max2Norm * NumTraits<RealScalar>::epsilon();
} else { } else {
pivotThreshold = m_threshold; pivotThreshold = m_threshold;
} }
@@ -497,7 +497,7 @@ void SparseQR<MatrixType, OrderingType>::factorize(const MatrixType& mat) {
} // End update current column } // End update current column
Scalar tau = RealScalar(0); Scalar tau = RealScalar(0);
RealScalar beta = 0; RealScalar beta = RealScalar(0);
if (nonzeroCol < diagSize) { if (nonzeroCol < diagSize) {
// Compute the Householder reflection that eliminate the current column // Compute the Householder reflection that eliminate the current column
@@ -505,16 +505,16 @@ void SparseQR<MatrixType, OrderingType>::factorize(const MatrixType& mat) {
Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0); Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
// First, the squared norm of Q((col+1):m, col) // First, the squared norm of Q((col+1):m, col)
RealScalar sqrNorm = 0.; RealScalar sqrNorm = RealScalar(0.);
for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
if (sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) { if (sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) {
beta = numext::real(c0); beta = numext::real(c0);
tval(Qidx(0)) = 1; tval(Qidx(0)) = Scalar(1);
} else { } else {
using std::sqrt; using std::sqrt;
beta = sqrt(numext::abs2(c0) + sqrNorm); beta = sqrt(numext::abs2(c0) + sqrNorm);
if (numext::real(c0) >= RealScalar(0)) beta = -beta; if (numext::real(c0) >= RealScalar(0)) beta = -beta;
tval(Qidx(0)) = 1; tval(Qidx(0)) = Scalar(1);
for (Index itq = 1; itq < nzcolQ; ++itq) tval(Qidx(itq)) /= (c0 - beta); for (Index itq = 1; itq < nzcolQ; ++itq) tval(Qidx(itq)) /= (c0 - beta);
tau = numext::conj((beta - c0) / beta); tau = numext::conj((beta - c0) / beta);
} }