@@ -1872,6 +1872,11 @@ i_t basis_update_mpf_t<i_t, f_t>::update(const std::vector<f_t>& utilde,
18721872 S_.x .data () + S_.col_start [S_start + 1 ],
18731873 v_nz);
18741874
1875+ if (std::abs (mu) < 1e-13 ) {
1876+ // Force a refactor. Otherwise we will get numerical issues when dividing by mu.
1877+ return 1 ;
1878+ }
1879+
18751880#ifdef CHECK_MU
18761881 const f_t mu_check = 1.0 + dot<i_t , f_t >(etilde, u);
18771882 printf (" Update: mu %e mu_check %e diff %e\n " , mu, mu_check, std::abs (mu - mu_check));
@@ -1924,12 +1929,17 @@ i_t basis_update_mpf_t<i_t, f_t>::update(const sparse_vector_t<i_t, f_t>& utilde
19241929 S_.append_column (etilde);
19251930
19261931 // Compute mu = 1 + v^T * u
1927- mu_values_.push_back (1.0 + sparse_dot (S_.i .data () + S_.col_start [S_start],
1928- S_.x .data () + S_.col_start [S_start],
1929- S_.col_start [S_start + 1 ] - S_.col_start [S_start],
1930- S_.i .data () + S_.col_start [S_start + 1 ],
1931- S_.x .data () + S_.col_start [S_start + 1 ],
1932- S_.col_start [S_start + 2 ] - S_.col_start [S_start + 1 ]));
1932+ const f_t mu = 1.0 + sparse_dot (S_.i .data () + S_.col_start [S_start],
1933+ S_.x .data () + S_.col_start [S_start],
1934+ S_.col_start [S_start + 1 ] - S_.col_start [S_start],
1935+ S_.i .data () + S_.col_start [S_start + 1 ],
1936+ S_.x .data () + S_.col_start [S_start + 1 ],
1937+ S_.col_start [S_start + 2 ] - S_.col_start [S_start + 1 ]);
1938+ if (std::abs (mu) < 1e-13 ) {
1939+ // Force a refactor. Otherwise we will get numerical issues when dividing by mu.
1940+ return 1 ;
1941+ }
1942+ mu_values_.push_back (mu);
19331943 // Clear the workspace
19341944 for (i_t k = 0 ; k < nz; ++k) {
19351945 const i_t i = xi_workspace_[m + k];
0 commit comments