adding weights through the code (linear system and resids)

Author: cnpetra

src/Optimization/KktLinSysLowRank.cpp

@@ -100,6 +100,7 @@ bool KktLinSysLowRank::update(const hiopIterate* iter,
   Dx_->setToZero();
   Dx_->axdzpy_w_pattern(1.0, *iter_->zl, *iter_->sxl, nlp_->get_ixl());
   Dx_->axdzpy_w_pattern(1.0, *iter_->zu, *iter_->sxu, nlp_->get_ixu());
+  Dx_->componentMult(*nlp_->inner_prod()->M_lumped());
   nlp_->log->write("Dx in KKT", *Dx_, hovMatrices);
 
   hess_low_rank->update_logbar_diag(*Dx_);
@@ -218,7 +219,7 @@ bool KktLinSysLowRank::solveCompressed(hiopVector& rx,
   // some outputing
   nlp_->log->write("KKT Low rank: solve compressed SOL", hovIteration);
   nlp_->log->write("  dx: ", dx, hovIteration);
-  nlp_->log->write(" dyc: ", dyc, hovIteration);
+  nlp_->log->write(" dyc: ", dyc, hovWarning);
   nlp_->log->write(" dyd: ", dyd, hovIteration);
   delete r;
 #endif

src/Optimization/hiopIterate.cpp

@@ -621,7 +621,8 @@ void hiopIterate::addLinearDampingTermToGrad_x(const double& mu,
   assert(x->get_local_size() == grad_x.get_local_size());
 
   const double ct = kappa_d * mu * beta;
-  grad_x.addLinearDampingTerm(nlp->get_ixl(), nlp->get_ixu(), 1.0, ct);
+  //grad_x.addLinearDampingTerm(nlp->get_ixl(), nlp->get_ixu(), 1.0, ct);
+  nlp->inner_prod()->add_linear_damping_term(nlp->get_ixl(), nlp->get_ixu(), ct, grad_x);
 }
 
 void hiopIterate::addLinearDampingTermToGrad_d(const double& mu,

src/Optimization/hiopKKTLinSys.cpp

@@ -568,7 +568,7 @@ bool hiopKKTLinSysCompressedXYcYd::update(const hiopIterate* iter,
   Dx_->axdzpy_w_pattern(1.0, *iter_->zl, *iter_->sxl, nlp_->get_ixl());
   Dx_->axdzpy_w_pattern(1.0, *iter_->zu, *iter_->sxu, nlp_->get_ixu());
   nlp_->log->write("Dx in KKT", *Dx_, hovMatrices);
-
+  assert(false);
   // Dd=(Sdl)^{-1}Vu + (Sdu)^{-1}Vu
   Dd_->setToZero();
   Dd_->axdzpy_w_pattern(1.0, *iter_->vl, *iter_->sdl, nlp_->get_idl());
@@ -595,31 +595,58 @@ bool hiopKKTLinSysCompressedXYcYd::computeDirections(const hiopResidual* resid,
 
   /***********************************************************************
    * perform the reduction to the compressed linear system
-   * rx_tilde  = rx+Sxl^{-1}*[rszl-Zl*rxl] - Sxu^{-1}*(rszu-Zu*rxu)
+   * rx_tilde  = rx+ M_lumped* Sxl^{-1}*[rszl-Zl*rxl] - M_lumped* Sxu^{-1}*(rszu-Zu*rxu)
    * ryd_tilde = ryd + [(Sdl^{-1}Vl+Sdu^{-1}Vu)]^{-1}*
    *                     [rd + Sdl^{-1}*(rsvl-Vl*rdl)-Sdu^{-1}(rsvu-Vu*rdu)]
    * rd_tilde = rd + Sdl^{-1}*(rsvl-Vl*rdl)-Sdu^{-1}(rsvu-Vu*rdu)
    * Dd_inv = [(Sdl^{-1}Vl+Sdu^{-1}Vu)]^{-1}
    * yd_tilde = ryd + Dd_inv*rd_tilde
    */
   rx_tilde_->copyFrom(*r.rx);
+  
+  // if(nlp_->n_low_local() > 0) {
+  //   // rl:=rszl-Zl*rxl (using dir->x as working buffer)
+  //   hiopVector& rl = *(dir->x);  // temporary working buffer
+  //   rl.copyFrom(*r.rszl);
+  //   rl.axzpy(-1.0, *iter_->zl, *r.rxl);
+  //   // rx_tilde = rx+Sxl^{-1}*rl
+  //   rx_tilde_->axdzpy_w_pattern(1.0, rl, *iter_->sxl, nlp_->get_ixl());
+  // }
+  // if(nlp_->n_upp_local() > 0) {
+  //   // ru:=rszu-Zu*rxu (using dir->x as working buffer)
+  //   hiopVector& ru = *(dir->x);  // temporary working buffer
+  //   ru.copyFrom(*r.rszu);
+  //   ru.axzpy(-1.0, *iter_->zu, *r.rxu);
+  //   // rx_tilde = rx_tilde - Sxu^{-1}*ru
+  //   rx_tilde_->axdzpy_w_pattern(-1.0, ru, *iter_->sxu, nlp_->get_ixu());
+  // }
+
+  hiopVector& rx2 = *(dir->sxl);  // temporary working buffer
+  rx2.setToZero();
   if(nlp_->n_low_local() > 0) {
-    // rl:=rszl-Zl*rxl (using dir->x as working buffer)
     hiopVector& rl = *(dir->x);  // temporary working buffer
     rl.copyFrom(*r.rszl);
     rl.axzpy(-1.0, *iter_->zl, *r.rxl);
-    // rx_tilde = rx+Sxl^{-1}*rl
-    rx_tilde_->axdzpy_w_pattern(1.0, rl, *iter_->sxl, nlp_->get_ixl());
+    //// rx_tilde = rx+Sxl^{-1}*rl
+    //rx_tilde_->axdzpy_w_pattern(1.0, rl, *iter_->sxl, nlp_->get_ixl());
+    rx2.axdzpy_w_pattern(1.0, rl, *iter_->sxl, nlp_->get_ixl()); 
   }
+  
   if(nlp_->n_upp_local() > 0) {
-    // ru:=rszu-Zu*rxu (using dir->x as working buffer)
-    hiopVector& ru = *(dir->x);  // temporary working buffer
+    // ru:=rszu-Zu*rxu
+    hiopVector& ru = *(dir->x);
     ru.copyFrom(*r.rszu);
     ru.axzpy(-1.0, *iter_->zu, *r.rxu);
-    // rx_tilde = rx_tilde - Sxu^{-1}*ru
-    rx_tilde_->axdzpy_w_pattern(-1.0, ru, *iter_->sxu, nlp_->get_ixu());
+    //// rx_tilde = rx_tilde - Sxu^{-1}*ru
+    ////rx_tilde_->axdzpy_w_pattern(-1.0, ru, *iter_->sxu, nlp_->get_ixu());
+    rx2.axdzpy_w_pattern(-1.0, ru, *iter_->sxu, nlp_->get_ixu());
+  }
+  if(nlp_->n_low_local() > 0 || nlp_->n_upp_local() > 0) {
+    rx2.componentMult(*nlp_->inner_prod()->M_lumped());
+    rx_tilde_->axpy(1.0, rx2);
   }
 
+  
   // for ryd_tilde:
   ryd_tilde_->copyFrom(*r.ryd);
   // 1. the diag (Sdl^{-1}Vl+Sdu^{-1}Vu)^{-1} has already computed in Dd_inv in 'update'

src/Optimization/hiopResidual.cpp

@@ -146,7 +146,7 @@ double hiopResidual::compute_nlp_infeasib_onenorm(const hiopIterate& it, const h
     rdu->axpy(-1.0, *it.d);
     rdu->selectPattern(nlp->get_idu());
   }
-
+  //! the ifs above seem to do nothing for the return value
   nlp->runStats.tmSolverInternal.stop();
   return nrmOne_infeasib;
 }
@@ -192,6 +192,8 @@ int hiopResidual::update(const hiopIterate& it,
   jac_d.transTimesVec(1.0, *rx, 1.0, *it.yd);
   rx->axpy(1.0, grad);
 
+  //nlp->log->write("rx w/o log terms:", *rx, hovWarning);
+  
   //buf = rx->infnorm_local();
   buf = nlp->inner_prod()->norm_stationarity(*rx);
   nrmInf_nlp_optim = fmax(nrmInf_nlp_optim, buf);
@@ -205,6 +207,8 @@ int hiopResidual::update(const hiopIterate& it,
   nrmInf_bar_optim = fmax(nrmInf_bar_optim, nlp->inner_prod()->norm_stationarity(*rx));
   nrmOne_bar_optim += nlp->inner_prod()->norm_M_one(*rx);
 
+  //nlp->log->write("rx with log terms:", *rx, hovWarning);
+  
   //~ done with rx
   // rd
   rd->copyFrom(*it.yd);
@@ -470,7 +474,7 @@ void hiopResidual::update_soc(const hiopIterate& it,
   nrmInf_bar_optim = nrmInf_bar_feasib = nrmInf_bar_complem = 0;
   nrmOne_nlp_feasib = nrmOne_bar_feasib = 0.;
   nrmOne_nlp_optim = nrmOne_bar_optim = 0.;
-
+  assert(false);
   size_type nx_loc = rx->get_local_size();
   const double& mu = logprob.mu;
   double buf;