@@ -103,35 +103,30 @@ VectorXd C3Plus::SolveSingleProjection(const MatrixXd& U,
103103 const int & warm_start_index) {
104104 VectorXd delta_proj = delta_c;
105105
106- // Handle complementarity constraints for each lambda-eta pair
107- for (int i = 0 ; i < m_; ++i) {
108- double w_eta = std::abs (U (n_ + m_ + k_ + i, n_ + m_ + k_ + i));
109- double w_lambda = std::abs (U (n_ + i, n_ + i));
110-
111- double lambda_val = delta_c (n_ + i);
112- double eta_val = delta_c (n_ + m_ + k_ + i);
113-
114- if (lambda_val <= 0 ) {
115- delta_proj (n_ + i) = 0 ;
116- delta_proj (n_ + m_ + k_ + i) = std::max (0.0 , eta_val);
117- } else {
118- if (eta_val <= 0 ) {
119- delta_proj (n_ + i) = lambda_val;
120- delta_proj (n_ + m_ + k_ + i) = 0 ;
121- } else {
122- // If point (lambda, eta) is above the slope sqrt(w_lambda/w_eta), set
123- // lambda to 0 and keep eta Otherwise, set lambda to lambda and set eta
124- // to 0
125- if (eta_val * std::sqrt (w_eta) > lambda_val * std::sqrt (w_lambda)) {
126- delta_proj (n_ + i) = 0 ;
127- delta_proj (n_ + m_ + k_ + i) = eta_val;
128- } else {
129- delta_proj (n_ + i) = lambda_val;
130- delta_proj (n_ + m_ + k_ + i) = 0 ;
131- }
132- }
133- }
134- }
106+ // Extract the weight vectors for lambda and eta from the diagonal of the cost
107+ // matrix U.
108+ // Use absolute values to ensure numerical safety when taking square roots,
109+ // in case the user inadvertently supplies negative weights.
110+ VectorXd w_eta_vec =
111+ U.block (n_ + m_ + k_, n_ + m_ + k_, m_, m_).diagonal ().cwiseAbs ();
112+ VectorXd w_lambda_vec = U.block (n_, n_, m_, m_).diagonal ().cwiseAbs ();
113+
114+ VectorXd lambda_c = delta_c.segment (n_, m_);
115+ VectorXd eta_c = delta_c.segment (n_ + m_ + k_, m_);
116+
117+ // Set the smaller of lambda and eta to zero
118+ Eigen::Array<bool , Eigen::Dynamic, 1 > eta_larger =
119+ eta_c.array () * w_eta_vec.array ().sqrt () >
120+ lambda_c.array () * w_lambda_vec.array ().sqrt ();
121+
122+ delta_proj.segment (n_, m_) = eta_larger.select (VectorXd::Zero (m_), lambda_c);
123+ delta_proj.segment (n_ + m_ + k_, m_) =
124+ eta_larger.select (eta_c, VectorXd::Zero (m_));
125+
126+ // Clip lambda and eta at 0
127+ delta_proj.segment (n_, m_) = delta_proj.segment (n_, m_).cwiseMax (0 );
128+ delta_proj.segment (n_ + m_ + k_, m_) =
129+ delta_proj.segment (n_ + m_ + k_, m_).cwiseMax (0 );
135130
136131 return delta_proj;
137132}
0 commit comments