Skip to content

Commit cac9358

Browse files
Wentzellclaude
andcommitted
Fix gelss rcond parameter for MKL compatibility
Use rcond=-1.0 (machine precision) instead of 0.0 in gelss calls. MKL interprets rcond=0.0 as "use all singular values" which causes numerical instability for rank-deficient systems, while NETLIB treats it as machine precision. Using -1.0 explicitly requests machine precision and works correctly with both implementations. Co-Authored-By: Claude <noreply@anthropic.com>
1 parent 897acf5 commit cac9358

1 file changed

Lines changed: 4 additions & 4 deletions

File tree

c++/cppdlr2d/dlr2d.cpp

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -464,7 +464,7 @@ vals2coefs(int r, fmatrix cf2if, nda::vector_const_view<dcomplex> vals,
464464

465465
auto s = nda::vector<double>(m); // Singular values (not needed)
466466
int rank = 0; // Rank (not needed)
467-
nda::lapack::gelss(cf2if, tmp, s, 0.0, rank);
467+
nda::lapack::gelss(cf2if, tmp, s, -1.0, rank);
468468

469469
auto coefreg = nda::zeros<dcomplex>(3, r, r);
470470
auto coefsng = nda::zeros<dcomplex>(r);
@@ -494,7 +494,7 @@ vals2coefs_many(int r, fmatrix cf2if,
494494

495495
auto s = nda::vector<double>(m); // Singular values (not needed)
496496
int rank = 0; // Rank (not needed)
497-
nda::lapack::gelss(cf2if, tmp, s, 0.0, rank);
497+
nda::lapack::gelss(cf2if, tmp, s, -1.0, rank);
498498

499499
auto coefreg = nda::zeros<dcomplex>(nrhs, 3, r, r);
500500
auto coefsng = nda::zeros<dcomplex>(nrhs, r);
@@ -523,7 +523,7 @@ vals2coefs_if_3term(fmatrix cf2if, nda::vector_const_view<dcomplex> vals,
523523

524524
auto s = nda::vector<double>(m); // Singular values (not needed)
525525
int rank = 0; // Rank (not needed)
526-
nda::lapack::gelss(cf2if, tmp, s, 0.0, rank);
526+
nda::lapack::gelss(cf2if, tmp, s, -1.0, rank);
527527

528528
auto coefreg = nda::array<dcomplex, 3>(2, r, r);
529529
auto coefsng = nda::array<dcomplex, 1>(r);
@@ -546,7 +546,7 @@ vals2coefs_if_many_3term(fmatrix cf2if,
546546

547547
auto s = nda::vector<double>(m); // Singular values (not needed)
548548
int rank = 0; // Rank (not needed)
549-
nda::lapack::gelss(cf2if, tmp, s, 0.0, rank);
549+
nda::lapack::gelss(cf2if, tmp, s, -1.0, rank);
550550

551551
auto coefreg = nda::array<dcomplex, 4>(nrhs, 2, r, r);
552552
auto coefsng = nda::array<dcomplex, 2>(nrhs, r);

0 commit comments

Comments
 (0)