diff --git a/src/Optimization/hiopHessianLowRank.cpp b/src/Optimization/hiopHessianLowRank.cpp index 4f03fa336..38e501b1e 100644 --- a/src/Optimization/hiopHessianLowRank.cpp +++ b/src/Optimization/hiopHessianLowRank.cpp @@ -46,6 +46,13 @@ // Lawrence Livermore National Security, LLC, and shall not be used for advertising or // product endorsement purposes. +/** + * @file hiopHessianLowRank.cpp + * + * @author Cosmin G. Petra , LLNL + * + */ + #include "hiopHessianLowRank.hpp" #include "LinAlgFactory.hpp" #include "hiopVectorPar.hpp" @@ -75,123 +82,145 @@ using namespace std; namespace hiop { -hiopHessianLowRank::hiopHessianLowRank(hiopNlpDenseConstraints* nlp_, int max_mem_len) - : l_max(max_mem_len), l_curr(-1), sigma(1.), sigma0(1.), nlp(nlp_), matrixChanged(false) +hiopHessianLowRank::hiopHessianLowRank(hiopNlpDenseConstraints* nlp, int max_mem_len) + : l_max_(max_mem_len), + l_curr_(-1), + sigma_(1.), + sigma0_(1.), + nlp_(nlp), + matrixChanged_(false) { - DhInv = nlp->alloc_primal_vec(); - St = nlp->alloc_multivector_primal(0,l_max); - Yt = St->alloc_clone(); //faster than nlp->alloc_multivector_primal(...); + DhInv_ = nlp->alloc_primal_vec(); + St_ = nlp->alloc_multivector_primal(0, l_max_); + Yt_ = St_->alloc_clone(); //faster than nlp->alloc_multivector_primal(...); //these are local - L = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); - D = LinearAlgebraFactory::create_vector("DEFAULT", 0); - V = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); + L_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); + D_ = LinearAlgebraFactory::create_vector("DEFAULT", 0); + V_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); - //the previous iteration's objects are set to NULL - _it_prev=NULL; _grad_f_prev=NULL; _Jac_c_prev=NULL; _Jac_d_prev=NULL; + //the previous iteration's objects are set to nullptr + it_prev_ = nullptr; + grad_f_prev_ = nullptr; + Jac_c_prev_ = nullptr; + Jac_d_prev_ = nullptr; //internal buffers for memory pool (none of them should be in n) #ifdef HIOP_USE_MPI - _buff_kxk = new double[nlp->m() * nlp->m()]; - _buff_2lxk = new double[nlp->m() * 2*l_max]; - _buff1_lxlx3 = new double[3*l_max*l_max]; - _buff2_lxlx3 = new double[3*l_max*l_max]; + buff_kxk_ = new double[nlp->m() * nlp->m()]; + buff_2lxk_ = new double[nlp->m() * 2 * l_max_]; + buff1_lxlx3_ = new double[3 * l_max_ * l_max_]; + buff2_lxlx3_ = new double[3 * l_max_ * l_max_]; #else //not needed in non-MPI mode - _buff_kxk = NULL; - _buff_2lxk = NULL; - _buff1_lxlx3 = _buff2_lxlx3 = NULL; + buff_kxk_ = nullptr; + buff_2lxk_ = nullptr; + buff1_lxlx3_ = nullptr; + buff2_lxlx3_ = nullptr; #endif //auxiliary objects/buffers - _S1=_Y1=NULL; - _lxl_mat1=_kxl_mat1=_kx2l_mat1=NULL; - _l_vec1 = _l_vec2 = _2l_vec1 = NULL; - _n_vec1 = DhInv->alloc_clone(); - _n_vec2 = DhInv->alloc_clone(); - - _V_work_vec=LinearAlgebraFactory::create_vector("DEFAULT", 0); - _V_ipiv_vec=NULL; _V_ipiv_size=-1; - - - sigma0 = nlp->options->GetNumeric("sigma0"); - sigma=sigma0; + S1_ = Y1_ = nullptr; + lxl_mat1_ = nullptr; + kxl_mat1_ = nullptr; + kx2l_mat1_ = nullptr; + l_vec1_ = nullptr; + l_vec2_ = nullptr; + twol_vec1_ = nullptr; + n_vec1_ = DhInv_->alloc_clone(); + n_vec2_ = DhInv_->alloc_clone(); + + V_work_vec_ = LinearAlgebraFactory::create_vector("DEFAULT", 0); + V_ipiv_vec_ = nullptr; + V_ipiv_size_ = -1; + + sigma0_ = nlp->options->GetNumeric("sigma0"); + sigma_ = sigma0_; string sigma_strategy = nlp->options->GetString("sigma_update_strategy"); transform(sigma_strategy.begin(), sigma_strategy.end(), sigma_strategy.begin(), ::tolower); - sigma_update_strategy = SIGMA_STRATEGY3; - if(sigma_strategy=="sty") - sigma_update_strategy=SIGMA_STRATEGY1; - else if(sigma_strategy=="sty_inv") - sigma_update_strategy=SIGMA_STRATEGY2; - else if(sigma_strategy=="snrm_ynrm") - sigma_update_strategy=SIGMA_STRATEGY3; - else if(sigma_strategy=="sty_srnm_ynrm") - sigma_update_strategy=SIGMA_STRATEGY4; - else if(sigma_strategy=="sigma0") - sigma_update_strategy=SIGMA_CONSTANT; - else assert(false && "sigma_update_strategy option not recognized"); - - sigma_safe_min=1e-8; - sigma_safe_max=1e+8; - nlp->log->printf(hovScalars, "Hessian Low Rank: initial sigma is %g\n", sigma); - nlp->log->printf(hovScalars, "Hessian Low Rank: sigma update strategy is %d [%s]\n", sigma_update_strategy, sigma_strategy.c_str()); - - _Dx = DhInv->alloc_clone(); + sigma_update_strategy_ = SIGMA_STRATEGY3; + if(sigma_strategy=="sty") { + sigma_update_strategy_ = SIGMA_STRATEGY1; + } else if(sigma_strategy=="sty_inv") { + sigma_update_strategy_ = SIGMA_STRATEGY2; + } else if(sigma_strategy=="snrm_ynrm") { + sigma_update_strategy_ = SIGMA_STRATEGY3; + } else if(sigma_strategy=="sty_srnm_ynrm") { + sigma_update_strategy_ = SIGMA_STRATEGY4; + } else if(sigma_strategy=="sigma0") { + sigma_update_strategy_ = SIGMA_CONSTANT; + } else { + assert(false && "sigma_update_strategy option not recognized"); + } + + sigma_safe_min_ = 1e-8; + sigma_safe_max_ = 1e+8; + nlp->log->printf(hovScalars, "Hessian Low Rank: initial sigma is %g\n", sigma_); + nlp->log->printf(hovScalars, "Hessian Low Rank: sigma update strategy is %d [%s]\n", sigma_update_strategy_, sigma_strategy.c_str()); + + Dx_ = DhInv_->alloc_clone(); #ifdef HIOP_DEEPCHECKS - _Vmat = V->alloc_clone(); + Vmat_ = V_->alloc_clone(); #endif - yk = nlp->alloc_primal_vec(); - sk = nlp->alloc_primal_vec(); - + yk_ = nlp->alloc_primal_vec(); + sk_ = nlp->alloc_primal_vec(); } hiopHessianLowRank::~hiopHessianLowRank() { - if(DhInv) delete DhInv; - delete _Dx; - - if(St) delete St; - if(Yt) delete Yt; - if(L) delete L; - if(D) delete D; - if(V) delete V; - if(yk) delete yk; - if(sk) delete sk; + delete DhInv_; + delete Dx_; + delete St_; + delete Yt_; + delete L_; + delete D_; + delete V_; + delete yk_; + delete sk_; #ifdef HIOP_DEEPCHECKS - delete _Vmat; + delete Vmat_; #endif + delete it_prev_; + delete grad_f_prev_; + delete Jac_c_prev_; + delete Jac_d_prev_; - if(_it_prev) delete _it_prev; - if(_grad_f_prev)delete _grad_f_prev; - if(_Jac_c_prev) delete _Jac_c_prev; - if(_Jac_d_prev) delete _Jac_d_prev; - - if(_buff_kxk) delete[] _buff_kxk; - if(_buff_2lxk) delete[] _buff_2lxk; - if(_buff1_lxlx3) delete[] _buff1_lxlx3; - if(_buff2_lxlx3) delete[] _buff2_lxlx3; - - if(_S1) delete _S1; - if(_Y1) delete _Y1; - if(_lxl_mat1) delete _lxl_mat1; - if(_kxl_mat1) delete _kxl_mat1; - if(_kx2l_mat1) delete _kx2l_mat1; + if(buff_kxk_) { + delete[] buff_kxk_; + } + if(buff_2lxk_) { + delete[] buff_2lxk_; + } + if(buff1_lxlx3_) { + delete[] buff1_lxlx3_; + } + if(buff2_lxlx3_) { + delete[] buff2_lxlx3_; + } - if(_l_vec1) delete _l_vec1; - if(_l_vec2) delete _l_vec2; - if(_n_vec1) delete _n_vec1; - if(_n_vec2) delete _n_vec2; - if(_2l_vec1) delete _2l_vec1; - if(_V_ipiv_vec) delete[] _V_ipiv_vec; - if(_V_work_vec) delete _V_work_vec; + delete S1_; + delete Y1_; + delete lxl_mat1_; + delete kxl_mat1_; + delete kx2l_mat1_; + + delete l_vec1_; + delete l_vec2_; + delete n_vec1_; + delete n_vec2_; + delete twol_vec1_; + if(V_ipiv_vec_) { + delete[] V_ipiv_vec_; + } + delete V_work_vec_; - for(auto* it: a) { + for(auto* it: a_) { delete it; } - for(auto* it: b) { + for(auto* it: b_) { delete it; } } @@ -199,15 +228,15 @@ hiopHessianLowRank::~hiopHessianLowRank() bool hiopHessianLowRank::updateLogBarrierDiagonal(const hiopVector& Dx) { - DhInv->setToConstant(sigma); - DhInv->axpy(1.0,Dx); - _Dx->copyFrom(Dx); + DhInv_->setToConstant(sigma_); + DhInv_->axpy(1.0, Dx); + Dx_->copyFrom(Dx); #ifdef HIOP_DEEPCHECKS - assert(DhInv->allPositive()); + assert(DhInv_->allPositive()); #endif - DhInv->invert(); - nlp->log->write("hiopHessianLowRank: inverse diag DhInv:", *DhInv, hovMatrices); - matrixChanged=true; + DhInv_->invert(); + nlp_->log->write("hiopHessianLowRank: inverse diag DhInv:", *DhInv_, hovMatrices); + matrixChanged_ = true; return true; } @@ -216,147 +245,163 @@ void hiopHessianLowRank::print(FILE* f, hiopOutVerbosity v, const char* msg) con { fprintf(f, "%s\n", msg); #ifdef HIOP_DEEPCHECKS - nlp->log->write("Dx", *_Dx, v); + nlp_->log->write("Dx", *Dx_, v); #else fprintf(f, "Dx is not stored in this class, but it can be computed from Dx=DhInv^(1)-sigma"); #endif - nlp->log->printf(v, "sigma=%22.16f;\n", sigma); - nlp->log->write("DhInv", *DhInv, v); - nlp->log->write("S_trans", *St, v); - nlp->log->write("Y_trans", *Yt, v); + nlp_->log->printf(v, "sigma=%22.16f;\n", sigma_); + nlp_->log->write("DhInv_", *DhInv_, v); + nlp_->log->write("S_trans", *St_, v); + nlp_->log->write("Y_trans", *Yt_, v); fprintf(f, " [[Internal representation]]\n"); #ifdef HIOP_DEEPCHECKS - nlp->log->write("V", *_Vmat, v); + nlp_->log->write("V", *Vmat_, v); #else fprintf(f, "V matrix is available at this point (only its LAPACK factorization). Print it in updateInternalBFGSRepresentation() instead, before factorizeV()\n"); #endif - nlp->log->write("L", *L, v); - nlp->log->write("D", *D, v); + nlp_->log->write("L", *L_, v); + nlp_->log->write("D", *D_, v); } #endif #include -bool hiopHessianLowRank::update(const hiopIterate& it_curr, const hiopVector& grad_f_curr_, - const hiopMatrix& Jac_c_curr_, const hiopMatrix& Jac_d_curr_) +bool hiopHessianLowRank::update(const hiopIterate& it_curr, + const hiopVector& grad_f_curr, + const hiopMatrix& Jac_c_curr_in, + const hiopMatrix& Jac_d_curr_in) { - nlp->runStats.tmSolverInternal.start(); + nlp_->runStats.tmSolverInternal.start(); - const hiopMatrixDense& Jac_c_curr = dynamic_cast(Jac_c_curr_); - const hiopMatrixDense& Jac_d_curr = dynamic_cast(Jac_d_curr_); + const hiopMatrixDense& Jac_c_curr = dynamic_cast(Jac_c_curr_in); + const hiopMatrixDense& Jac_d_curr = dynamic_cast(Jac_d_curr_in); #ifdef HIOP_DEEPCHECKS - assert(it_curr.zl->matchesPattern(nlp->get_ixl())); - assert(it_curr.zu->matchesPattern(nlp->get_ixu())); - assert(it_curr.sxl->matchesPattern(nlp->get_ixl())); - assert(it_curr.sxu->matchesPattern(nlp->get_ixu())); + assert(it_curr.zl->matchesPattern(nlp_->get_ixl())); + assert(it_curr.zu->matchesPattern(nlp_->get_ixu())); + assert(it_curr.sxl->matchesPattern(nlp_->get_ixl())); + assert(it_curr.sxu->matchesPattern(nlp_->get_ixu())); #endif - //on first call l_curr=-1 - if(l_curr>=0) { - size_type n=grad_f_curr_.get_size(); + //on first call l_curr_ = -1 + if(l_curr_ >= 0) { + size_type n = grad_f_curr.get_size(); //compute s_new = x_curr-x_prev - hiopVector& s_new = new_n_vec1(n); s_new.copyFrom(*it_curr.x); s_new.axpy(-1.,*_it_prev->x); - double s_infnorm=s_new.infnorm(); - if(s_infnorm>=100*std::numeric_limits::epsilon()) { //norm of s not too small + hiopVector& s_new = new_n_vec1(n); + s_new.copyFrom(*it_curr.x); + s_new.axpy(-1., *it_prev_->x); + double s_infnorm = s_new.infnorm(); + if(s_infnorm >= 100 * std::numeric_limits::epsilon()) { //norm of s not too small //compute y_new = \grad J(x_curr,\lambda_curr) - \grad J(x_prev, \lambda_curr) (yes, J(x_prev, \lambda_curr)) // = graf_f_curr-grad_f_prev + (Jac_c_curr-Jac_c_prev)yc_curr+ (Jac_d_curr-Jac_c_prev)yd_curr - zl_curr*s_new + zu_curr*s_new hiopVector& y_new = new_n_vec2(n); - y_new.copyFrom(grad_f_curr_); - y_new.axpy(-1., *_grad_f_prev); + y_new.copyFrom(grad_f_curr); + y_new.axpy(-1., *grad_f_prev_); Jac_c_curr.transTimesVec (1.0, y_new, 1.0, *it_curr.yc); - _Jac_c_prev->transTimesVec(1.0, y_new,-1.0, *it_curr.yc); //!opt if nlp->Jac_c_isLinear no need for the multiplications + Jac_c_prev_->transTimesVec(1.0, y_new,-1.0, *it_curr.yc); //!opt if nlp_->Jac_c_isLinear no need for the multiplications Jac_d_curr.transTimesVec (1.0, y_new, 1.0, *it_curr.yd); //!opt same here - _Jac_d_prev->transTimesVec(1.0, y_new,-1.0, *it_curr.yd); + Jac_d_prev_->transTimesVec(1.0, y_new,-1.0, *it_curr.yd); - double sTy = s_new.dotProductWith(y_new), s_nrm2=s_new.twonorm(), y_nrm2=y_new.twonorm(); + double sTy = s_new.dotProductWith(y_new); + double s_nrm2 = s_new.twonorm(); + double y_nrm2 = y_new.twonorm(); #ifdef HIOP_DEEPCHECKS - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianLowRank: s^T*y=%20.14e ||s||=%20.14e ||y||=%20.14e\n", sTy, s_nrm2, y_nrm2); - nlp->log->write("hiopHessianLowRank s_new",s_new, hovIteration); - nlp->log->write("hiopHessianLowRank y_new",y_new, hovIteration); + nlp_->log->printf(hovLinAlgScalarsVerb, "hiopHessianLowRank: s^T*y=%20.14e ||s||=%20.14e ||y||=%20.14e\n", sTy, s_nrm2, y_nrm2); + nlp_->log->write("hiopHessianLowRank s_new", s_new, hovIteration); + nlp_->log->write("hiopHessianLowRank y_new", y_new, hovIteration); #endif - if(sTy>s_nrm2*y_nrm2*sqrt(std::numeric_limits::epsilon())) { //sTy far away from zero + if(sTy > s_nrm2 * y_nrm2 * sqrt(std::numeric_limits::epsilon())) { //sTy far away from zero - if(l_max>0) { + if(l_max_ > 0) { //compute the new row in L, update S and Y (either augment them or shift cols and add s_new and y_new) - hiopVector& YTs = new_l_vec1(l_curr); - Yt->timesVec(0.0, YTs, 1.0, s_new); + hiopVector& YTs = new_l_vec1(l_curr_); + Yt_->timesVec(0.0, YTs, 1.0, s_new); //update representation - if(l_currappendRow(s_new); - Yt->appendRow(y_new); - growL(l_curr, l_max, YTs); - growD(l_curr, l_max, sTy); - l_curr++; + St_->appendRow(s_new); + Yt_->appendRow(y_new); + growL(l_curr_, l_max_, YTs); + growD(l_curr_, l_max_, sTy); + l_curr_++; } else { //shift - St->shiftRows(-1); - Yt->shiftRows(-1); - St->replaceRow(l_max-1, s_new); - Yt->replaceRow(l_max-1, y_new); - updateL(YTs,sTy); + St_->shiftRows(-1); + Yt_->shiftRows(-1); + St_->replaceRow(l_max_-1, s_new); + Yt_->replaceRow(l_max_-1, y_new); + updateL(YTs, sTy); updateD(sTy); - l_curr=l_max; + l_curr_ = l_max_; } - } //end of l_max>0 + } //end of l_max_>0 #ifdef HIOP_DEEPCHECKS - nlp->log->printf(hovMatrices, "\nhiopHessianLowRank: these are L and D from the BFGS compact representation\n"); - nlp->log->write("L", *L, hovMatrices); - nlp->log->write("D", *D, hovMatrices); - nlp->log->printf(hovMatrices, "\n"); + nlp_->log->printf(hovMatrices, "\nhiopHessianLowRank: these are L and D from the BFGS compact representation\n"); + nlp_->log->write("L", *L_, hovMatrices); + nlp_->log->write("D", *D_, hovMatrices); + nlp_->log->printf(hovMatrices, "\n"); #endif //update B0 (i.e., sigma) - switch (sigma_update_strategy ) { - case SIGMA_STRATEGY1: - sigma=sTy/(s_nrm2*s_nrm2); - break; - case SIGMA_STRATEGY2: - sigma=y_nrm2*y_nrm2/sTy; - break; - case SIGMA_STRATEGY3: - sigma=sqrt(s_nrm2*s_nrm2 / y_nrm2 / y_nrm2); - break; - case SIGMA_STRATEGY4: - sigma=0.5*(sTy/(s_nrm2*s_nrm2)+y_nrm2*y_nrm2/sTy); - break; - case SIGMA_CONSTANT: - sigma=sigma0; - break; - default: - assert(false && "Option value for sigma_update_strategy was not recognized."); - break; + switch (sigma_update_strategy_ ) { + case SIGMA_STRATEGY1: + sigma_ = sTy / (s_nrm2 * s_nrm2); + break; + case SIGMA_STRATEGY2: + sigma_ = y_nrm2 * y_nrm2 / sTy; + break; + case SIGMA_STRATEGY3: + sigma_ = sqrt(s_nrm2 * s_nrm2 / y_nrm2 / y_nrm2); + break; + case SIGMA_STRATEGY4: + sigma_ = 0.5*(sTy / (s_nrm2 * s_nrm2) + y_nrm2 * y_nrm2 / sTy); + break; + case SIGMA_CONSTANT: + sigma_ = sigma0_; + break; + default: + assert(false && "Option value for sigma_update_strategy was not recognized."); + break; } // else of the switch //safe guard it - sigma=fmax(fmin(sigma_safe_max, sigma), sigma_safe_min); - nlp->log->printf(hovLinAlgScalars, "hiopHessianLowRank: sigma was updated to %22.16e\n", sigma); + sigma_ = fmax(fmin(sigma_safe_max_, sigma_), sigma_safe_min_); + nlp_->log->printf(hovLinAlgScalars, "hiopHessianLowRank: sigma was updated to %22.16e\n", sigma_); } else { //sTy is too small or negative -> skip - nlp->log->printf(hovLinAlgScalars, "hiopHessianLowRank: s^T*y=%12.6e not positive enough... skipping the Hessian update\n", sTy); + nlp_->log->printf(hovLinAlgScalars, "hiopHessianLowRank: s^T*y=%12.6e not positive enough... skipping the Hessian update\n", sTy); } } else {// norm of s_new is too small -> skip - nlp->log->printf(hovLinAlgScalars, "hiopHessianLowRank: ||s_new||=%12.6e too small... skipping the Hessian update\n", s_infnorm); + nlp_->log->printf(hovLinAlgScalars, "hiopHessianLowRank: ||s_new||=%12.6e too small... skipping the Hessian update\n", s_infnorm); } //save this stuff for next update - _it_prev->copyFrom(it_curr); _grad_f_prev->copyFrom(grad_f_curr_); - _Jac_c_prev->copyFrom(Jac_c_curr); _Jac_d_prev->copyFrom(Jac_d_curr); - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianLowRank: storing the iteration info as 'previous'\n", s_infnorm); + it_prev_->copyFrom(it_curr); + grad_f_prev_->copyFrom(grad_f_curr); + Jac_c_prev_->copyFrom(Jac_c_curr); + Jac_d_prev_->copyFrom(Jac_d_curr); + nlp_->log->printf(hovLinAlgScalarsVerb, "hiopHessianLowRank: storing the iteration info as 'previous'\n", s_infnorm); } else { //this is the first optimization iterate, just save the iterate and exit - if(NULL==_it_prev) _it_prev = it_curr.new_copy(); - if(NULL==_grad_f_prev) _grad_f_prev = grad_f_curr_.new_copy(); - if(NULL==_Jac_c_prev) _Jac_c_prev = Jac_c_curr.new_copy(); - if(NULL==_Jac_d_prev) _Jac_d_prev = Jac_d_curr.new_copy(); + if(nullptr == it_prev_) { + it_prev_ = it_curr.new_copy(); + } + if(nullptr == grad_f_prev_) { + grad_f_prev_ = grad_f_curr.new_copy(); + } + if(nullptr == Jac_c_prev_) { + Jac_c_prev_ = Jac_c_curr.new_copy(); + } + if(nullptr == Jac_d_prev_) { + Jac_d_prev_ = Jac_d_curr.new_copy(); + } - nlp->log->printf(hovLinAlgScalarsVerb, "HessianLowRank on first update, just saving iteration\n"); + nlp_->log->printf(hovLinAlgScalarsVerb, "HessianLowRank on first update, just saving iteration\n"); - l_curr++; + l_curr_++; } - nlp->runStats.tmSolverInternal.stop(); + nlp_->runStats.tmSolverInternal.stop(); return true; } @@ -372,79 +417,91 @@ bool hiopHessianLowRank::update(const hiopIterate& it_curr, const hiopVector& gr */ void hiopHessianLowRank::updateInternalBFGSRepresentation() { - size_type n=St->n(), l=St->m(); + size_type n = St_->n(); + size_type l = St_->m(); //grow L,D, andV if needed - if(L->m()!=l) { delete L; L=LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, l);} - if(D->get_size()!=l) { delete D; D=LinearAlgebraFactory::create_vector("DEFAULT", l); } - if(V->m()!=2*l) {delete V; V=LinearAlgebraFactory::create_matrix_dense("DEFAULT", 2*l, 2*l); } + if(L_->m() != l) { + delete L_; + L_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, l); + } + if(D_->get_size() != l) { + delete D_; + D_ = LinearAlgebraFactory::create_vector("DEFAULT", l); + } + if(V_->m() != 2 * l) { + delete V_; + V_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 2 * l, 2 * l); + } //-- block (2,2) hiopMatrixDense& DpYtDhInvY = new_lxl_mat1(l); - symmMatTimesDiagTimesMatTrans_local(0.0, DpYtDhInvY, 1.0,*Yt,*DhInv); + symmMatTimesDiagTimesMatTrans_local(0.0, DpYtDhInvY, 1.0, *Yt_, *DhInv_); #ifdef HIOP_USE_MPI - const size_t buffsize=l*l*sizeof(double); - memcpy(_buff1_lxlx3, DpYtDhInvY.local_data(), buffsize); + const size_t buffsize = l * l * sizeof(double); + memcpy(buff1_lxlx3_, DpYtDhInvY.local_data(), buffsize); #else - DpYtDhInvY.addDiagonal(1., *D); - V->copyBlockFromMatrix(l,l,DpYtDhInvY); + DpYtDhInvY.addDiagonal(1., *D_); + V_->copyBlockFromMatrix(l, l, DpYtDhInvY); #endif //-- block (1,2) hiopMatrixDense& StB0DhInvYmL = DpYtDhInvY; //just a rename hiopVector& B0DhInv = new_n_vec1(n); - B0DhInv.copyFrom(*DhInv); B0DhInv.scale(sigma); - matTimesDiagTimesMatTrans_local(StB0DhInvYmL, *St, B0DhInv, *Yt); + B0DhInv.copyFrom(*DhInv_); + B0DhInv.scale(sigma_); + matTimesDiagTimesMatTrans_local(StB0DhInvYmL, *St_, B0DhInv, *Yt_); #ifdef HIOP_USE_MPI - memcpy(_buff1_lxlx3+l*l, StB0DhInvYmL.local_data(), buffsize); + memcpy(buff1_lxlx3_ + l * l, StB0DhInvYmL.local_data(), buffsize); #else //substract L - StB0DhInvYmL.addMatrix(-1.0, *L); + StB0DhInvYmL.addMatrix(-1.0, *L_); // (1,2) block in V - V->copyBlockFromMatrix(0,l,StB0DhInvYmL); + V_->copyBlockFromMatrix(0, l, StB0DhInvYmL); #endif //-- block (2,2) hiopVector& theDiag = B0DhInv; //just a rename, also reuses values theDiag.addConstant(-1.0); //at this point theDiag=DhInv*B0-I - theDiag.scale(sigma); + theDiag.scale(sigma_); hiopMatrixDense& StDS = DpYtDhInvY; //a rename - symmMatTimesDiagTimesMatTrans_local(0.0, StDS, 1.0, *St, theDiag); + symmMatTimesDiagTimesMatTrans_local(0.0, StDS, 1.0, *St_, theDiag); #ifdef HIOP_USE_MPI - memcpy(_buff1_lxlx3+2*l*l, DpYtDhInvY.local_data(), buffsize); + memcpy(buff1_lxlx3_+2*l*l, DpYtDhInvY.local_data(), buffsize); #else - V->copyBlockFromMatrix(0,0,StDS); + V_->copyBlockFromMatrix(0, 0, StDS); #endif #ifdef HIOP_USE_MPI int ierr; - ierr = MPI_Allreduce(_buff1_lxlx3, _buff2_lxlx3, 3*l*l, MPI_DOUBLE, MPI_SUM, nlp->get_comm()); assert(ierr==MPI_SUCCESS); + ierr = MPI_Allreduce(buff1_lxlx3_, buff2_lxlx3_, 3 * l * l, MPI_DOUBLE, MPI_SUM, nlp_->get_comm()); + assert(ierr == MPI_SUCCESS); // - block (2,2) - DpYtDhInvY.copyFrom(_buff2_lxlx3); - DpYtDhInvY.addDiagonal(1., *D); - V->copyBlockFromMatrix(l,l,DpYtDhInvY); + DpYtDhInvY.copyFrom(buff2_lxlx3_); + DpYtDhInvY.addDiagonal(1., *D_); + V_->copyBlockFromMatrix(l, l, DpYtDhInvY); // - block (1,2) - StB0DhInvYmL.copyFrom(_buff2_lxlx3+l*l); - StB0DhInvYmL.addMatrix(-1.0, *L); - V->copyBlockFromMatrix(0,l,StB0DhInvYmL); + StB0DhInvYmL.copyFrom(buff2_lxlx3_ + l * l); + StB0DhInvYmL.addMatrix(-1.0, *L_); + V_->copyBlockFromMatrix(0, l, StB0DhInvYmL); // - block (1,1) - StDS.copyFrom(_buff2_lxlx3+2*l*l); - V->copyBlockFromMatrix(0,0,StDS); + StDS.copyFrom(buff2_lxlx3_ + 2 * l * l); + V_->copyBlockFromMatrix(0, 0, StDS); #endif #ifdef HIOP_DEEPCHECKS - delete _Vmat; - _Vmat = V->new_copy(); - _Vmat->overwriteLowerTriangleWithUpper(); + delete Vmat_; + Vmat_ = V_->new_copy(); + Vmat_->overwriteLowerTriangleWithUpper(); #endif //finally, factorize V factorizeV(); - matrixChanged=false; + matrixChanged_ = false; } /* Solves this*x = res as x = this^{-1}*res @@ -457,42 +514,48 @@ void hiopHessianLowRank::updateInternalBFGSRepresentation() */ void hiopHessianLowRank::solve(const hiopVector& rhsx, hiopVector& x) { - if(matrixChanged) updateInternalBFGSRepresentation(); + if(matrixChanged_) { + updateInternalBFGSRepresentation(); + } - size_type n=St->n(), l=St->m(); + size_type n = St_->n(); + size_type l = St_->m(); #ifdef HIOP_DEEPCHECKS - assert(rhsx.get_size()==n); - assert(x.get_size()==n); - assert(DhInv->get_size()==n); - assert(DhInv->isfinite_local() && "inf or nan entry detected"); + assert(rhsx.get_size() == n); + assert(x.get_size() == n); + assert(DhInv_->get_size() == n); + assert(DhInv_->isfinite_local() && "inf or nan entry detected"); assert(rhsx.isfinite_local() && "inf or nan entry detected in rhs"); #endif //1. x = DhInv*res x.copyFrom(rhsx); - x.componentMult(*DhInv); + x.componentMult(*DhInv_); //2. stx= S^T*B0*DhInv*res and ytx=Y^T*DhInv*res - hiopVector&stx=new_l_vec1(l), &ytx=new_l_vec2(l); - stx.setToZero(); ytx.setToZero(); - Yt->timesVec(0.0,ytx,1.0,x); + hiopVector& stx = new_l_vec1(l); + hiopVector& ytx = new_l_vec2(l); + stx.setToZero(); + ytx.setToZero(); + Yt_->timesVec(0.0, ytx, 1.0, x); hiopVector& B0DhInvx = new_n_vec1(n); B0DhInvx.copyFrom(x); //it contains DhInv*res - B0DhInvx.scale(sigma); //B0*(DhInv*res) - St->timesVec(0.0,stx,1.0,B0DhInvx); + B0DhInvx.scale(sigma_); //B0*(DhInv*res) + St_->timesVec(0.0, stx, 1.0, B0DhInvx); //3. solve with V - hiopVector& spart=stx; hiopVector& ypart=ytx; - solveWithV(spart,ypart); + hiopVector& spart = stx; + hiopVector& ypart = ytx; + solveWithV(spart, ypart); //4. multiply with DhInv*[B0*S Y], namely // result = DhInv*(B0*S*spart + Y*ypart) - hiopVector& result = new_n_vec1(n); - St->transTimesVec(0.0, result, 1.0, spart); - result.scale(sigma); - Yt->transTimesVec(1.0, result, 1.0, ypart); - result.componentMult(*DhInv); + hiopVector& result = new_n_vec1(n); + St_->transTimesVec(0.0, result, 1.0, spart); + result.scale(sigma_); + Yt_->transTimesVec(1.0, result, 1.0, ypart); + result.componentMult(*DhInv_); //5. x = first term - second term = x_computed_in_1 - result x.axpy(-1.0,result); @@ -510,53 +573,60 @@ void hiopHessianLowRank::solve(const hiopVector& rhsx, hiopVector& x) * X is kxn */ void hiopHessianLowRank:: -symMatTimesInverseTimesMatTrans(double beta, hiopMatrixDense& W, - double alpha, const hiopMatrixDense& X) +symMatTimesInverseTimesMatTrans(double beta, hiopMatrixDense& W, double alpha, const hiopMatrixDense& X) { - if(matrixChanged) updateInternalBFGSRepresentation(); + if(matrixChanged_) { + updateInternalBFGSRepresentation(); + } - size_type n=St->n(), l=St->m(); - size_type k=W.m(); - assert(X.m()==k); - assert(X.n()==n); + size_type n = St_->n(); + size_type l = St_->m(); + size_type k = W.m(); + assert(X.m() == k); + assert(X.n() == n); #ifdef HIOP_DEEPCHECKS - nlp->log->write("symMatTimesInverseTimesMatTrans: X is: ", X, hovMatrices); + nlp_->log->write("symMatTimesInverseTimesMatTrans: X is: ", X, hovMatrices); #endif //1. compute W=beta*W + alpha*X*DhInv*X' #ifdef HIOP_USE_MPI - if(0==nlp->get_rank()) - symmMatTimesDiagTimesMatTrans_local(beta,W,alpha,X,*DhInv); - else - symmMatTimesDiagTimesMatTrans_local(0.0, W,alpha,X,*DhInv); + if(0 == nlp_->get_rank()) { + symmMatTimesDiagTimesMatTrans_local(beta, W, alpha, X, *DhInv_); + } else { + symmMatTimesDiagTimesMatTrans_local(0.0, W, alpha, X, *DhInv_); + } //W will be MPI_All_reduced later #else - symmMatTimesDiagTimesMatTrans_local(beta,W,alpha,X,*DhInv); + symmMatTimesDiagTimesMatTrans_local(beta, W, alpha, X, *DhInv_); #endif //2. compute S1=X*DhInv*B0*S and Y1=X*DhInv*Y - hiopMatrixDense &S1=new_S1(X,*St), &Y1=new_Y1(X,*Yt); //both are kxl + hiopMatrixDense& S1 = new_S1(X, *St_); + hiopMatrixDense& Y1 = new_Y1(X, *Yt_); //both are kxl hiopVector& B0DhInv = new_n_vec1(n); - B0DhInv.copyFrom(*DhInv); B0DhInv.scale(sigma); - matTimesDiagTimesMatTrans_local(S1, X, B0DhInv, *St); - matTimesDiagTimesMatTrans_local(Y1, X, *DhInv, *Yt); + B0DhInv.copyFrom(*DhInv_); + B0DhInv.scale(sigma_); + matTimesDiagTimesMatTrans_local(S1, X, B0DhInv, *St_); + matTimesDiagTimesMatTrans_local(Y1, X, *DhInv_, *Yt_); //3. reduce W, S1, and Y1 (dimensions: kxk, kxl, kxl) - hiopMatrixDense& S2Y2 = new_kx2l_mat1(k,l); //Initialy S2Y2 = [Y1 S1] - S2Y2.copyBlockFromMatrix(0,0,S1); - S2Y2.copyBlockFromMatrix(0,l,Y1); + hiopMatrixDense& S2Y2 = new_kx2l_mat1(k, l); //Initialy S2Y2 = [Y1 S1] + S2Y2.copyBlockFromMatrix(0, 0, S1); + S2Y2.copyBlockFromMatrix(0, l, Y1); #ifdef HIOP_USE_MPI int ierr; - ierr = MPI_Allreduce(S2Y2.local_data(), _buff_2lxk, 2*l*k, MPI_DOUBLE, MPI_SUM, nlp->get_comm()); assert(ierr==MPI_SUCCESS); - ierr = MPI_Allreduce(W.local_data(), _buff_kxk, k*k, MPI_DOUBLE, MPI_SUM, nlp->get_comm()); assert(ierr==MPI_SUCCESS); - S2Y2.copyFrom(_buff_2lxk); - W.copyFrom(_buff_kxk); + ierr = MPI_Allreduce(S2Y2.local_data(), buff_2lxk_, 2 * l * k, MPI_DOUBLE, MPI_SUM, nlp_->get_comm()); + assert(ierr == MPI_SUCCESS); + ierr = MPI_Allreduce(W.local_data(), buff_kxk_, k * k, MPI_DOUBLE, MPI_SUM, nlp_->get_comm()); + assert(ierr == MPI_SUCCESS); + S2Y2.copyFrom(buff_2lxk_); + W.copyFrom(buff_kxk_); //also copy S1 and Y1 - S1.copyFromMatrixBlock(S2Y2, 0,0); - Y1.copyFromMatrixBlock(S2Y2, 0,l); + S1.copyFromMatrixBlock(S2Y2, 0, 0); + Y1.copyFromMatrixBlock(S2Y2, 0, l); #endif #ifdef HIOP_DEEPCHECKS - nlp->log->write("symMatTimesInverseTimesMatTrans: W first term is: ", W, hovMatrices); + nlp_->log->write("symMatTimesInverseTimesMatTrans: W first term is: ", W, hovMatrices); #endif //4. [S2] = V \ [S1^T] // [Y2] [Y1^T] @@ -568,159 +638,183 @@ symMatTimesInverseTimesMatTrans(double beta, hiopMatrixDense& W, //5. W = W-alpha*[S1 Y1]*[S2^T] // [Y2^T] S2Y2 = RHS_fortran; - alpha = 0-alpha; - hiopMatrixDense& S2=new_kxl_mat1(k,l); + alpha = 0 - alpha; + hiopMatrixDense& S2 = new_kxl_mat1(k,l); S2.copyFromMatrixBlock(S2Y2, 0, 0); S1.timesMatTrans_local(1.0, W, alpha, S2); - hiopMatrixDense& Y2=S2; + hiopMatrixDense& Y2 = S2; Y2.copyFromMatrixBlock(S2Y2, 0, l); Y1.timesMatTrans_local(1.0, W, alpha, Y2); - //nlp->log->write("symMatTimesInverseTimesMatTrans: Y1 is : ", Y1, hovMatrices); - //nlp->log->write("symMatTimesInverseTimesMatTrans: Y2 is : ", Y2, hovMatrices); - //nlp->log->write("symMatTimesInverseTimesMatTrans: W is : ", W, hovMatrices); + //nlp_->log->write("symMatTimesInverseTimesMatTrans: Y1 is : ", Y1, hovMatrices); + //nlp_->log->write("symMatTimesInverseTimesMatTrans: Y2 is : ", Y2, hovMatrices); + //nlp_->log->write("symMatTimesInverseTimesMatTrans: W is : ", W, hovMatrices); //we're done here #ifdef HIOP_DEEPCHECKS - nlp->log->write("symMatTimesInverseTimesMatTrans: final matrix is : ", W, hovMatrices); + nlp_->log->write("symMatTimesInverseTimesMatTrans: final matrix is : ", W, hovMatrices); #endif } - void hiopHessianLowRank::factorizeV() { - int N=V->n(), lda=N, info; - if(N==0) return; + int N = V_->n(); + int lda = N; + int info; + if(N == 0) { + return; + } #ifdef HIOP_DEEPCHECKS - nlp->log->write("factorizeV: V is ", *V, hovMatrices); + nlp_->log->write("factorizeV: V is ", *V_, hovMatrices); #endif - char uplo='L'; //V is upper in C++ so it's lower in fortran + char uplo = 'L'; //V is upper in C++ so it's lower in fortran - if(_V_ipiv_vec==NULL) _V_ipiv_vec=new int[N]; - else if(_V_ipiv_size!=N) { delete[] _V_ipiv_vec; _V_ipiv_vec=new int[N]; _V_ipiv_size=N; } + if(V_ipiv_vec_ == nullptr) { + V_ipiv_vec_ = new int[N]; + } else if(V_ipiv_size_ != N) { + delete[] V_ipiv_vec_; + V_ipiv_vec_ = new int[N]; + V_ipiv_size_ = N; + } - int lwork=-1;//inquire sizes + int lwork = -1;//inquire sizes double Vwork_tmp; - DSYTRF(&uplo, &N, V->local_data(), &lda, _V_ipiv_vec, &Vwork_tmp, &lwork, &info); - assert(info==0); + DSYTRF(&uplo, &N, V_->local_data(), &lda, V_ipiv_vec_, &Vwork_tmp, &lwork, &info); + assert(info == 0); - lwork=(int)Vwork_tmp; - if(lwork != _V_work_vec->get_size()) { - if(_V_work_vec!=NULL) { - delete _V_work_vec; + lwork = (int)Vwork_tmp; + if(lwork != V_work_vec_->get_size()) { + if(V_work_vec_ != nullptr) { + delete V_work_vec_; } - _V_work_vec = LinearAlgebraFactory::create_vector("DEFAULT", lwork); - } else assert(_V_work_vec); + V_work_vec_ = LinearAlgebraFactory::create_vector("DEFAULT", lwork); + } else { + assert(V_work_vec_); + } - DSYTRF(&uplo, &N, V->local_data(), &lda, _V_ipiv_vec, _V_work_vec->local_data(), &lwork, &info); + DSYTRF(&uplo, &N, V_->local_data(), &lda, V_ipiv_vec_, V_work_vec_->local_data(), &lwork, &info); - if(info<0) - nlp->log->printf(hovError, "hiopHessianLowRank::factorizeV error: %d argument to dsytrf has an illegal value\n", -info); - else if(info>0) - nlp->log->printf(hovError, "hiopHessianLowRank::factorizeV error: %d entry in the factorization's diagonal is exactly zero. Division by zero will occur if it a solve is attempted.\n", info); - assert(info==0); + if(info < 0) { + nlp_->log->printf(hovError, "hiopHessianLowRank::factorizeV error: %d argument to dsytrf has an illegal value\n", -info); + } else if(info > 0) { + nlp_->log->printf(hovError, "hiopHessianLowRank::factorizeV error: %d entry in the factorization's diagonal is exactly zero. Division by zero will occur if it a solve is attempted.\n", info); + } + assert(info == 0); #ifdef HIOP_DEEPCHECKS - nlp->log->write("factorizeV: factors of V: ", *V, hovMatrices); + nlp_->log->write("factorizeV: factors of V: ", *V_, hovMatrices); #endif } void hiopHessianLowRank::solveWithV(hiopVector& rhs_s, hiopVector& rhs_y) { - int N=V->n(); - if(N==0) return; - - int l=rhs_s.get_size(); + int N = V_->n(); + if(N == 0) { + return; + } + int l = rhs_s.get_size(); #ifdef HIOP_DEEPCHECKS - nlp->log->write("hiopHessianLowRank::solveWithV: RHS IN 's' part: ", rhs_s, hovMatrices); - nlp->log->write("hiopHessianLowRank::solveWithV: RHS IN 'y' part: ", rhs_y, hovMatrices); - hiopVector* rhs_saved= LinearAlgebraFactory::create_vector("DEFAULT", rhs_s.get_size()+rhs_y.get_size()); + nlp_->log->write("hiopHessianLowRank::solveWithV: RHS IN 's' part: ", rhs_s, hovMatrices); + nlp_->log->write("hiopHessianLowRank::solveWithV: RHS IN 'y' part: ", rhs_y, hovMatrices); + hiopVector* rhs_saved = LinearAlgebraFactory::create_vector("DEFAULT", rhs_s.get_size()+rhs_y.get_size()); rhs_saved->copyFromStarting(0, rhs_s); rhs_saved->copyFromStarting(l, rhs_y); #endif - int lda=N, one=1, info; - char uplo='L'; + int lda = N; + int one = 1; + int info; + char uplo ='L'; #ifdef HIOP_DEEPCHECKS - assert(N==rhs_s.get_size()+rhs_y.get_size()); + assert(N == rhs_s.get_size() + rhs_y.get_size()); #endif - hiopVector& rhs=new_2l_vec1(l); + hiopVector& rhs = new_2l_vec1(l); rhs.copyFromStarting(0, rhs_s); rhs.copyFromStarting(l, rhs_y); - DSYTRS(&uplo, &N, &one, V->local_data(), &lda, _V_ipiv_vec, rhs.local_data(), &N, &info); + DSYTRS(&uplo, &N, &one, V_->local_data(), &lda, V_ipiv_vec_, rhs.local_data(), &N, &info); - if(info<0) nlp->log->printf(hovError, "hiopHessianLowRank::solveWithV error: %d argument to dsytrf has an illegal value\n", -info); - assert(info==0); + if(info < 0) { + nlp_->log->printf(hovError, "hiopHessianLowRank::solveWithV error: %d argument to dsytrf has an illegal value\n", -info); + } + assert(info == 0); //copy back the solution - rhs.copyToStarting(0,rhs_s); - rhs.copyToStarting(l,rhs_y); + rhs.copyToStarting(0, rhs_s); + rhs.copyToStarting(l, rhs_y); #ifdef HIOP_DEEPCHECKS - nlp->log->write("solveWithV: SOL OUT 's' part: ", rhs_s, hovMatrices); - nlp->log->write("solveWithV: SOL OUT 'y' part: ", rhs_y, hovMatrices); + nlp_->log->write("solveWithV: SOL OUT 's' part: ", rhs_s, hovMatrices); + nlp_->log->write("solveWithV: SOL OUT 'y' part: ", rhs_y, hovMatrices); //residual calculation - double nrmrhs=rhs_saved->infnorm(); - _Vmat->timesVec(1.0, *rhs_saved, -1.0, rhs); - double nrmres=rhs_saved->infnorm(); - //nlp->log->printf(hovLinAlgScalars, "hiopHessianLowRank::solveWithV 1rhs: rel resid norm=%g\n", nrmres/(1+nrmrhs)); - nlp->log->printf(hovScalars, "hiopHessianLowRank::solveWithV 1rhs: rel resid norm=%g\n", nrmres/(1+nrmrhs)); - if(nrmres>1e-8) - nlp->log->printf(hovWarning, "hiopHessianLowRank::solveWithV large residual=%g\n", nrmres); + double nrmrhs = rhs_saved->infnorm(); + Vmat_->timesVec(1.0, *rhs_saved, -1.0, rhs); + double nrmres = rhs_saved->infnorm(); + //nlp_->log->printf(hovLinAlgScalars, "hiopHessianLowRank::solveWithV 1rhs: rel resid norm=%g\n", nrmres/(1+nrmrhs)); + nlp_->log->printf(hovScalars, "hiopHessianLowRank::solveWithV 1rhs: rel resid norm=%g\n", nrmres/(1+nrmrhs)); + if(nrmres > 1e-8) { + nlp_->log->printf(hovWarning, "hiopHessianLowRank::solveWithV large residual=%g\n", nrmres); + } delete rhs_saved; #endif - } void hiopHessianLowRank::solveWithV(hiopMatrixDense& rhs) { - int N=V->n(); - if(0==N) return; - + int N = V_->n(); + if(0 == N) { + return; + } #ifdef HIOP_DEEPCHECKS - nlp->log->write("solveWithV: RHS IN: ", rhs, hovMatrices); + nlp_->log->write("solveWithV: RHS IN: ", rhs, hovMatrices); hiopMatrixDense* rhs_saved = rhs.new_copy(); #endif //rhs is transpose in C++ - char uplo='L'; - int lda=N, ldb=N, nrhs=rhs.m(), info; + char uplo = 'L'; + int lda = N; + int ldb = N; + int nrhs = rhs.m(); + int info; #ifdef HIOP_DEEPCHECKS - assert(N==rhs.n()); + assert(N == rhs.n()); #endif - DSYTRS(&uplo, &N, &nrhs, V->local_data(), &lda, _V_ipiv_vec, rhs.local_data(), &ldb, &info); + DSYTRS(&uplo, &N, &nrhs, V_->local_data(), &lda, V_ipiv_vec_, rhs.local_data(), &ldb, &info); - if(info<0) nlp->log->printf(hovError, "hiopHessianLowRank::solveWithV error: %d argument to dsytrf has an illegal value\n", -info); - assert(info==0); + if(info < 0) { + nlp_->log->printf(hovError, "hiopHessianLowRank::solveWithV error: %d argument to dsytrf has an illegal value\n", -info); + } + assert(info == 0); #ifdef HIOP_DEEPCHECKS - nlp->log->write("solveWithV: SOL OUT: ", rhs, hovMatrices); + nlp_->log->write("solveWithV: SOL OUT: ", rhs, hovMatrices); hiopMatrixDense& sol = rhs; //matrix of solutions - /// TODO: get rid of these uses of specific hiopVector implementation hiopVector* x = LinearAlgebraFactory::create_vector("DEFAULT", rhs.n()); //again, keep in mind rhs is transposed hiopVector* r = LinearAlgebraFactory::create_vector("DEFAULT", rhs.n()); - double resnorm=0.0; - for(int k=0; kgetRow(k, *r); - sol.getRow(k,*x); - double nrmrhs = r->infnorm();//nrmrhs=.0; - _Vmat->timesVec(1.0, *r, -1.0, *x); + sol.getRow(k, *x); + double nrmrhs = r->infnorm(); //nrmrhs=.0; + Vmat_->timesVec(1.0, *r, -1.0, *x); double nrmres = r->infnorm(); - if(nrmres>1e-8) - nlp->log->printf(hovWarning, "hiopHessianLowRank::solveWithV mult-rhs: rhs number %d has large resid norm=%g\n", k, nrmres); - if(nrmres/(nrmrhs+1)>resnorm) resnorm=nrmres/(nrmrhs+1); + if(nrmres > 1e-8) { + nlp_->log->printf(hovWarning, "hiopHessianLowRank::solveWithV mult-rhs: rhs number %d has large resid norm=%g\n", k, nrmres); + } + if(nrmres / (nrmrhs + 1) > resnorm) { + resnorm = nrmres / (nrmrhs + 1); + } } - nlp->log->printf(hovLinAlgScalars, "hiopHessianLowRank::solveWithV mult-rhs: rel resid norm=%g\n", resnorm); + nlp_->log->printf(hovLinAlgScalars, "hiopHessianLowRank::solveWithV mult-rhs: rel resid norm=%g\n", resnorm); delete x; delete r; delete rhs_saved; @@ -730,188 +824,214 @@ void hiopHessianLowRank::solveWithV(hiopMatrixDense& rhs) void hiopHessianLowRank::growL(const int& lmem_curr, const int& lmem_max, const hiopVector& YTs) { - int l=L->m(); + int l = L_->m(); #ifdef HIOP_DEEPCHECKS - assert(l==L->n()); - assert(lmem_curr==l); - assert(lmem_max>=l); + assert(l == L_->n()); + assert(lmem_curr == l); + assert(lmem_max >= l); #endif //newL = [ L 0] // [ Y^T*s 0] hiopMatrixDense* newL = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l+1, l+1); assert(newL); //copy from L to newL - newL->copyBlockFromMatrix(0,0, *L); + newL->copyBlockFromMatrix(0, 0, *L_); - double* newL_mat=newL->local_data(); //doing the rest here - const double* YTs_vec=YTs.local_data_const(); + double* newL_mat = newL->local_data(); //doing the rest here + const double* YTs_vec = YTs.local_data_const(); //for(int j=0; jget_size(); - assert(l==lmem_curr); - assert(lmem_max>=l); + int l = D_->get_size(); + assert(l == lmem_curr); + assert(lmem_max >= l); - hiopVector* Dnew=LinearAlgebraFactory::create_vector("DEFAULT", l+1); - double* Dnew_vec=Dnew->local_data(); - memcpy(Dnew_vec, D->local_data_const(), l*sizeof(double)); - Dnew_vec[l]=sTy; + hiopVector* Dnew = LinearAlgebraFactory::create_vector("DEFAULT", l + 1); + double* Dnew_vec = Dnew->local_data(); + memcpy(Dnew_vec, D_->local_data_const(), l * sizeof(double)); + Dnew_vec[l] = sTy; - delete D; - D=Dnew; + delete D_; + D_ = Dnew; } -/* L_{ij} = s_{i-1}^T y_{j-1}, if i>j, otherwise zero. Here i,j = 0,1,...,l_curr-1 +/* L_{ij} = s_{i-1}^T y_{j-1}, if i>j, otherwise zero. Here i,j = 0,1,...,l_curr_-1 * L_new = lift and shift L to the left; replace last row with [Yts;0] */ void hiopHessianLowRank::updateL(const hiopVector& YTs, const double& sTy) { - int l=YTs.get_size(); - assert(l==L->m()); - assert(l==L->n()); + int l = YTs.get_size(); + assert(l == L_->m()); + assert(l == L_->n()); #ifdef HIOP_DEEPCHECKS - assert(l_curr==l); - assert(l_curr==l_max); + assert(l_curr_ == l); + assert(l_curr_ == l_max_); #endif - const int lm1=l-1; - double* L_mat=L->local_data(); - const double* yts_vec=YTs.local_data_const(); - for(int i=1; ilocal_data(); + const double* yts_vec = YTs.local_data_const(); + for(int i = 1; i < lm1; i++) { + for(int j = 0; j < i; j++) { //L_mat[i][j] = L_mat[i+1][j+1]; - L_mat[i*l+j] = L_mat[(i+1)*l+j+1]; + L_mat[i * l + j] = L_mat[(i + 1) * l + j + 1]; } } //is this really needed? - //for(int i=0; iget_size(); - double* D_vec = D->local_data(); - for(int i=0; iget_size(); + double* D_vec = D_->local_data(); + for(int i = 0; i < l - 1; i++) + D_vec[i] = D_vec[i + 1]; + D_vec[l - 1] = sTy; } - -hiopVector& hiopHessianLowRank::new_l_vec1(int l) +hiopVector& hiopHessianLowRank::new_l_vec1(int l) { - if(_l_vec1!=NULL && _l_vec1->get_size()==l) return *_l_vec1; - - if(_l_vec1!=NULL) { - delete _l_vec1; + if(l_vec1_ != nullptr && l_vec1_->get_size() == l) { + return *l_vec1_; + } + + if(l_vec1_ != nullptr) { + delete l_vec1_; } - _l_vec1= LinearAlgebraFactory::create_vector("DEFAULT", l); - return *_l_vec1; + l_vec1_ = LinearAlgebraFactory::create_vector("DEFAULT", l); + return *l_vec1_; } -hiopVector& hiopHessianLowRank::new_l_vec2(int l) +hiopVector& hiopHessianLowRank::new_l_vec2(int l) { - if(_l_vec2!=NULL && _l_vec2->get_size()==l) return *_l_vec2; - - if(_l_vec2!=NULL) { - delete _l_vec2; + if(l_vec2_ != nullptr && l_vec2_->get_size() == l) { + return *l_vec2_; + } + + if(l_vec2_ != nullptr) { + delete l_vec2_; } - _l_vec2= LinearAlgebraFactory::create_vector("DEFAULT", l); - return *_l_vec2; + l_vec2_ = LinearAlgebraFactory::create_vector("DEFAULT", l); + return *l_vec2_; } hiopMatrixDense& hiopHessianLowRank::new_lxl_mat1(int l) { - if(_lxl_mat1!=NULL) { - if( l==_lxl_mat1->m() ) { - return *_lxl_mat1; + if(lxl_mat1_ != nullptr) { + if(l == lxl_mat1_->m()) { + return *lxl_mat1_; } else { - delete _lxl_mat1; - _lxl_mat1=NULL; + delete lxl_mat1_; + lxl_mat1_ = nullptr; } } - _lxl_mat1 = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, l); - return *_lxl_mat1; + lxl_mat1_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, l); + return *lxl_mat1_; } + hiopMatrixDense& hiopHessianLowRank::new_kx2l_mat1(int k, int l) { - int twol=2*l; - if(NULL!=_kx2l_mat1) { - assert(_kx2l_mat1->m()==k); - if( twol==_kx2l_mat1->n() ) { - return *_kx2l_mat1; + int twol = 2 * l; + if(nullptr != kx2l_mat1_) { + assert(kx2l_mat1_->m() == k); + if( twol == kx2l_mat1_->n() ) { + return *kx2l_mat1_; } else { - delete _kx2l_mat1; - _kx2l_mat1=NULL; + delete kx2l_mat1_; + kx2l_mat1_ = nullptr; } } - _kx2l_mat1 = LinearAlgebraFactory::create_matrix_dense("DEFAULT", k, twol); + kx2l_mat1_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", k, twol); - return *_kx2l_mat1; + return *kx2l_mat1_; } hiopMatrixDense& hiopHessianLowRank::new_kxl_mat1(int k, int l) { - if(_kxl_mat1!=NULL) { - assert(_kxl_mat1->m()==k); - if( l==_kxl_mat1->n() ) { - return *_kxl_mat1; + if(kxl_mat1_ != nullptr) { + assert(kxl_mat1_->m() == k); + if( l == kxl_mat1_->n() ) { + return *kxl_mat1_; } else { - delete _kxl_mat1; - _kxl_mat1=NULL; + delete kxl_mat1_; + kxl_mat1_ = nullptr; } } - _kxl_mat1 = LinearAlgebraFactory::create_matrix_dense("DEFAULT", k, l); - return *_kxl_mat1; + kxl_mat1_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", k, l); + return *kxl_mat1_; } hiopMatrixDense& hiopHessianLowRank::new_S1(const hiopMatrixDense& X, const hiopMatrixDense& St) { //S1 is X*some_diag*S (kxl). Here St=S^T is lxn and X is kxn (l BFGS memory size, k number of constraints) - size_type k=X.m(), n=St.n(), l=St.m(); + size_type k = X.m(); + size_type n = St.n(); + size_type l = St.m(); #ifdef HIOP_DEEPCHECKS - assert(n==X.n()); - if(_S1!=NULL) - assert(_S1->m()==k); + assert(n == X.n()); + if(S1_ != nullptr) { + assert(S1_->m() == k); + } #endif - if(NULL!=_S1 && _S1->n()!=l) { delete _S1; _S1=NULL; } + if(nullptr != S1_ && S1_->n() != l) { + delete S1_; + S1_ = nullptr; + } - if(NULL==_S1) _S1=LinearAlgebraFactory::create_matrix_dense("DEFAULT", k, l); + if(nullptr == S1_) { + S1_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", k, l); + } - return *_S1; + return *S1_; } hiopMatrixDense& hiopHessianLowRank::new_Y1(const hiopMatrixDense& X, const hiopMatrixDense& Yt) { //Y1 is X*somediag*Y (kxl). Here Yt=Y^T is lxn, X is kxn - size_type k=X.m(), n=Yt.n(), l=Yt.m(); + size_type k = X.m(); + size_type n = Yt.n(); + size_type l = Yt.m(); #ifdef HIOP_DEEPCHECKS - assert(X.n()==n); - if(_Y1!=NULL) assert(_Y1->m()==k); + assert(X.n() == n); + if(Y1_ != nullptr) { + assert(Y1_->m() == k); + } #endif - if(NULL!=_Y1 && _Y1->n()!=l) { delete _Y1; _Y1=NULL; } + if(nullptr != Y1_ && Y1_->n() != l) { + delete Y1_; + Y1_ = nullptr; + } - if(NULL==_Y1) _Y1=LinearAlgebraFactory::create_matrix_dense("DEFAULT", k, l); - return *_Y1; + if(nullptr == Y1_) { + Y1_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", k, l); + } + return *Y1_; } #ifdef HIOP_DEEPCHECKS @@ -925,85 +1045,87 @@ void hiopHessianLowRank::timesVec_noLogBarrierTerm(double beta, hiopVector& y, d void hiopHessianLowRank::timesVecCmn(double beta, hiopVector& y, double alpha, const hiopVector& x, bool addLogTerm) const { - size_type n=St->n(); - assert(l_curr==St->m()); - assert(y.get_size()==n); - assert(St->get_local_size_n() == Yt->get_local_size_n()); + size_type n = St_->n(); + assert(l_curr_ == St_->m()); + assert(y.get_size() == n); + assert(St_->get_local_size_n() == Yt_->get_local_size_n()); //we have B+=B-B*s*B*s'/(s'*B*s)+yy'/(y'*s) - //B0 is sigma*I. There is an additional diagonal log-barrier term _Dx + //B0 is sigma*I. There is an additional diagonal log-barrier term Dx_ - bool print=false; + bool print = false; if(print) { - nlp->log->printf(hovMatrices, "---hiopHessianLowRank::timesVec \n"); - nlp->log->write("S=", *St, hovMatrices); - nlp->log->write("Y=", *Yt, hovMatrices); - nlp->log->write("DhInv=", *DhInv, hovMatrices); - nlp->log->printf(hovMatrices, "sigma=%22.16e; addLogTerm=%d;\n", sigma, addLogTerm); - if(addLogTerm) - nlp->log->write("Dx=", *_Dx, hovMatrices); - nlp->log->printf(hovMatrices, "y=beta*y + alpha*this*x : beta=%g alpha=%g\n", beta, alpha); - nlp->log->write("x_in=", x, hovMatrices); - nlp->log->write("y_in=", y, hovMatrices); + nlp_->log->printf(hovMatrices, "---hiopHessianLowRank::timesVec \n"); + nlp_->log->write("S=", *St_, hovMatrices); + nlp_->log->write("Y=", *Yt_, hovMatrices); + nlp_->log->write("DhInv=", *DhInv_, hovMatrices); + nlp_->log->printf(hovMatrices, "sigma=%22.16e; addLogTerm=%d;\n", sigma_, addLogTerm); + if(addLogTerm) { + nlp_->log->write("Dx=", *Dx_, hovMatrices); + } + nlp_->log->printf(hovMatrices, "y=beta*y + alpha*this*x : beta=%g alpha=%g\n", beta, alpha); + nlp_->log->write("x_in=", x, hovMatrices); + nlp_->log->write("y_in=", y, hovMatrices); } //allocate and compute a_k and b_k - a.resize(l_curr, nullptr); - b.resize(l_curr, nullptr); - int n_local = Yt->get_local_size_n(); - for(int k=0; kget_local_size_n(); + for(int k = 0; k < l_curr_; k++) { //bk=yk/sqrt(yk'*sk) - yk->copyFrom(Yt->local_data() + k*n_local); - sk->copyFrom(St->local_data() + k*n_local); - double skTyk=yk->dotProductWith(*sk); + yk_->copyFrom(Yt_->local_data() + k * n_local); + sk_->copyFrom(St_->local_data() + k * n_local); + double skTyk = yk_->dotProductWith(*sk_); if(skTyk < std::numeric_limits::epsilon()) { - nlp->log->printf(hovLinAlgScalars, "hiopHessianLowRank: ||s_k^T*y_k||=%12.6e too small..." + nlp_->log->printf(hovLinAlgScalars, "hiopHessianLowRank: ||s_k^T*y_k||=%12.6e too small..." " set it to machine precision = %12.6e \n", skTyk, std::numeric_limits::epsilon()); skTyk = std::numeric_limits::epsilon(); } - if(a[k] == nullptr && b[k] == nullptr) { - b[k] = nlp->alloc_primal_vec(); - a[k] = nlp->alloc_primal_vec(); + if(a_[k] == nullptr && b_[k] == nullptr) { + b_[k] = nlp_->alloc_primal_vec(); + a_[k] = nlp_->alloc_primal_vec(); } - b[k]->copyFrom(*yk); - b[k]->scale(1/sqrt(skTyk)); + b_[k]->copyFrom(*yk_); + b_[k]->scale(1/sqrt(skTyk)); //compute ak by an inner loop - a[k]->copyFrom(*sk); - a[k]->scale(sigma); - - for(int i=0; idotProductWith(*sk); - a[k]->axpy(+biTsk, *b[i]); - double aiTsk = a[i]->dotProductWith(*sk); - a[k]->axpy(-aiTsk, *a[i]); + a_[k]->copyFrom(*sk_); + a_[k]->scale(sigma_); + + for(int i = 0; i < k; i++) { + double biTsk = b_[i]->dotProductWith(*sk_); + a_[k]->axpy(+biTsk, *b_[i]); + double aiTsk = a_[i]->dotProductWith(*sk_); + a_[k]->axpy(-aiTsk, *a_[i]); } - double skTak = a[k]->dotProductWith(*sk); - a[k]->scale(1/sqrt(skTak)); + double skTak = a_[k]->dotProductWith(*sk_); + a_[k]->scale(1 / sqrt(skTak)); } - //now we have B= B_0 + sum{ bk bk' - ak ak' : k=0,1,...,l_curr-1} + //now we have B= B_0 + sum{ bk bk' - ak ak' : k=0,1,...,l_curr_-1} //compute the product with x - //y = beta*y+alpha*(B0+Dx)*x + alpha* sum { bk'x bk - ak'x ak : k=0,1,...,l_curr-1} + //y = beta*y+alpha*(B0+Dx)*x + alpha* sum { bk'x bk - ak'x ak : k=0,1,...,l_curr_-1} y.scale(beta); - if(addLogTerm) - y.axzpy(alpha,x,*_Dx); + if(addLogTerm) { + y.axzpy(alpha, x, *Dx_); + } - y.axpy(alpha*sigma, x); + y.axpy(alpha * sigma_, x); - for(int k=0; kdotProductWith(x); - double akTx = a[k]->dotProductWith(x); + for(int k = 0; k < l_curr_; k++) { + double bkTx = b_[k]->dotProductWith(x); + double akTx = a_[k]->dotProductWith(x); - y.axpy( alpha*bkTx, *b[k]); - y.axpy(-alpha*akTx, *a[k]); + y.axpy( alpha * bkTx, *b_[k]); + y.axpy(-alpha * akTx, *a_[k]); } if(print) { - nlp->log->write("y_out=", y, hovMatrices); + nlp_->log->write("y_out=", y, hovMatrices); } } @@ -1027,79 +1149,86 @@ void hiopHessianLowRank::timesVec(double beta, hiopVector& y, double alpha, cons * The ops are perform locally. The reduce is done separately/externally to decrease comm */ void hiopHessianLowRank:: -symmMatTimesDiagTimesMatTrans_local(double beta, hiopMatrixDense& W, - double alpha, const hiopMatrixDense& X, - const hiopVector& d) +symmMatTimesDiagTimesMatTrans_local(double beta, + hiopMatrixDense& W, + double alpha, + const hiopMatrixDense& X, + const hiopVector& d) { - size_type k=W.m(); - size_type n=X.n(); - size_t n_local=X.get_local_size_n(); + size_type k = W.m(); + size_type n = X.n(); + size_t n_local = X.get_local_size_n(); - assert(X.m()==k); + assert(X.m() == k); #ifdef HIOP_DEEPCHECKS - assert(W.n()==k); - assert(d.get_size()==n); - assert(d.get_local_size()==n_local); + assert(W.n() == k); + assert(d.get_size() == n); + assert(d.get_local_size() == n_local); #endif //#define chunk 512; //!opt const double *xi, *xj; double acc; - double *Wdata=W.local_data(); - const double *Xdata=X.local_data_const(); - const double* dd=d.local_data_const(); - for(int i=0; i(nlp_); - // assert(nlp==NULL && "only NLPs with a small number of constraints are supported by HessianLowRank"); + const hiopNlpDenseConstraints* nlp = dynamic_cast(nlp_in); + // assert(nlp= = nullptr && "only NLPs with a small number of constraints are supported by HessianLowRank"); - H0 = nlp->alloc_primal_vec(); - St = nlp->alloc_multivector_primal(0,l_max); - Yt = St->alloc_clone(); //faster than nlp->alloc_multivector_primal(...); + H0_ = nlp->alloc_primal_vec(); + St_ = nlp->alloc_multivector_primal(0, l_max_); + Yt_ = St_->alloc_clone(); //faster than nlp->alloc_multivector_primal(...); //these are local - R = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); - D = LinearAlgebraFactory::create_vector("DEFAULT", 0); - - //the previous iteration's objects are set to NULL - _it_prev=NULL; _grad_f_prev=NULL; _Jac_c_prev=NULL; _Jac_d_prev=NULL; + R_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); + D_ = LinearAlgebraFactory::create_vector("DEFAULT", 0); + //the previous iteration's objects are set to nullptr + it_prev_ = nullptr; + grad_f_prev_ = nullptr; + Jac_c_prev_ = nullptr; + Jac_d_prev_ = nullptr; //internal buffers for memory pool (none of them should be in n) #ifdef HIOP_USE_MPI - _buff_kxk = new double[nlp->m() * nlp->m()]; - _buff_lxk = new double[nlp->m() * l_max]; - _buff_lxl = new double[l_max*l_max]; + buff_kxk_ = new double[nlp->m() * nlp->m()]; + buff_lxk_ = new double[nlp->m() * l_max_]; + buff_lxl_ = new double[l_max_ * l_max_]; #else //not needed in non-MPI mode - _buff_kxk = NULL; - _buff_lxk = NULL; - _buff_lxl = NULL; + buff_kxk_ = nullptr; + buff_lxk_ = nullptr; + buff_lxl_ = nullptr; #endif //auxiliary objects/buffers - _S1=_Y1=_S3=NULL; - _DpYtH0Y=NULL; - _l_vec1 = _l_vec2 = NULL; - _n_vec1 = H0->alloc_clone(); - _n_vec2 = H0->alloc_clone(); - //H0->setToConstant(sigma); - - sigma=sigma0; - sigma_update_strategy = SIGMA_STRATEGY1; - sigma_safe_min=1e-8; - sigma_safe_max=1e+8; - nlp->log->printf(hovScalars, "Hessian Low Rank: initial sigma is %g\n", sigma); + S1_ = nullptr; + Y1_ = nullptr; + S3_ = nullptr; + DpYtH0Y_ = nullptr; + l_vec1_ = nullptr; + l_vec2_ = nullptr; + n_vec1_ = H0_->alloc_clone(); + n_vec2_ = H0_->alloc_clone(); + //H0_->setToConstant(sigma_); + + sigma_ = sigma0_; + sigma_update_strategy_ = SIGMA_STRATEGY1; + sigma_safe_min_ = 1e-8; + sigma_safe_max_ = 1e+8; + nlp->log->printf(hovScalars, "Hessian Low Rank: initial sigma is %g\n", sigma_); } + hiopHessianInvLowRank_obsolette::~hiopHessianInvLowRank_obsolette() { - if(H0) delete H0; - if(St) delete St; - if(Yt) delete Yt; - if(R) delete R; - if(D) delete D; - - if(_it_prev) delete _it_prev; - if(_grad_f_prev)delete _grad_f_prev; - if(_Jac_c_prev) delete _Jac_c_prev; - if(_Jac_d_prev) delete _Jac_d_prev; - - if(_buff_kxk) delete[] _buff_kxk; - if(_buff_lxk) delete[] _buff_lxk; - if(_buff_lxl) delete[] _buff_lxl; - if(_S1) delete _S1; - if(_Y1) delete _Y1; - if(_DpYtH0Y) delete _DpYtH0Y; - if(_S3) delete _S3; - if(_l_vec1) delete _l_vec1; - if(_l_vec2) delete _l_vec2; - if(_n_vec1) delete _n_vec1; - if(_n_vec2) delete _n_vec2; + delete H0_; + delete St_; + delete Yt_; + delete R_; + delete D_; + + delete it_prev_; + delete grad_f_prev_; + delete Jac_c_prev_; + delete Jac_d_prev_; + + if(buff_kxk_) { + delete[] buff_kxk_; + } + if(buff_lxk_) { + delete[] buff_lxk_; + } + if(buff_lxl_) { + delete[] buff_lxl_; + } + delete S1_; + delete Y1_; + delete DpYtH0Y_; + delete S3_; + delete l_vec1_; + delete l_vec2_; + delete n_vec1_; + delete n_vec2_; } #include -bool hiopHessianInvLowRank_obsolette:: -update(const hiopIterate& it_curr, const hiopVector& grad_f_curr_, - const hiopMatrix& Jac_c_curr_, const hiopMatrix& Jac_d_curr_) +bool hiopHessianInvLowRank_obsolette::update(const hiopIterate& it_curr, + const hiopVector& grad_f_curr, + const hiopMatrix& Jac_c_curr_in, + const hiopMatrix& Jac_d_curr_in) { - const hiopMatrixDense& Jac_c_curr = dynamic_cast(Jac_c_curr_); - const hiopMatrixDense& Jac_d_curr = dynamic_cast(Jac_d_curr_); + const hiopMatrixDense& Jac_c_curr = dynamic_cast(Jac_c_curr_in); + const hiopMatrixDense& Jac_d_curr = dynamic_cast(Jac_d_curr_in); #ifdef HIOP_DEEPCHECKS - assert(it_curr.zl->matchesPattern(nlp->get_ixl())); - assert(it_curr.zu->matchesPattern(nlp->get_ixu())); - assert(it_curr.sxl->matchesPattern(nlp->get_ixl())); - assert(it_curr.sxu->matchesPattern(nlp->get_ixu())); + assert(it_curr.zl->matchesPattern(nlp_->get_ixl())); + assert(it_curr.zu->matchesPattern(nlp_->get_ixu())); + assert(it_curr.sxl->matchesPattern(nlp_->get_ixl())); + assert(it_curr.sxu->matchesPattern(nlp_->get_ixu())); #endif - if(l_curr>0) { - size_type n=grad_f_curr_.get_size(); + if(l_curr_ > 0) { + size_type n = grad_f_curr.get_size(); //compute s_new = x_curr-x_prev - hiopVector& s_new = new_n_vec1(n); s_new.copyFrom(*it_curr.x); s_new.axpy(-1.,*_it_prev->x); - double s_infnorm=s_new.infnorm(); - if(s_infnorm>=100*std::numeric_limits::epsilon()) { //norm of s not too small + hiopVector& s_new = new_n_vec1(n); + s_new.copyFrom(*it_curr.x); + s_new.axpy(-1., *it_prev_->x); + double s_infnorm = s_new.infnorm(); + if(s_infnorm >= 100 * std::numeric_limits::epsilon()) { //norm of s not too small //compute y_new = \grad J(x_curr,\lambda_curr) - \grad J(x_prev, \lambda_curr) (yes, J(x_prev, \lambda_curr)) // = graf_f_curr-grad_f_prev + (Jac_c_curr-Jac_c_prev)yc_curr+ (Jac_d_curr-Jac_c_prev)yd_curr - zl_curr*s_new + zu_curr*s_new hiopVector& y_new = new_n_vec2(n); - y_new.copyFrom(grad_f_curr_); - y_new.axpy(-1., *_grad_f_prev); + y_new.copyFrom(grad_f_curr); + y_new.axpy(-1., *grad_f_prev_); Jac_c_curr.transTimesVec (1.0, y_new, 1.0, *it_curr.yc); - _Jac_c_prev->transTimesVec(1.0, y_new,-1.0, *it_curr.yc); //!opt if nlp->Jac_c_isLinear no need for the multiplications + Jac_c_prev_->transTimesVec(1.0, y_new,-1.0, *it_curr.yc); //!opt if nlp_->Jac_c_isLinear no need for the multiplications Jac_d_curr.transTimesVec (1.0, y_new, 1.0, *it_curr.yd); //!opt same here - _Jac_d_prev->transTimesVec(1.0, y_new,-1.0, *it_curr.yd); + Jac_d_prev_->transTimesVec(1.0, y_new,-1.0, *it_curr.yd); //y_new.axzpy(-1.0, s_new, *it_curr.zl); //y_new.axzpy( 1.0, s_new, *it_curr.zu); - double sTy = s_new.dotProductWith(y_new), s_nrm2=s_new.twonorm(), y_nrm2=y_new.twonorm(); - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: s^T*y=%20.14e ||s||=%20.14e ||y||=%20.14e\n", sTy, s_nrm2, y_nrm2); - nlp->log->write("hiopHessianInvLowRank_obsolette s_new",s_new, hovIteration); - nlp->log->write("hiopHessianInvLowRank_obsolette y_new",y_new, hovIteration); - - if(sTy>s_nrm2*y_nrm2*std::numeric_limits::epsilon()) { //sTy far away from zero - //compute the new column in R, update S and Y (either augment them or shift cols and add s_new and y_new) - hiopVector& STy = new_l_vec1(l_curr-1); - St->timesVec(0.0, STy, 1.0, y_new); - //update representation - if(l_currappendRow(s_new); - Yt->appendRow(y_new); - growR(l_curr, l_max, STy, sTy); - growD(l_curr, l_max, sTy); - l_curr++; - } else { - //shift - St->shiftRows(-1); - Yt->shiftRows(-1); - St->replaceRow(l_max-1, s_new); - Yt->replaceRow(l_max-1, y_new); - updateR(STy,sTy); - updateD(sTy); - l_curr=l_max; - } - - //update B0 (i.e., sigma) - switch (sigma_update_strategy ) { - case SIGMA_STRATEGY1: - sigma=sTy/(s_nrm2*s_nrm2); - break; - case SIGMA_STRATEGY2: - sigma=y_nrm2*y_nrm2/sTy; - break; - case SIGMA_STRATEGY3: - sigma=sqrt(s_nrm2*s_nrm2 / y_nrm2 / y_nrm2); - break; - case SIGMA_STRATEGY4: - sigma=0.5*(sTy/(s_nrm2*s_nrm2)+y_nrm2*y_nrm2/sTy); - break; - case SIGMA_CONSTANT: - sigma=sigma0; - break; - default: - assert(false && "Option value for sigma_update_strategy was not recognized."); - break; - } // else of the switch - //safe guard it - sigma=fmax(fmin(sigma_safe_max, sigma), sigma_safe_min); - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: sigma was updated to %16.10e\n", sigma); + double sTy = s_new.dotProductWith(y_new); + double s_nrm2 = s_new.twonorm(); + double y_nrm2 = y_new.twonorm(); + nlp_->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: s^T*y=%20.14e ||s||=%20.14e ||y||=%20.14e\n", sTy, s_nrm2, y_nrm2); + nlp_->log->write("hiopHessianInvLowRank_obsolette s_new",s_new, hovIteration); + nlp_->log->write("hiopHessianInvLowRank_obsolette y_new",y_new, hovIteration); + + if(sTy > s_nrm2 * y_nrm2 * std::numeric_limits::epsilon()) { //sTy far away from zero + //compute the new column in R, update S and Y (either augment them or shift cols and add s_new and y_new) + hiopVector& STy = new_l_vec1(l_curr_ - 1); + St_->timesVec(0.0, STy, 1.0, y_new); + //update representation + if(l_curr_ < l_max_) { + //just grow/augment the matrices + St_->appendRow(s_new); + Yt_->appendRow(y_new); + growR(l_curr_, l_max_, STy, sTy); + growD(l_curr_, l_max_, sTy); + l_curr_++; + } else { + //shift + St_->shiftRows(-1); + Yt_->shiftRows(-1); + St_->replaceRow(l_max_ - 1, s_new); + Yt_->replaceRow(l_max_ - 1, y_new); + updateR(STy, sTy); + updateD(sTy); + l_curr_ = l_max_; + } + + //update B0 (i.e., sigma) + switch (sigma_update_strategy_ ) { + case SIGMA_STRATEGY1: + sigma_ = sTy / (s_nrm2 * s_nrm2); + break; + case SIGMA_STRATEGY2: + sigma_ = y_nrm2 * y_nrm2 / sTy; + break; + case SIGMA_STRATEGY3: + sigma_ = sqrt(s_nrm2 * s_nrm2 / y_nrm2 / y_nrm2); + break; + case SIGMA_STRATEGY4: + sigma_ = 0.5 * (sTy / (s_nrm2 * s_nrm2) + y_nrm2 * y_nrm2 / sTy); + break; + case SIGMA_CONSTANT: + sigma_ = sigma0_; + break; + default: + assert(false && "Option value for sigma_update_strategy_ was not recognized."); + break; + } // else of the switch + + //safe guard it + sigma_ = fmax(fmin(sigma_safe_max_, sigma_), sigma_safe_min_); + nlp_->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: sigma was updated to %16.10e\n", sigma_); } else { //sTy is too small or negative -> skip - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: s^T*y=%12.6e not positive enough... skipping the Hessian update\n", sTy); + nlp_->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: s^T*y=%12.6e not positive enough... skipping the Hessian update\n", sTy); } } else {// norm of s_new is too small -> skip - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: ||s_new||=%12.6e too small... skipping the Hessian update\n", s_infnorm); + nlp_->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: ||s_new||=%12.6e too small... skipping the Hessian update\n", s_infnorm); } //save this stuff for next update - _it_prev->copyFrom(it_curr); _grad_f_prev->copyFrom(grad_f_curr_); - _Jac_c_prev->copyFrom(Jac_c_curr); _Jac_d_prev->copyFrom(Jac_d_curr); - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: storing the iteration info as 'previous'\n", s_infnorm); + it_prev_->copyFrom(it_curr); + grad_f_prev_->copyFrom(grad_f_curr); + Jac_c_prev_->copyFrom(Jac_c_curr); + Jac_d_prev_->copyFrom(Jac_d_curr); + nlp_->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: storing the iteration info as 'previous'\n", s_infnorm); } else { //this is the first optimization iterate, just save the iterate and exit - if(NULL==_it_prev) _it_prev = it_curr.new_copy(); - if(NULL==_grad_f_prev) _grad_f_prev = grad_f_curr_.new_copy(); - if(NULL==_Jac_c_prev) _Jac_c_prev = Jac_c_curr.new_copy(); - if(NULL==_Jac_d_prev) _Jac_d_prev = Jac_d_curr.new_copy(); - - nlp->log->printf(hovLinAlgScalarsVerb, "HessianInvLowRank on first update, just saving iteration\n"); + if(nullptr == it_prev_) { + it_prev_ = it_curr.new_copy(); + } + if(nullptr == grad_f_prev_) { + grad_f_prev_ = grad_f_curr.new_copy(); + }if(nullptr == Jac_c_prev_) { + Jac_c_prev_ = Jac_c_curr.new_copy(); + } + if(nullptr == Jac_d_prev_) { + Jac_d_prev_ = Jac_d_curr.new_copy(); + } + nlp_->log->printf(hovLinAlgScalarsVerb, "HessianInvLowRank on first update, just saving iteration\n"); - l_curr++; + l_curr_++; } return true; } bool hiopHessianInvLowRank_obsolette::updateLogBarrierDiagonal(const hiopVector& Dx) { - H0->setToConstant(sigma); - H0->axpy(1.0,Dx); + H0_->setToConstant(sigma_); + H0_->axpy(1.0, Dx); #ifdef HIOP_DEEPCHECKS - assert(H0->allPositive()); + assert(H0_->allPositive()); #endif - H0->invert(); - nlp->log->write("H0 InvHes", *H0, hovMatrices); + H0_->invert(); + nlp_->log->write("H0 InvHes", *H0_, hovMatrices); return true; } - /* Y = beta*Y + alpha*this*X * where 'this' is * M^{-1} = H0 + [S HO*Y] [ R^{-T}*(D+Y^T*H0*Y)*R^{-1} -R^{-T} ] [ S^T ] @@ -1312,45 +1467,49 @@ bool hiopHessianInvLowRank_obsolette::updateLogBarrierDiagonal(const hiopVector& */ void hiopHessianInvLowRank_obsolette::apply(double beta, hiopVector& y, double alpha, const hiopVector& x) { - size_type n=St->n(), l=St->m(); + size_type n = St_->n(); + size_type l = St_->m(); #ifdef HIOP_DEEPCHECKS - assert(y.get_size()==n); - assert(x.get_size()==n); - assert(H0->get_size()==n); + assert(y.get_size() == n); + assert(x.get_size() == n); + assert(H0_->get_size() == n); #endif //0. y = beta*y + alpha*H0*y y.scale(beta); - y.axzpy(alpha,*H0,x); + y.axzpy(alpha, *H0_,x); //1. stx = S^T*x and ytx = Y^T*H0*x - hiopVector &stx=new_l_vec1(l), &ytx=new_l_vec2(l), &H0x=new_n_vec1(n); - St->timesVec(0.0,stx,1.0,x); - H0x.copyFrom(x); H0x.componentMult(*H0); - Yt->timesVec(0.0,ytx,1.0,H0x); + hiopVector& stx = new_l_vec1(l); + hiopVector& ytx = new_l_vec2(l); + hiopVector& H0x = new_n_vec1(n); + St_->timesVec(0.0, stx, 1.0, x); + H0x.copyFrom(x); + H0x.componentMult(*H0_); + Yt_->timesVec(0.0, ytx, 1.0, H0x); //2.ytx = R^{-T}* [(D+Y^T*H0*Y)*R^{-1} stx - ytx ] // stx = -R^{-1}*stx //2.1. stx = R^{-1}*stx - triangularSolve(*R, stx); + triangularSolve(*R_, stx); //2.2 ytx = -ytx + (D+Y^T*H0*Y)*R^{-1} stx - hiopMatrixDense& DpYtH0Y = *_DpYtH0Y; // ! this is computed in symmetricTimesMat - DpYtH0Y.timesVec(-1.0,ytx, 1.0, stx); + hiopMatrixDense& DpYtH0Y = *DpYtH0Y_; // ! this is computed in symmetricTimesMat + DpYtH0Y.timesVec(-1.0, ytx, 1.0, stx); //2.3 ytx = R^{-T}* [(D+Y^T*H0*Y)*R^{-1} stx - ytx ] and stx = -R^{-1}*stx - triangularSolveTrans(*R,ytx); + triangularSolveTrans(*R_, ytx); //3. add alpha(S*ytx + H0*Y*stx) to y - St->transTimesVec(1.0,y,alpha,ytx); - hiopVector& H0Ystx=new_n_vec1(n); + St_->transTimesVec(1.0, y, alpha, ytx); + hiopVector& H0Ystx = new_n_vec1(n); //H0Ystx=H0*Y*stx - Yt->transTimesVec(0.0, H0Ystx, -1.0, stx); //-1.0 since we haven't negated stx in 2.3 - H0Ystx.componentMult(*H0); - y.axpy(alpha,H0Ystx); + Yt_->transTimesVec(0.0, H0Ystx, -1.0, stx); //-1.0 since we haven't negated stx in 2.3 + H0Ystx.componentMult(*H0_); + y.axpy(alpha, H0Ystx); } + void hiopHessianInvLowRank_obsolette::apply(double beta, hiopMatrix& Y, double alpha, const hiopMatrix& X) { assert(false && "not yet implemented"); } - /* W = beta*W + alpha*X*this*X^T * where 'this' is * M^{-1} = H0 + [S HO*Y] [ R^{-T}*(D+Y^T*H0*Y)*R^{-1} -R^{-T} ] [ S^T ] @@ -1359,53 +1518,64 @@ void hiopHessianInvLowRank_obsolette::apply(double beta, hiopMatrix& Y, double a * W is kxk, H0 is nxn, S,Y are nxl, R is upper triangular lxl, and X is kxn * Remember we store Yt=Y^T and St=S^T */ -void hiopHessianInvLowRank_obsolette:: -symmetricTimesMat(double beta, hiopMatrixDense& W, - double alpha, const hiopMatrixDense& X) +void hiopHessianInvLowRank_obsolette::symmetricTimesMat(double beta, + hiopMatrixDense& W, + double alpha, + const hiopMatrixDense& X) { - size_type n=St->n(), l=St->m(), k=W.m(); - assert(n==H0->get_size()); - assert(k==W.n()); - assert(l==Yt->m()); - assert(n==Yt->n()); assert(n==St->n()); - assert(k==X.m()); assert(n==X.n()); + size_type n = St_->n(); + size_type l = St_->m(); + size_type k = W.m(); + assert(n == H0_->get_size()); + assert(k == W.n()); + assert(l == Yt_->m()); + assert(n == Yt_->n()); + assert(n == St_->n()); + assert(k == X.m()); + assert(n == X.n()); //1.--compute W = beta*W + alpha*X*HO*X^T by calling symmMatTransTimesDiagTimesMat_local #ifdef HIOP_USE_MPI int myrank, ierr; - ierr=MPI_Comm_rank(nlp->get_comm(),&myrank); assert(MPI_SUCCESS==ierr); - if(0==myrank) - symmMatTimesDiagTimesMatTrans_local(beta,W,alpha,X,*H0); - else - symmMatTimesDiagTimesMatTrans_local(0.0,W,alpha,X,*H0); + ierr = MPI_Comm_rank(nlp_->get_comm(), &myrank); + assert(MPI_SUCCESS == ierr); + if(0 == myrank) { + symmMatTimesDiagTimesMatTrans_local(beta, W, alpha, X, *H0_); + } else { + symmMatTimesDiagTimesMatTrans_local(0.0, W, alpha, X, *H0_); + } //W will be MPI_All_reduced later #else - symmMatTimesDiagTimesMatTrans_local(beta,W,alpha,X,*H0); + symmMatTimesDiagTimesMatTrans_local(beta, W, alpha, X, *H0_); #endif //2.--compute S1=S^T*X^T=St*X^T and Y1=Y^T*H0*X^T=Yt*H0*X^T - hiopMatrixDense& S1 = new_S1(*St,X); - St->timesMatTrans_local(0.0,S1,1.0,X); - hiopMatrixDense& Y1 = new_Y1(*Yt,X); - matTimesDiagTimesMatTrans_local(Y1,*Yt,*H0,X); + hiopMatrixDense& S1 = new_S1(*St_,X); + St_->timesMatTrans_local(0.0, S1, 1.0, X); + hiopMatrixDense& Y1 = new_Y1(*Yt_, X); + matTimesDiagTimesMatTrans_local(Y1, *Yt_, *H0_, X); //3.--compute Y^T*H0*Y from D+Y^T*H0*Y - hiopMatrixDense& DpYtH0Y = new_DpYtH0Y(*Yt); - symmMatTimesDiagTimesMatTrans_local(0.0,DpYtH0Y, 1.0,*Yt,*H0); + hiopMatrixDense& DpYtH0Y = new_DpYtH0Y(*Yt_); + symmMatTimesDiagTimesMatTrans_local(0.0, DpYtH0Y, 1.0, *Yt_, *H0_); #ifdef HIOP_USE_MPI //!opt - use one buffer and one reduce call - ierr=MPI_Allreduce(S1.local_data(), _buff_lxk,l*k, MPI_DOUBLE,MPI_SUM,nlp->get_comm()); assert(ierr==MPI_SUCCESS); - S1.copyFrom(_buff_lxk); - ierr=MPI_Allreduce(Y1.local_data(), _buff_lxk,l*k, MPI_DOUBLE,MPI_SUM,nlp->get_comm()); assert(ierr==MPI_SUCCESS); - Y1.copyFrom(_buff_lxk); - ierr=MPI_Allreduce(W.local_data(), _buff_kxk,k*k, MPI_DOUBLE,MPI_SUM,nlp->get_comm()); assert(ierr==MPI_SUCCESS); - W.copyFrom(_buff_kxk); - - ierr=MPI_Allreduce(DpYtH0Y.local_data(), _buff_lxl,l*l, MPI_DOUBLE,MPI_SUM,nlp->get_comm()); assert(ierr==MPI_SUCCESS); - DpYtH0Y.copyFrom(_buff_lxl); + ierr = MPI_Allreduce(S1.local_data(), buff_lxk_, l * k, MPI_DOUBLE,MPI_SUM,nlp_->get_comm()); + assert(ierr == MPI_SUCCESS); + S1.copyFrom(buff_lxk_); + ierr = MPI_Allreduce(Y1.local_data(), buff_lxk_, l * k, MPI_DOUBLE,MPI_SUM,nlp_->get_comm()); + assert(ierr == MPI_SUCCESS); + Y1.copyFrom(buff_lxk_); + ierr = MPI_Allreduce(W.local_data(), buff_kxk_, k * k, MPI_DOUBLE,MPI_SUM,nlp_->get_comm()); + assert(ierr == MPI_SUCCESS); + W.copyFrom(buff_kxk_); + + ierr = MPI_Allreduce(DpYtH0Y.local_data(), buff_lxl_, l * l, MPI_DOUBLE,MPI_SUM,nlp_->get_comm()); + assert(ierr == MPI_SUCCESS); + DpYtH0Y.copyFrom(buff_lxl_); #endif //add D to finish calculating D+Y^T*H0*Y - DpYtH0Y.addDiagonal(1., *D); + DpYtH0Y.addDiagonal(1., *D_); //now we have W = beta*W+alpha*X*HO*X^T and the remaining term in the form // [S1^T Y1^T] [ R^{-T}*DpYtH0Y*R^{-1} -R^{-T} ] [ S1 ] @@ -1415,20 +1585,20 @@ symmetricTimesMat(double beta, hiopMatrixDense& W, // or W += S2^T *DpYtH0Y* S2 - S2^T *Y1 - Y1^T*S2 //4.-- compute S1 = R\S1 - triangularSolve(*R,S1); + triangularSolve(*R_, S1); //5.-- W += S2^T *DpYtH0Y* S2 //S3=DpYtH0Y*S2 - hiopMatrix& S3=new_S3(DpYtH0Y, S1); - DpYtH0Y.timesMat(0.0,S3,1.0,S1); + hiopMatrix& S3 = new_S3(DpYtH0Y, S1); + DpYtH0Y.timesMat(0.0, S3, 1.0, S1); // W += S2^T * S3 - S1.transTimesMat(1.0,W,1.0,S3); + S1.transTimesMat(1.0, W, 1.0, S3); //6.-- W += - S2^T *Y1 - (S2^T *Y1)^T // W = W - S2^T*Y1 - S1.transTimesMat(1.0,W,-1.0,Y1); + S1.transTimesMat(1.0, W, -1.0, Y1); // W = W - Y1*S2^T - Y1.transTimesMat(1.0,W,-1.0,S1); + Y1.transTimesMat(1.0, W, -1.0, S1); // -- done -- #ifdef HIOP_DEEPCHECKS @@ -1443,52 +1613,58 @@ symmetricTimesMat(double beta, hiopMatrixDense& W, * The ops are perform locally. The reduce is done separately/externally to decrease comm */ void hiopHessianInvLowRank_obsolette:: -symmMatTimesDiagTimesMatTrans_local(double beta, hiopMatrixDense& W, - double alpha, const hiopMatrixDense& X, - const hiopVector& d) +symmMatTimesDiagTimesMatTrans_local(double beta, + hiopMatrixDense& W, + double alpha, + const hiopMatrixDense& X, + const hiopVector& d) { - size_type k=W.m(); - size_type n=X.n(); - size_t n_local=X.get_local_size_n(); + size_type k = W.m(); + size_type n = X.n(); + size_t n_local = X.get_local_size_n(); - assert(X.m()==k); + assert(X.m() == k); #ifdef HIOP_DEEPCHECKS - assert(W.n()==k); - assert(d.get_size()==n); - assert(d.get_local_size()==n_local); + assert(W.n() == k); + assert(d.get_size() == n); + assert(d.get_local_size() == n_local); #endif //#define chunk 512; //!opt - const double *xi, *xj; + const double* xi; + const double* xj; double acc; - double* Wdata=W.local_data(); - const double* Xdata=X.local_data_const(); - const double* dd=d.local_data_const(); - for(int i=0; im(); + int l = R_->m(); #ifdef HIOP_DEEPCHECKS - assert(l==R->n()); - assert(lmem_curr==l); - assert(lmem_max>l); + assert(l == R_->n()); + assert(lmem_curr == l); + assert(lmem_max > l); #endif //newR = [ R S^T*y ] // [ 0 s^T*y ] - hiopMatrixDense* newR = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l+1, l+1); + hiopMatrixDense* newR = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l + 1, l + 1); assert(newR); //copy from R to newR - newR->copyBlockFromMatrix(0,0, *R); - - double* newR_mat=newR->local_data(); //doing the rest here - const double* STy_vec=STy.local_data_const(); - //for(int j=0; jcopyBlockFromMatrix(0, 0, *R_); + + double* newR_mat = newR->local_data(); //doing the rest here + const double* STy_vec = STy.local_data_const(); + //for(int j = 0; j < l; j++) { + // newR_mat[j][l] = STy_vec[j]; + //} + for(int j = 0; j < l; j++) { + newR_mat[j * (l + 1) + l] = STy_vec[j]; + } //newR_mat[l][l] = sTy; - newR_mat[l*(l+1)+l] = sTy; + newR_mat[l * (l + 1) + l] = sTy; //and the zero entries on the last row //for(int i=0; iget_size(); - assert(l==lmem_curr); - assert(lmem_max>l); + int l = D_->get_size(); + assert(l == lmem_curr); + assert(lmem_max > l); - hiopVector* Dnew=LinearAlgebraFactory::create_vector("DEFAULT", l+1); - double* Dnew_vec=Dnew->local_data(); - memcpy(Dnew_vec, D->local_data_const(), l*sizeof(double)); - Dnew_vec[l]=sTy; + hiopVector* Dnew = LinearAlgebraFactory::create_vector("DEFAULT", l + 1); + double* Dnew_vec = Dnew->local_data(); + memcpy(Dnew_vec, D_->local_data_const(), l * sizeof(double)); + Dnew_vec[l] = sTy; - delete D; - D=Dnew; + delete D_; + D_ = Dnew; } void hiopHessianInvLowRank_obsolette::updateR(const hiopVector& STy, const double& sTy) { - int l=STy.get_size(); - assert(l==R->m()); + int l = STy.get_size(); + assert(l == R_->m()); #ifdef HIOP_DEEPCHECKS - assert(l==R->n()); + assert(l == R_->n()); #endif - const int lm1=l-1; - double* R_mat=R->local_data(); - const double* sTy_vec=STy.local_data_const(); - for(int i=0; ilocal_data(); + const double* sTy_vec = STy.local_data_const(); + for(int i = 0; i < lm1; i++) { + for(int j = i; j < lm1; j++) { + //R_mat[i][j] = R_mat[i + 1][j + 1]; + R_mat[i * l + j] = R_mat[(i + 1) * l + j + 1]; + } + } + for(int j = 0; j < lm1; j++) { + //R_mat[lm1][j] = 0.0; + R_mat[lm1 * l + j] = 0.0; + } + for(int i = 0; i < lm1; i++) { + //R_mat[i][lm1] = sTy_vec[i]; + R_mat[i * l + lm1] = sTy_vec[i]; + } //R_mat[lm1][lm1]=sTy; - R_mat[lm1*l+lm1]=sTy; + R_mat[lm1 * l + lm1] = sTy; } + void hiopHessianInvLowRank_obsolette::updateD(const double& sTy) { - int l=D->get_size(); - double* D_vec = D->local_data(); - for(int i=0; iget_size(); + double* D_vec = D_->local_data(); + for(int i = 0; i < l-1; i++) { + D_vec[i] = D_vec[i + 1]; + } + D_vec[l - 1] = sTy; } hiopMatrixDense& hiopHessianInvLowRank_obsolette::new_S1(const hiopMatrixDense& St, const hiopMatrixDense& X) { //S1 is St*X^T (lxk), where St=S^T is lxn and X is kxn (l BFGS memory size, k number of constraints) - size_type k=X.m(), n=St.n(), l=St.m(); + size_type k = X.m(); + size_type n = St.n(); + size_type l = St.m(); #ifdef HIOP_DEEPCHECKS - assert(n==X.n()); - if(_S1!=NULL) - assert(_S1->n()==k); + assert(n == X.n()); + if(S1_ != nullptr) + assert(S1_->n() == k); #endif - if(NULL!=_S1 && _S1->n()!=l) { delete _S1; _S1=NULL; } + if(nullptr != S1_ && S1_->n() != l) { + delete S1_; + S1_ = nullptr; + } - if(NULL==_S1) _S1=LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, k); - - return *_S1; + if(nullptr == S1_) { + S1_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, k); + } + return *S1_; } hiopMatrixDense& hiopHessianInvLowRank_obsolette::new_Y1(const hiopMatrixDense& Yt, const hiopMatrixDense& X) { //Y1 is Yt*H0*X^T = Y^T*H0*X^T, where Y^T is lxn, H0 is diag nxn, X is kxn - size_type k=X.m(), n=Yt.n(), l=Yt.m(); + size_type k = X.m(); + size_type n = Yt.n(); + size_type l = Yt.m(); #ifdef HIOP_DEEPCHECKS - assert(X.n()==n); - if(_Y1!=NULL) assert(_Y1->n()==k); + assert(X.n() == n); + if(Y1_ != nullptr) { + assert(Y1_->n() == k); + } #endif - if(NULL!=_Y1 && _Y1->n()!=l) { delete _Y1; _Y1=NULL; } + if(nullptr != Y1_ && Y1_->n() != l) { + delete Y1_; + Y1_ = nullptr; + } - if(NULL==_Y1) _Y1=LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, k); - return *_Y1; + if(nullptr == Y1_) { + Y1_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, k); + } + return *Y1_; } + hiopMatrixDense& hiopHessianInvLowRank_obsolette::new_DpYtH0Y(const hiopMatrixDense& Yt) { size_type l = Yt.m(); #ifdef HIOP_DEEPCHECKS - if(_DpYtH0Y!=NULL) assert(_DpYtH0Y->m()==_DpYtH0Y->n()); + if(DpYtH0Y_ != nullptr) { + assert(DpYtH0Y_->m() == DpYtH0Y_->n()); + } #endif - if(_DpYtH0Y!=NULL && _DpYtH0Y->m()!=l) {delete _DpYtH0Y; _DpYtH0Y=NULL;} - if(_DpYtH0Y==NULL) _DpYtH0Y=LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, l); - return *_DpYtH0Y; + if(DpYtH0Y_ != nullptr && DpYtH0Y_->m() != l) { + delete DpYtH0Y_; + DpYtH0Y_ = nullptr; + } + if(DpYtH0Y_ == nullptr) { + DpYtH0Y_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, l); + } + return *DpYtH0Y_; } /* S3 = DpYtH0H * S2, where S2=R\S1. DpYtH0H is symmetric (!opt) lxl and S2 is lxk */ hiopMatrixDense& hiopHessianInvLowRank_obsolette::new_S3(const hiopMatrixDense& Left, const hiopMatrixDense& Right) { - int l=Left.m(), k=Right.n(); + int l = Left.m(); + int k = Right.n(); #ifdef HIOP_DEEPCHECKS - assert(Right.m()==l); - assert(Left.n()==l); - if(_S3!=NULL) assert(_S3->m()<=l); //< when the representation grows, = when it doesn't + assert(Right.m() == l); + assert(Left.n() == l); + if(S3_ != nullptr) { + assert(S3_->m() <= l); //< when the representation grows, = when it doesn't + } #endif - if(_S3!=NULL && _S3->m()!=l) { - delete _S3; - _S3=NULL; + if(S3_ != nullptr && S3_->m() != l) { + delete S3_; + S3_ = nullptr; } - if(_S3==NULL) { - _S3 = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, k); + if(S3_ == nullptr) { + S3_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, k); } - return *_S3; + return *S3_; } hiopVector& hiopHessianInvLowRank_obsolette::new_l_vec1(int l) { - if(_l_vec1!=NULL && _l_vec1->get_size()==l) { - return *_l_vec1; + if(l_vec1_ != nullptr && l_vec1_->get_size() == l) { + return *l_vec1_; } - if(_l_vec1!=NULL) { - delete _l_vec1; + if(l_vec1_ != nullptr) { + delete l_vec1_; } - _l_vec1 = LinearAlgebraFactory::create_vector("DEFAULT", l); - return *_l_vec1; + l_vec1_ = LinearAlgebraFactory::create_vector("DEFAULT", l); + return *l_vec1_; } hiopVector& hiopHessianInvLowRank_obsolette::new_l_vec2(int l) { - if(_l_vec2!=NULL && _l_vec2->get_size()==l) { - return *_l_vec2; + if(l_vec2_ != nullptr && l_vec2_->get_size() == l) { + return *l_vec2_; } - if(_l_vec2!=NULL) { - delete _l_vec2; + if(l_vec2_ != nullptr) { + delete l_vec2_; } - _l_vec2 = LinearAlgebraFactory::create_vector("DEFAULT", l); - return *_l_vec2; + l_vec2_ = LinearAlgebraFactory::create_vector("DEFAULT", l); + return *l_vec2_; } hiopVector& hiopHessianInvLowRank_obsolette::new_l_vec3(int l) { - if(_l_vec3!=NULL && _l_vec3->get_size()==l) return *_l_vec3; - - if(_l_vec3 != NULL) { - delete _l_vec3; + if(l_vec3_ != nullptr && l_vec3_->get_size() == l) { + return *l_vec3_; + } + if(l_vec3_ != nullptr) { + delete l_vec3_; } - _l_vec3 = LinearAlgebraFactory::create_vector("DEFAULT", l); - return *_l_vec3; + l_vec3_ = LinearAlgebraFactory::create_vector("DEFAULT", l); + return *l_vec3_; } - //#ifdef HIOP_DEEPCHECKS +//#ifdef HIOP_DEEPCHECKS // #include // using namespace std; // void hiopHessianInvLowRank_obsolette::timesVecCmn(double beta, hiopVector& y, double alpha, const hiopVector& x, bool addLogTerm) // { -// long long n=St->n(); -// assert(l_curr-1==St->m()); -// assert(y.get_size()==n); +// long long n=St_->n(); +// assert(l_curr_-1==St_->m()); +// assert(y.get_size() == n); // //we have B+=B-B*s*B*s'/(s'*B*s)+yy'/(y'*s) -// //B0 is sigma*I (and is NOT this->H0, since this->H0=(B0+Dx)^{-1}) +// //B0 is sigma*I (and is NOT this->H0_, since this->H0_ = (B0+Dx)^{-1}) // bool print=true; // if(print) { -// nlp->log->printf(hovMatrices, "---hiopHessianInvLowRank_obsolette::timesVec \n"); -// nlp->log->write("S':", *St, hovMatrices); -// nlp->log->write("Y':", *Yt, hovMatrices); -// nlp->log->write("H0:", *H0, hovMatrices); -// nlp->log->printf(hovMatrices, "sigma=%22.16e addLogTerm=%d\n", sigma, addLogTerm); -// nlp->log->printf(hovMatrices, "y=beta*y + alpha*this*x : beta=%g alpha=%g\n", beta, alpha); -// nlp->log->write("x_in:", x, hovMatrices); -// nlp->log->write("y_in:", y, hovMatrices); +// nlp_->log->printf(hovMatrices, "---hiopHessianInvLowRank_obsolette::timesVec \n"); +// nlp_->log->write("S':", *St_, hovMatrices); +// nlp_->log->write("Y':", *Yt_, hovMatrices); +// nlp_->log->write("H0_:", *H0_, hovMatrices); +// nlp_->log->printf(hovMatrices, "sigma=%22.16e addLogTerm=%d\n", sigma_, addLogTerm); +// nlp_->log->printf(hovMatrices, "y=beta*y + alpha*this*x : beta=%g alpha=%g\n", beta, alpha); +// nlp_->log->write("x_in:", x, hovMatrices); +// nlp_->log->write("y_in:", y, hovMatrices); // } -// hiopVectorPar *yk=dynamic_cast(nlp->alloc_primal_vec()); -// hiopVectorPar *sk=dynamic_cast(nlp->alloc_primal_vec()); +// hiopVectorPar *yk=dynamic_cast(nlp_->alloc_primal_vec()); +// hiopVectorPar *sk=dynamic_cast(nlp_->alloc_primal_vec()); // //allocate and compute a_k and b_k -// vector a(l_curr),b(l_curr); -// for(int k=0; k a(l_curr_),b(l_curr_); +// for(int k=0; kcopyFrom(Yt->local_data()[k]); -// sk->copyFrom(St->local_data()[k]); +// yk->copyFrom(Yt_->local_data()[k]); +// sk->copyFrom(St_->local_data()[k]); // double skTyk=yk->dotProductWith(*sk); // assert(skTyk>0); -// b[k]=dynamic_cast(nlp->alloc_primal_vec()); +// b[k]=dynamic_cast(nlp_->alloc_primal_vec()); // b[k]->copyFrom(*yk); // b[k]->scale(1/sqrt(skTyk)); -// a[k]=dynamic_cast(nlp->alloc_primal_vec()); +// a[k]=dynamic_cast(nlp_->alloc_primal_vec()); // //compute ak by an inner loop // a[k]->copyFrom(*sk); // if(addLogTerm) -// a[k]->componentDiv(*H0); +// a[k]->componentDiv(*H0_); // else -// a[k]->scale(sigma); +// a[k]->scale(sigma_); // for(int i=0; idotProductWith(*sk); // a[k]->axpy(biTsk, *b[i]); @@ -1811,15 +2057,15 @@ hiopVector& hiopHessianInvLowRank_obsolette::new_l_vec3(int l) // a[k]->scale(1/sqrt(skTak)); // } -// //new we have B= Dx+B_0 + sum{ bk bk' - ak ak' : k=0,1,...,l_curr-1} (H0=(Dx+B0)^{-1}) +// //new we have B= Dx+B_0 + sum{ bk bk' - ak ak' : k=0,1,...,l_curr_-1} (H0=(Dx+B0)^{-1}) // //compute the product with x -// //y = beta*y+alpha*H0_inv*x + alpha* sum { bk'x bk - ak'x ak : k=0,1,...,l_curr-1} +// //y = beta*y+alpha*H0_inv*x + alpha* sum { bk'x bk - ak'x ak : k=0,1,...,l_curr_-1} // y.scale(beta); // if(addLogTerm) -// y.axdzpy(alpha,x,*H0); +// y.axdzpy(alpha,x,*H0_); // else -// y.axpy(alpha*sigma, x); -// for(int k=0; kdotProductWith(x); // double akTx = a[k]->dotProductWith(x); @@ -1828,7 +2074,7 @@ hiopVector& hiopHessianInvLowRank_obsolette::new_l_vec3(int l) // } // if(print) { -// nlp->log->write("y_out:", y, hovMatrices); +// nlp_->log->write("y_out:", y, hovMatrices); // } // for(vector::iterator it=a.begin(); it!=a.end(); ++it) diff --git a/src/Optimization/hiopHessianLowRank.hpp b/src/Optimization/hiopHessianLowRank.hpp index b487ea086..710099237 100644 --- a/src/Optimization/hiopHessianLowRank.hpp +++ b/src/Optimization/hiopHessianLowRank.hpp @@ -46,6 +46,13 @@ // Lawrence Livermore National Security, LLC, and shall not be used for advertising or // product endorsement purposes. +/** + * @file hiopHessianLowRank.hpp + * + * @author Cosmin G. Petra , LLNL + * + */ + #ifndef HIOP_HESSIANLOWRANK #define HIOP_HESSIANLOWRANK @@ -91,12 +98,14 @@ namespace hiop class hiopHessianLowRank : public hiopMatrix { public: - hiopHessianLowRank(hiopNlpDenseConstraints* nlp_, int max_memory_length); + hiopHessianLowRank(hiopNlpDenseConstraints* nlp, int max_memory_length); virtual ~hiopHessianLowRank(); /* return false if the update destroys hereditary positive definitness and the BFGS update is not taken*/ - virtual bool update(const hiopIterate& x_curr, const hiopVector& grad_f_curr, - const hiopMatrix& Jac_c_curr, const hiopMatrix& Jac_d_curr); + virtual bool update(const hiopIterate& x_curr, + const hiopVector& grad_f_curr, + const hiopMatrix& Jac_c_curr_in, + const hiopMatrix& Jac_d_curr_in); /* updates the logBar diagonal term from the representation */ virtual bool updateLogBarrierDiagonal(const hiopVector& Dx); @@ -106,8 +115,10 @@ class hiopHessianLowRank : public hiopMatrix /* W = beta*W + alpha*X*inverse(this)*X^T (a more efficient version of solve) * This is performed as W = beta*W + alpha*X*(this\X^T) */ - virtual void symMatTimesInverseTimesMatTrans(double beta, hiopMatrixDense& W_, - double alpha, const hiopMatrixDense& X); + virtual void symMatTimesInverseTimesMatTrans(double beta, + hiopMatrixDense& W, + double alpha, + const hiopMatrixDense& X); #ifdef HIOP_DEEPCHECKS /* same as above but without the Dx term in H */ virtual void timesVec_noLogBarrierTerm(double beta, hiopVector& y, double alpha, const hiopVector&x); @@ -121,56 +132,69 @@ class hiopHessianLowRank : public hiopMatrix virtual void timesVec(double beta, hiopVector& y, double alpha, const hiopVector&x); /* code shared by the above two methods*/ - virtual void timesVecCmn(double beta, hiopVector& y, double alpha, const hiopVector&x, bool addLogBarTerm = false) const; + virtual void timesVecCmn(double beta, + hiopVector& y, + double alpha, + const hiopVector&x, + bool addLogBarTerm = false) const; protected: - int l_max; //max memory size - int l_curr; //number of pairs currently stored - double sigma; //initial scaling factor of identity - double sigma0; //default scaling factor of identity - int sigma_update_strategy; - double sigma_safe_min, sigma_safe_max; //min and max safety thresholds for sigma - hiopNlpDenseConstraints* nlp; - mutable std::vector a; - mutable std::vector b; - hiopVector* yk; - hiopVector* sk; + int l_max_; //max memory size + int l_curr_; //number of pairs currently stored + double sigma_; //initial scaling factor of identity + double sigma0_; //default scaling factor of identity + int sigma_update_strategy_; + double sigma_safe_min_; //min and max safety thresholds for sigma_ + double sigma_safe_max_; //min and max safety thresholds for sigma_ + hiopNlpDenseConstraints* nlp_; + mutable std::vector a_; + mutable std::vector b_; + hiopVector* yk_; + hiopVector* sk_; private: - hiopVector* DhInv; //(B0+Dk)^{-1} + hiopVector* DhInv_; //(B0+Dk)^{-1} // needed in timesVec (for residual checking in solveCompressed. // can be recomputed from DhInv decided to store it instead to avoid round-off errors - hiopVector* _Dx; - bool matrixChanged; + hiopVector* Dx_; + bool matrixChanged_; //these are matrices from the compact representation; they are updated at each iteration. // more exactly Bk=B0-[B0*St' Yt']*[St*B0*St' L]*[St*B0] // [ L' -D] [Yt ] - hiopMatrixDense *St,*Yt; //we store the transpose to easily access columns in S and T - hiopMatrixDense *L; //lower triangular from the compact representation - hiopVector* D; //diag + hiopMatrixDense* St_; + hiopMatrixDense* Yt_; //we store the transpose to easily access columns in S and T + hiopMatrixDense* L_; //lower triangular from the compact representation + hiopVector* D_; //diag //these are matrices from the representation of the inverse - hiopMatrixDense* V; + hiopMatrixDense* V_; #ifdef HIOP_DEEPCHECKS //copy of the matrix - needed to check the residual - hiopMatrixDense* _Vmat; + hiopMatrixDense* Vmat_; #endif void growL(const int& lmem_curr, const int& lmem_max, const hiopVector& YTs); - void growD(const int& l_curr, const int& l_max, const double& sTy); + void growD(const int& lmem_curr, const int& lmem_max, const double& sTy); void updateL(const hiopVector& STy, const double& sTy); void updateD(const double& sTy); //also stored are the iterate, gradient obj, and Jacobians at the previous optimization iteration - hiopIterate *_it_prev; - hiopVector *_grad_f_prev; - hiopMatrixDense *_Jac_c_prev, *_Jac_d_prev; + hiopIterate* it_prev_; + hiopVector* grad_f_prev_; + hiopMatrixDense* Jac_c_prev_; + hiopMatrixDense* Jac_d_prev_; //internal helpers void updateInternalBFGSRepresentation(); //internals buffers, mostly for MPIAll_reduce - double* _buff_kxk; // size = num_constraints^2 - double* _buff_2lxk; // size = 2 x q-Newton mem size x num_constraints - double *_buff1_lxlx3, *_buff2_lxlx3; + double* buff_kxk_; // size = num_constraints^2 + double* buff_2lxk_; // size = 2 x q-Newton mem size x num_constraints + double* buff1_lxlx3_; + double* buff2_lxlx3_; //auxiliary objects - hiopMatrixDense *_S1, *_Y1, *_lxl_mat1, *_kx2l_mat1, *_kxl_mat1; //preallocated matrices + hiopMatrixDense* S1_; + hiopMatrixDense* Y1_; + hiopMatrixDense* lxl_mat1_; + hiopMatrixDense* kx2l_mat1_; + hiopMatrixDense* kxl_mat1_; //preallocated matrices + //holds X*D*S hiopMatrixDense& new_S1(const hiopMatrixDense& X, const hiopMatrixDense& St); //holds X*D*Y @@ -179,62 +203,79 @@ class hiopHessianLowRank : public hiopMatrix hiopMatrixDense& new_kxl_mat1 (int k, int l); hiopMatrixDense& new_kx2l_mat1(int k, int l); - hiopVector *_l_vec1, *_l_vec2, *_n_vec1, *_n_vec2, *_2l_vec1; + hiopVector* l_vec1_; + hiopVector* l_vec2_; + hiopVector* n_vec1_; + hiopVector* n_vec2_; + hiopVector* twol_vec1_; hiopVector& new_l_vec1(int l); hiopVector& new_l_vec2(int l); inline hiopVector& new_n_vec1(size_type n) { #ifdef HIOP_DEEPCHECKS - assert(_n_vec1!=NULL); - assert(_n_vec1->get_size()==n); + assert(n_vec1_ != nullptr); + assert(n_vec1_->get_size()==n); #endif - return *_n_vec1; + return *n_vec1_; } inline hiopVector& new_n_vec2(size_type n) { #ifdef HIOP_DEEPCHECKS - assert(_n_vec2!=NULL); - assert(_n_vec2->get_size()==n); + assert(n_vec2_ != nullptr); + assert(n_vec2_->get_size()==n); #endif - return *_n_vec2; + return *n_vec2_; } inline hiopVector& new_2l_vec1(int l) { - if(_2l_vec1!=NULL && _2l_vec1->get_size()==2*l) return *_2l_vec1; - if(_2l_vec1!=NULL) delete _2l_vec1; - _2l_vec1=LinearAlgebraFactory::create_vector(nlp->options->GetString("mem_space"), 2*l); - return *_2l_vec1; + if(twol_vec1_ != nullptr && twol_vec1_->get_size() == 2 * l) { + return *twol_vec1_; + } + if(twol_vec1_ != nullptr) { + delete twol_vec1_; + } + twol_vec1_ = LinearAlgebraFactory::create_vector(nlp_->options->GetString("mem_space"), 2 * l); + return *twol_vec1_; } private: //utilities /* symmetric multiplication W = beta*W + alpha*X*Diag*X^T */ - static void symmMatTimesDiagTimesMatTrans_local(double beta, hiopMatrixDense& W_, - double alpha, const hiopMatrixDense& X_, - const hiopVector& d); + static void symmMatTimesDiagTimesMatTrans_local(double beta, + hiopMatrixDense& W, + double alpha, + const hiopMatrixDense& X, + const hiopVector& d); /* W=S*Diag*X^T */ - static void matTimesDiagTimesMatTrans_local(hiopMatrixDense& W, const hiopMatrixDense& S, - const hiopVector& d, const hiopMatrixDense& X); + static void matTimesDiagTimesMatTrans_local(hiopMatrixDense& W, + const hiopMatrixDense& S, + const hiopVector& d, + const hiopMatrixDense& X); + /* members and utilities related to V matrix: factorization and solve */ - hiopVector *_V_work_vec; - int _V_ipiv_size; int* _V_ipiv_vec; + hiopVector* V_work_vec_; + int V_ipiv_size_; + int* V_ipiv_vec_; + void factorizeV(); void solveWithV(hiopVector& rhs_s, hiopVector& rhs_y); void solveWithV(hiopMatrixDense& rhs); private: hiopHessianLowRank() {}; hiopHessianLowRank(const hiopHessianLowRank&) {}; - hiopHessianLowRank& operator=(const hiopHessianLowRank&) {return *this;}; + hiopHessianLowRank& operator=(const hiopHessianLowRank&) { + return *this; + } public: /* methods that need to be implemented as the class inherits from hiopMatrix*/ public: virtual hiopMatrix* alloc_clone() const { assert(false && "not provided because it is not needed"); - return NULL; + return nullptr; } virtual hiopMatrix* new_copy() const { assert(false && "not provided because it is not needed"); - return NULL; + return nullptr; } virtual void setToZero() @@ -249,8 +290,10 @@ class hiopHessianLowRank : public hiopMatrix void timesVec(double beta, hiopVector& y, double alpha, const hiopVector&x) const; /** y = beta * y + alpha * this^T * x */ - virtual void transTimesVec(double beta, hiopVector& y, - double alpha, const hiopVector& x ) const + virtual void transTimesVec(double beta, + hiopVector& y, + double alpha, + const hiopVector& x ) const { assert(false && "not provided because it is not needed"); } @@ -270,7 +313,7 @@ class hiopHessianLowRank : public hiopMatrix { assert(false && "not provided because it is not needed"); } - virtual void addDiagonal(const double& alpha, const hiopVector& d_) + virtual void addDiagonal(const double& alpha, const hiopVector& d) { assert(false && "not provided because it is not needed"); } @@ -278,15 +321,18 @@ class hiopHessianLowRank : public hiopMatrix { assert(false && "not provided because it is not needed"); } - virtual void addSubDiagonal(const double& alpha, index_type start, const hiopVector& d_) + virtual void addSubDiagonal(const double& alpha, index_type start, const hiopVector& d) { assert(false && "not provided because it is not needed"); } /* add to the diagonal of 'this' (destination) starting at 'start_on_dest_diag' elements of - * 'd_' (source) starting at index 'start_on_src_vec'. The number of elements added is 'num_elems' - * when num_elems>=0, or the remaining elems on 'd_' starting at 'start_on_src_vec'. */ - virtual void addSubDiagonal(int start_on_dest_diag, const double& alpha, - const hiopVector& d_, int start_on_src_vec, int num_elems=-1) + * 'd_src' (source) starting at index 'start_on_src_vec'. The number of elements added is 'num_elems' + * when num_elems>=0, or the remaining elems on 'd_src' starting at 'start_on_src_vec'. */ + virtual void addSubDiagonal(int start_on_dest_diag, + const double& alpha, + const hiopVector& d_src, + int start_on_src_vec, + int num_elems=-1) { assert(false && "not needed / implemented"); } @@ -310,8 +356,7 @@ class hiopHessianLowRank : public hiopMatrix { assert(false && "not needed; should not be used"); } - void addUpperTriangleToSymDenseMatrixUpperTriangle(int diag_start, - double alpha, hiopMatrixDense& W) const + void addUpperTriangleToSymDenseMatrixUpperTriangle(int diag_start, double alpha, hiopMatrixDense& W) const { assert(false && "not needed; should not be used"); } @@ -350,7 +395,7 @@ class hiopHessianLowRank : public hiopMatrix * given by the value of 'maxRows' will be printed. If this value is negative, all * elements will be printed. */ - virtual void print(FILE* f=NULL, const char* msg=NULL, int maxRows=-1, int maxCols=-1, int rank=-1) const + virtual void print(FILE* f=nullptr, const char* msg=nullptr, int maxRows=-1, int maxCols=-1, int rank=-1) const { assert(false && "not provided because it is not needed"); } @@ -382,8 +427,10 @@ class hiopHessianInvLowRank_obsolette : public hiopHessianLowRank { public: hiopHessianInvLowRank_obsolette(hiopNlpDenseConstraints* nlp, int max_memory_length); - virtual bool update(const hiopIterate& x_curr, const hiopVector& grad_f_curr, - const hiopMatrix& Jac_c_curr, const hiopMatrix& Jac_d_curr); + virtual bool update(const hiopIterate& x_curr, + const hiopVector& grad_f_curr, + const hiopMatrix& Jac_c_curr_in, + const hiopMatrix& Jac_d_curr_in); virtual bool updateLogBarrierDiagonal(const hiopVector& Dx); @@ -397,19 +444,23 @@ class hiopHessianInvLowRank_obsolette : public hiopHessianLowRank /* W = beta*W + alpha*X*this*X^T * ! make sure this is called before 'apply' */ - virtual void symmetricTimesMat(double beta, hiopMatrixDense& W, - double alpha, const hiopMatrixDense& X); + virtual void symmetricTimesMat(double beta, hiopMatrixDense& W, double alpha, const hiopMatrixDense& X); virtual ~hiopHessianInvLowRank_obsolette(); private: //internal methods /* symmetric multiplication W = beta*W + alpha*X*Diag*X^T */ - static void symmMatTimesDiagTimesMatTrans_local(double beta, hiopMatrixDense& W_, - double alpha, const hiopMatrixDense& X_, - const hiopVector& d); + static void symmMatTimesDiagTimesMatTrans_local(double beta, + hiopMatrixDense& W, + double alpha, + const hiopMatrixDense& X, + const hiopVector& d); /* W=S*Diag*X^T */ - static void matTimesDiagTimesMatTrans_local(hiopMatrixDense& W, const hiopMatrixDense& S, const hiopVector& d, const hiopMatrixDense& X); + static void matTimesDiagTimesMatTrans_local(hiopMatrixDense& W, + const hiopMatrixDense& S, + const hiopVector& d, + const hiopMatrixDense& X); /* rhs = R \ rhs, where R is upper triangular lxl and rhs is lx */ static void triangularSolve(const hiopMatrixDense& R, hiopMatrixDense& rhs); @@ -422,50 +473,59 @@ class hiopHessianInvLowRank_obsolette : public hiopHessianLowRank void updateR(const hiopVector& STy, const double& sTy); void updateD(const double& sTy); private: - hiopVector* H0; - hiopMatrixDense *St,*Yt; //we store the transpose to easily access columns in S and T - hiopMatrixDense *R; - hiopVector* D; + hiopVector* H0_; + hiopMatrixDense* St_; + hiopMatrixDense* Yt_; //we store the transpose to easily access columns in S and T + hiopMatrixDense* R_; + hiopVector* D_; - int sigma_update_strategy; - double sigma_safe_min, sigma_safe_max; + int sigma_update_strategy_; + double sigma_safe_min_; + double sigma_safe_max_; //also stored are the iterate, gradient obj, and Jacobians at the previous iterations - hiopIterate *_it_prev; - hiopVector *_grad_f_prev; - hiopMatrixDense *_Jac_c_prev, *_Jac_d_prev; + hiopIterate* it_prev_; + hiopVector* grad_f_prev_; + hiopMatrixDense* Jac_c_prev_; + hiopMatrixDense* Jac_d_prev_; //internals buffers - double* _buff_kxk; // size = num_constraints^2 - double* _buff_lxk; // size = q-Newton mem size x num_constraints - double* _buff_lxl; + double* buff_kxk_; // size = num_constraints^2 + double* buff_lxk_; // size = q-Newton mem size x num_constraints + double* buff_lxl_; //auxiliary objects - hiopMatrixDense *_S1, *_Y1, *_DpYtH0Y; //aux matrices to hold St*X, Yt*H0*X, and D+Y^T*H0*Y in symmetricTimesMat + hiopMatrixDense* S1_; //aux matrices to hold St*X, Yt*H0*X, and D+Y^T*H0*Y in symmetricTimesMat + hiopMatrixDense* Y1_; + hiopMatrixDense* DpYtH0Y_; hiopMatrixDense& new_S1(const hiopMatrixDense& St, const hiopMatrixDense& X); hiopMatrixDense& new_Y1(const hiopMatrixDense& Yt, const hiopMatrixDense& X); hiopMatrixDense& new_DpYtH0Y(const hiopMatrixDense& Yt); //similar for S3=DpYtH0Y*S2 - hiopMatrixDense *_S3; + hiopMatrixDense* S3_; hiopMatrixDense& new_S3(const hiopMatrixDense& Left, const hiopMatrixDense& Right); - hiopVector *_l_vec1, *_l_vec2, *_l_vec3, *_n_vec1, *_n_vec2; + hiopVector* l_vec1_; + hiopVector* l_vec2_; + hiopVector* l_vec3_; + hiopVector* n_vec1_; + hiopVector* n_vec2_; hiopVector& new_l_vec1(int l); hiopVector& new_l_vec2(int l); hiopVector& new_l_vec3(int l); inline hiopVector& new_n_vec1(size_type n) { #ifdef HIOP_DEEPCHECKS - assert(_n_vec1!=NULL); - assert(_n_vec1->get_size()==n); + assert(n_vec1_ != nullptr); + assert(n_vec1_->get_size()==n); #endif - return *_n_vec1; + return *n_vec1_; } inline hiopVector& new_n_vec2(size_type n) { #ifdef HIOP_DEEPCHECKS - assert(_n_vec2!=NULL); - assert(_n_vec2->get_size()==n); + assert(n_vec2_ != nullptr); + assert(n_vec2_->get_size()==n); #endif - return *_n_vec2; + return *n_vec2_; } };