diff --git a/src/Optimization/CMakeLists.txt b/src/Optimization/CMakeLists.txt index 30b7bdc6..22fc3226 100644 --- a/src/Optimization/CMakeLists.txt +++ b/src/Optimization/CMakeLists.txt @@ -5,7 +5,8 @@ set(hiopOptimization_SRC hiopResidual.cpp hiopFilter.cpp hiopAlgFilterIPM.cpp - hiopKKTLinSys.cpp + hiopKKTLinSys.cpp + KktLinSysLowRank.cpp hiopKKTLinSysMDS.cpp hiopHessianLowRank.cpp hiopDualsUpdater.cpp diff --git a/src/Optimization/KktLinSysLowRank.cpp b/src/Optimization/KktLinSysLowRank.cpp new file mode 100644 index 00000000..bbb9c3bb --- /dev/null +++ b/src/Optimization/KktLinSysLowRank.cpp @@ -0,0 +1,505 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC. +// Produced at the Lawrence Livermore National Laboratory (LLNL). +// LLNL-CODE-742473. All rights reserved. +// +// This file is part of HiOp. For details, see https://github.com/LLNL/hiop. HiOp +// is released under the BSD 3-clause license (https://opensource.org/licenses/BSD-3-Clause). +// Please also read "Additional BSD Notice" below. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// i. Redistributions of source code must retain the above copyright notice, this list +// of conditions and the disclaimer below. +// ii. Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the disclaimer (as noted below) in the documentation and/or +// other materials provided with the distribution. +// iii. Neither the name of the LLNS/LLNL nor the names of its contributors may be used to +// endorse or promote products derived from this software without specific prior written +// permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES +// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +// SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED +// AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, +// EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Additional BSD Notice +// 1. This notice is required to be provided under our contract with the U.S. Department +// of Energy (DOE). This work was produced at Lawrence Livermore National Laboratory under +// Contract No. DE-AC52-07NA27344 with the DOE. +// 2. Neither the United States Government nor Lawrence Livermore National Security, LLC +// nor any of their employees, makes any warranty, express or implied, or assumes any +// liability or responsibility for the accuracy, completeness, or usefulness of any +// information, apparatus, product, or process disclosed, or represents that its use would +// not infringe privately-owned rights. +// 3. Also, reference herein to any specific commercial products, process, or services by +// trade name, trademark, manufacturer or otherwise does not necessarily constitute or +// imply its endorsement, recommendation, or favoring by the United States Government or +// Lawrence Livermore National Security, LLC. The views and opinions of authors expressed +// herein do not necessarily state or reflect those of the United States Government or +// Lawrence Livermore National Security, LLC, and shall not be used for advertising or +// product endorsement purposes. + +/** + * @file KktLinSysLowRank.cpp + * + * @author Cosmin G. Petra , LLNL + * + */ + +#include "KktLinSysLowRank.hpp" + +namespace hiop +{ + +KktLinSysLowRank::KktLinSysLowRank(hiopNlpFormulation* nlp) + : hiopKKTLinSysCompressedXYcYd(nlp) +{ + nlpD = dynamic_cast(nlp_); + + kxn_mat_ = nlpD->alloc_multivector_primal(nlpD->m()); + assert("DEFAULT" == toupper(nlpD->options->GetString("mem_space"))); + N = LinearAlgebraFactory::create_matrix_dense(nlpD->options->GetString("mem_space"), + nlpD->m(), + nlpD->m()); +#ifdef HIOP_DEEPCHECKS + Nmat=N->alloc_clone(); +#endif + k_vec1_ = dynamic_cast(nlpD->alloc_dual_vec()); +} + +KktLinSysLowRank::~KktLinSysLowRank() +{ + if(N) delete N; +#ifdef HIOP_DEEPCHECKS + if(Nmat) delete Nmat; +#endif + if(kxn_mat_) delete kxn_mat_; + if(k_vec1_) delete k_vec1_; +} + +bool KktLinSysLowRank::update(const hiopIterate* iter, + const hiopVector* grad_f, + const hiopMatrixDense* Jac_c, + const hiopMatrixDense* Jac_d, + hiopHessianLowRank* Hess) +{ + nlp_->runStats.tmSolverInternal.start(); + + iter_=iter; + grad_f_ = dynamic_cast(grad_f); + Jac_c_ = Jac_c; Jac_d_ = Jac_d; + //Hess = dynamic_cast(Hess_); + Hess_=HessLowRank=Hess; + + //compute the diagonals + //Dx=(Sxl)^{-1}Zl + (Sxu)^{-1}Zu + Dx_->setToZero(); + Dx_->axdzpy_w_pattern(1.0, *iter_->zl, *iter_->sxl, nlp_->get_ixl()); + Dx_->axdzpy_w_pattern(1.0, *iter_->zu, *iter_->sxu, nlp_->get_ixu()); + nlp_->log->write("Dx in KKT", *Dx_, hovMatrices); + + HessLowRank->updateLogBarrierDiagonal(*Dx_); + + //Dd=(Sdl)^{-1}Vu + (Sdu)^{-1}Vu + Dd_inv_->setToZero(); + Dd_inv_->axdzpy_w_pattern(1.0, *iter_->vl, *iter_->sdl, nlp_->get_idl()); + Dd_inv_->axdzpy_w_pattern(1.0, *iter_->vu, *iter_->sdu, nlp_->get_idu()); +#ifdef HIOP_DEEPCHECKS + assert(true==Dd_inv_->allPositive()); +#endif + Dd_->copyFrom(*Dd_inv_); + Dd_inv_->invert(); + + nlp_->runStats.tmSolverInternal.stop(); + + nlp_->log->write("Dd_inv in KKT", *Dd_inv_, hovMatrices); + return true; +} + + +/* Solves the system corresponding to directions for x, yc, and yd, namely + * [ H_BFGS + Dx Jc^T Jd^T ] [ dx] [ rx ] + * [ Jc 0 0 ] [dyc] = [ ryc ] + * [ Jd 0 -Dd^{-1}] [dyd] [ ryd ] + * + * This is done by forming and solving + * [ Jc*(H+Dx)^{-1}*Jc^T Jc*(H+Dx)^{-1}*Jd^T ] [dyc] = [ Jc(H+Dx)^{-1} rx - ryc ] + * [ Jd*(H+Dx)^{-1}*Jc^T Jd*(H+Dx)^{-1}*Jd^T + Dd^{-1}] [dyd] [ Jd(H+dx)^{-1} rx - ryd ] + * and then solving for dx from + * dx = - (H+Dx)^{-1}*(Jc^T*dyc+Jd^T*dyd - rx) + * + * Note that ops H+Dx are provided by hiopHessianLowRank + */ +bool KktLinSysLowRank::solveCompressed(hiopVector& rx, + hiopVector& ryc, + hiopVector& ryd, + hiopVector& dx, + hiopVector& dyc, + hiopVector& dyd) +{ +#ifdef HIOP_DEEPCHECKS + //some outputing + nlp_->log->write("KKT Low rank: solve compressed RHS", hovIteration); + nlp_->log->write(" rx: ", rx, hovIteration); + nlp_->log->write(" ryc: ", ryc, hovIteration); + nlp_->log->write(" ryd: ", ryd, hovIteration); + nlp_->log->write(" Jc: ", *Jac_c_, hovMatrices); + nlp_->log->write(" Jd: ", *Jac_d_, hovMatrices); + nlp_->log->write(" Dd_inv: ", *Dd_inv_, hovMatrices); + assert(Dd_inv_->isfinite_local() && "Something bad happened: nan or inf value"); +#endif + + hiopMatrixDense& J = *kxn_mat_; + const hiopMatrixDense* Jac_c_de = dynamic_cast(Jac_c_); assert(Jac_c_de); + const hiopMatrixDense* Jac_d_de = dynamic_cast(Jac_d_); assert(Jac_d_de); + J.copyRowsFrom(*Jac_c_de, nlp_->m_eq(), 0); //!opt + J.copyRowsFrom(*Jac_d_de, nlp_->m_ineq(), nlp_->m_eq());//!opt + + //N = J*(Hess\J') + //Hess->symmetricTimesMat(0.0, *N, 1.0, J); + HessLowRank->symMatTimesInverseTimesMatTrans(0.0, *N, 1.0, J); + + //subdiag of N += 1., Dd_inv + N->addSubDiagonal(1., nlp_->m_eq(), *Dd_inv_); +#ifdef HIOP_DEEPCHECKS + assert(J.isfinite()); + nlp_->log->write("solveCompressed: N is", *N, hovMatrices); + nlp_->log->write("solveCompressed: rx is", rx, hovMatrices); + nlp_->log->printf(hovLinAlgScalars, "inf norm of Dd_inv is %g\n", Dd_inv_->infnorm()); + N->assertSymmetry(1e-10); +#endif + + //compute the rhs of the lin sys involving N + // 1. first compute (H+Dx)^{-1} rx_tilde and store it temporarily in dx + HessLowRank->solve(rx, dx); +#ifdef HIOP_DEEPCHECKS + assert(rx.isfinite_local() && "Something bad happened: nan or inf value"); + assert(dx.isfinite_local() && "Something bad happened: nan or inf value"); +#endif + + // 2 . then rhs = [ Jc(H+Dx)^{-1}*rx - ryc ] + // [ Jd(H+dx)^{-1}*rx - ryd ] + hiopVector& rhs=*k_vec1_; + rhs.copyFromStarting(0, ryc); + rhs.copyFromStarting(nlp_->m_eq(), ryd); + J.timesVec(-1.0, rhs, 1.0, dx); + +#ifdef HIOP_DEEPCHECKS + nlp_->log->write("solveCompressed: dx sol is", dx, hovMatrices); + nlp_->log->write("solveCompressed: rhs for N is", rhs, hovMatrices); + Nmat->copyFrom(*N); + hiopVector* r=rhs.new_copy(); //save the rhs to check the norm of the residual +#endif + + // + //solve N * dyc_dyd = rhs + // + int ierr = solveWithRefin(*N,rhs); + //int ierr = solve(*N,rhs); + + hiopVector& dyc_dyd= rhs; + dyc_dyd.copyToStarting(0, dyc); + dyc_dyd.copyToStarting(nlp_->m_eq(), dyd); + + //now solve for dx = - (H+Dx)^{-1}*(Jc^T*dyc+Jd^T*dyd - rx) + //first rx = -(Jc^T*dyc+Jd^T*dyd - rx) + J.transTimesVec(1.0, rx, -1.0, dyc_dyd); + //then dx = (H+Dx)^{-1} rx + HessLowRank->solve(rx, dx); + +#ifdef HIOP_DEEPCHECKS + //some outputing + nlp_->log->write("KKT Low rank: solve compressed SOL", hovIteration); + nlp_->log->write(" dx: ", dx, hovIteration); + nlp_->log->write(" dyc: ", dyc, hovIteration); + nlp_->log->write(" dyd: ", dyd, hovIteration); + delete r; +#endif + + return ierr==0; +} + +int KktLinSysLowRank::solveWithRefin(hiopMatrixDense& M, hiopVector& rhs) +{ + // 1. Solve dposvx (solve + equilibrating + iterative refinement + forward and backward error estimates) + // 2. Check the residual norm + // 3. If residual norm is not small enough, then perform iterative refinement. This is because dposvx + // does not always provide a small enough residual since it stops (possibly without refinement) based on + // the forward and backward estimates + + int N=M.n(); + if(N<=0) return 0; + + hiopMatrixDense* Aref = M.new_copy(); + hiopVector* rhsref = rhs.new_copy(); + + char FACT='E'; + char UPLO='L'; + + int NRHS=1; + double* A=M.local_data(); + int LDA=N; + double* AF=new double[N*N]; + int LDAF=N; + char EQUED='N'; //it is an output if FACT='E' + double* S = new double[N]; + double* B = rhs.local_data(); + int LDB=N; + double* X = new double[N]; + int LDX = N; + double RCOND, FERR, BERR; + double* WORK = new double[3*N]; + int* IWORK = new int[N]; + int INFO; + + // + // 1. solve + // + DPOSVX(&FACT, &UPLO, &N, &NRHS, A, &LDA, AF, &LDAF, &EQUED, S, B, &LDB, X, &LDX, &RCOND, &FERR, &BERR, WORK, IWORK, &INFO); + //printf("INFO ===== %d RCOND=%g FERR=%g BERR=%g EQUED=%c\n", INFO, RCOND, FERR, BERR, EQUED); + // + // 2. check residual + // + hiopVector* x = rhs.alloc_clone(); + hiopVector* dx = rhs.alloc_clone(); + hiopVector* resid = rhs.alloc_clone(); + int nIterRefin=0;double nrmResid; + int info; + const int MAX_ITER_REFIN=3; + while(true) { + x->copyFrom(X); + resid->copyFrom(*rhsref); + Aref->timesVec(1.0, *resid, -1.0, *x); + + nlp_->log->write("resid", *resid, hovLinAlgScalars); + + nrmResid= resid->infnorm(); + nlp_->log->printf(hovScalars, "KktLinSysLowRank::solveWithRefin iterrefin=%d residual norm=%g\n", nIterRefin, nrmResid); + + if(nrmResid<1e-8) break; + + if(nIterRefin>=MAX_ITER_REFIN) { + nlp_->log->write("N", *Aref, hovMatrices); + nlp_->log->write("sol", *x, hovMatrices); + nlp_->log->write("rhs", *rhsref, hovMatrices); + + nlp_->log->printf(hovWarning, + "KktLinSysLowRank::solveWithRefin reduced residual to ONLY (inf-norm) %g after %d iterative refinements\n", + nrmResid, + nIterRefin); + break; + //assert(false && "too many refinements"); + } + if(0) { //iter refin based on symmetric indefinite factorization+solve + + + int _V_ipiv_vec[1000]; + double _V_work_vec[1000]; + int lwork=1000; + M.copyFrom(*Aref); + DSYTRF(&UPLO, &N, M.local_data(), &LDA, _V_ipiv_vec, _V_work_vec, &lwork, &info); + assert(info==0); + DSYTRS(&UPLO, &N, &NRHS, M.local_data(), &LDA, _V_ipiv_vec, resid->local_data(), &LDB, &info); + assert(info==0); + } else { + //iter refin based on symmetric positive definite factorization+solve + M.copyFrom(*Aref); + DPOTRF(&UPLO, &N, M.local_data(), &LDA, &info); + if(info>0) { + nlp_->log->printf(hovError, + "KktLinSysLowRank::factorizeMat: dpotrf (Chol fact) detected %d minor being indefinite.\n", + info); + } else { + if(info<0) { + nlp_->log->printf(hovError, "KktLinSysLowRank::factorizeMat: dpotrf returned error %d\n", info); + } + } + + DPOTRS(&UPLO,&N, &NRHS, M.local_data(), &LDA, resid->local_data(), &LDA, &info); + if(info<0) { + nlp_->log->printf(hovError, "KktLinSysLowRank::solveWithFactors: dpotrs returned error %d\n", info); + } + } + + dx->copyFrom(*resid); + x->axpy(1., *dx); + + nIterRefin++; + } + rhs.copyFrom(*x); + delete[] AF; + delete[] S; + delete[] X; + delete[] WORK; + delete[] IWORK; + delete Aref; + delete rhsref; + delete x; + delete dx; + delete resid; + + return 0; +} + +int KktLinSysLowRank::solve(hiopMatrixDense& M, hiopVector& rhs) +{ + char FACT='E'; + char UPLO='L'; + int N=M.n(); + int NRHS=1; + double* A=M.local_data(); + int LDA=N; + double* AF=new double[N*N]; + int LDAF=N; + char EQUED='N'; //it is an output if FACT='E' + double* S = new double[N]; + double* B = rhs.local_data(); + int LDB=N; + double* X = new double[N]; + int LDX = N; + double RCOND, FERR, BERR; + double* WORK = new double[3*N]; + int* IWORK = new int[N]; + int INFO; + + DPOSVX(&FACT, &UPLO, &N, &NRHS, A, &LDA, AF, &LDAF, &EQUED, S, B, &LDB, X, &LDX, &RCOND, &FERR, &BERR, WORK, IWORK, &INFO); + + rhs.copyFrom(S); + nlp_->log->write("Scaling S", rhs, hovSummary); + + rhs.copyFrom(X); + delete [] AF; + delete [] S; + delete [] X; + delete [] WORK; + delete [] IWORK; + return 0; +} + +/* this code works fine but requires xblas +int KktLinSysLowRank::solveWithRefin(hiopMatrixDense& M, hiopVectorPar& rhs) +{ + char FACT='E'; + char UPLO='L'; + int N=M.n(); + int NRHS=1; + double* A=M.local_buffer(); + int LDA=N; + double* AF=new double[N*N]; + int LDAF=N; + char EQUED='N'; //it is an output if FACT='E' + double* S = new double[N]; + double* B = rhs.local_data(); + int LDB=N; + double* X = new double[N]; + int LDX = N; + double RCOND, BERR; + double RPVGRW; //Reciprocal pivot growth + int N_ERR_BNDS=3; + double* ERR_BNDS_NORM = new double[NRHS*N_ERR_BNDS]; + double* ERR_BNDS_COMP = new double[NRHS*N_ERR_BNDS]; + int NPARAMS=3; + double PARAMS[NPARAMS]; + PARAMS[0]=1.0; //Use the extra-precise refinement algorithm + PARAMS[1]=3.0; //Maximum number of residual computations allowed for refinement + PARAMS[2]=1.0; //attempt to find a solution with small componentwise + double* WORK = new double[4*N]; + int* IWORK = new int[N]; + int INFO; + + dposvxx_(&FACT, &UPLO, &N, &NRHS, + A, &LDA, + AF, &LDAF, + &EQUED, + S, + B, &LDB, + X, &LDX, + &RCOND, &RPVGRW, &BERR, + &N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, + &NPARAMS, PARAMS, + WORK, IWORK, + &INFO); + + //rhs.copyFrom(S); + //nlp_->log->write("Scaling S", rhs, hovSummary); + + //M.copyFrom(AF); + //nlp_->log->write("Factoriz ", M, hovSummary); + + printf("INFO ===== %d RCOND=%g RPVGRW=%g BERR=%g EQUED=%c\n", INFO, RCOND, RPVGRW, BERR, EQUED); + printf(" ERR_BNDS_NORM=%g %g %g ERR_BNDS_COMP=%g %g %g \n", ERR_BNDS_NORM[0], ERR_BNDS_NORM[1], ERR_BNDS_NORM[2], ERR_BNDS_COMP[0], ERR_BNDS_COMP[1], ERR_BNDS_COMP[2]); + printf(" PARAMS=%g %g %g \n", PARAMS[0], PARAMS[1], PARAMS[2]); + + + rhs.copyFrom(X); + delete [] AF; + delete [] S; + delete [] X; + delete [] ERR_BNDS_NORM; + delete [] ERR_BNDS_COMP; + delete [] WORK; + delete [] IWORK; + return 0; +} +*/ + +#ifdef HIOP_DEEPCHECKS + +double KktLinSysLowRank::errorCompressedLinsys(const hiopVector& rx, + const hiopVector& ryc, + const hiopVector& ryd, + const hiopVector& dx, + const hiopVector& dyc, + const hiopVector& dyd) +{ + nlp_->log->printf(hovLinAlgScalars, "KktLinSysLowRank::errorCompressedLinsys residuals norm:\n"); + + double derr = -1.; + double aux; + hiopVector* RX = rx.new_copy(); + //RX=rx-H*dx-J'c*dyc-J'*dyd + HessLowRank->timesVec(1.0, *RX, -1.0, dx); + //RX->axzpy(-1.0,*Dx,dx); + Jac_c_->transTimesVec(1.0, *RX, -1.0, dyc); + Jac_d_->transTimesVec(1.0, *RX, -1.0, dyd); + aux = RX->twonorm(); + derr = fmax(derr,aux); + nlp_->log->printf(hovLinAlgScalars, " >>> rx=%g\n", aux); + delete RX; + + hiopVector* RC = ryc.new_copy(); + Jac_c_->timesVec(1.0,*RC, -1.0,dx); + aux = RC->twonorm(); + derr = fmax(derr,aux); + nlp_->log->printf(hovLinAlgScalars, " >>> ryc=%g\n", aux); + delete RC; + + hiopVector* RD = ryd.new_copy(); + Jac_d_->timesVec(1.0,*RD, -1.0, dx); + RD->axzpy(1.0, *Dd_inv_, dyd); + aux = RD->twonorm(); + derr=fmax(derr,aux); + nlp_->log->printf(hovLinAlgScalars, " >>> ryd=%g\n", aux); + delete RD; + + return derr; +} + +double KktLinSysLowRank::solveError(const hiopMatrixDense& M, const hiopVector& x, hiopVector& rhs) +{ + double relError; + M.timesVec(1.0,rhs,-1.0,x); + double resnorm = rhs.infnorm(); + + relError=resnorm;// / (1+rhsnorm); + return relError; +} +#endif + +} //end namespace diff --git a/src/Optimization/KktLinSysLowRank.hpp b/src/Optimization/KktLinSysLowRank.hpp new file mode 100644 index 00000000..ab8f7d94 --- /dev/null +++ b/src/Optimization/KktLinSysLowRank.hpp @@ -0,0 +1,175 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC. +// Produced at the Lawrence Livermore National Laboratory (LLNL). +// Written by Cosmin G. Petra, petra1@llnl.gov. +// LLNL-CODE-742473. All rights reserved. +// +// This file is part of HiOp. For details, see https://github.com/LLNL/hiop. HiOp +// is released under the BSD 3-clause license (https://opensource.org/licenses/BSD-3-Clause). +// Please also read “Additional BSD Notice” below. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// i. Redistributions of source code must retain the above copyright notice, this list +// of conditions and the disclaimer below. +// ii. Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the disclaimer (as noted below) in the documentation and/or +// other materials provided with the distribution. +// iii. Neither the name of the LLNS/LLNL nor the names of its contributors may be used to +// endorse or promote products derived from this software without specific prior written +// permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES +// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +// SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED +// AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, +// EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Additional BSD Notice +// 1. This notice is required to be provided under our contract with the U.S. Department +// of Energy (DOE). This work was produced at Lawrence Livermore National Laboratory under +// Contract No. DE-AC52-07NA27344 with the DOE. +// 2. Neither the United States Government nor Lawrence Livermore National Security, LLC +// nor any of their employees, makes any warranty, express or implied, or assumes any +// liability or responsibility for the accuracy, completeness, or usefulness of any +// information, apparatus, product, or process disclosed, or represents that its use would +// not infringe privately-owned rights. +// 3. Also, reference herein to any specific commercial products, process, or services by +// trade name, trademark, manufacturer or otherwise does not necessarily constitute or +// imply its endorsement, recommendation, or favoring by the United States Government or +// Lawrence Livermore National Security, LLC. The views and opinions of authors expressed +// herein do not necessarily state or reflect those of the United States Government or +// Lawrence Livermore National Security, LLC, and shall not be used for advertising or +// product endorsement purposes. + +/** + * @file KktLinSysLowRank.cpp + * + * @author Cosmin G. Petra , LLNL + * + */ + +#ifndef HIOP_KKTLINSYSY_LOWRANK +#define HIOP_KKTLINSYSY_LOWRANK + +#include "hiopKKTLinSys.hpp" + +namespace hiop +{ + +/** + * @brief Encapsulates solves with the KKT system of IPM filter. + * + * This class is for problems where the Hessian of the Lagrangian is a or is approximated + * by low-rank matrix plus a multiple of identity and the number of the constraints is not + * too large. + * + * It works with Hessian being a HessianLowRank class and the constraints Jacobian being + * hiopMatrixDense. + * + * This class solves the XYcYd compression of the full KKT. See solveCompressed method + * for details on the approach used to solve the linear system. + */ + +class KktLinSysLowRank : public hiopKKTLinSysCompressedXYcYd +{ +public: + KktLinSysLowRank(hiopNlpFormulation* nlp); + virtual ~KktLinSysLowRank(); + + /// @brief Updates the KKT system with new info at current iteration + bool update(const hiopIterate* iter, + const hiopVector* grad_f, + const hiopMatrix* Jac_c, + const hiopMatrix* Jac_d, + hiopMatrix* Hess) + { + const hiopMatrixDense* Jac_c_ = dynamic_cast(Jac_c); + const hiopMatrixDense* Jac_d_ = dynamic_cast(Jac_d); + hiopHessianLowRank* Hess_ = dynamic_cast(Hess); + if(Jac_c_==nullptr || Jac_d_==nullptr || Hess_==nullptr) { + assert(false); + return false; + } + return update(iter, grad_f_, Jac_c_, Jac_d_, Hess_); + } + + /// @brief Updates the KKT system with new info at current iteration + virtual bool update(const hiopIterate* iter, + const hiopVector* grad_f, + const hiopMatrixDense* Jac_c, + const hiopMatrixDense* Jac_d, + hiopHessianLowRank* Hess); + + virtual bool build_kkt_matrix(const hiopPDPerturbation& pdreg) + { + assert(false && "not yet implemented"); + return false; + } + + /** + * Solves the compressed linear system, part of the KKT Linear System interface + * + * Solves the system corresponding to directions for x, yc, and yd, namely + * [ H_BFGS + Dx Jc^T Jd^T ] [ dx] [ rx ] + * [ Jc 0 0 ] [dyc] = [ ryc] + * [ Jd 0 -Dd^{-1}] [dyd] [ ryd] + * + * This is done by forming and solving + * [ Jc*(H+Dx)^{-1}*Jc^T Jc*(H+Dx)^{-1}*Jd^T ] [dyc] = [ Jc(H+Dx)^{-1} rx - ryc ] + * [ Jd*(H+Dx)^{-1}*Jc^T Jd*(H+Dx)^{-1}*Jd^T + Dd^{-1}] [dyd] [ Jd(H+dx)^{-1} rx - ryd ] + * and then solving for dx from + * dx = - (H+Dx)^{-1}*(Jc^T*dyc+Jd^T*dyd - rx) + * + */ + virtual bool solveCompressed(hiopVector& rx, + hiopVector& ryc, + hiopVector& ryd, + hiopVector& dx, + hiopVector& dyc, + hiopVector& dyd); + + //LAPACK wrappers + int solve(hiopMatrixDense& M, hiopVector& rhs); + int solveWithRefin(hiopMatrixDense& M, hiopVector& rhs); +#ifdef HIOP_DEEPCHECKS + static double solveError(const hiopMatrixDense& M, const hiopVector& x, hiopVector& rhs); + double errorCompressedLinsys(const hiopVector& rx, + const hiopVector& ryc, + const hiopVector& ryd, + const hiopVector& dx, + const hiopVector& dyc, + const hiopVector& dyd); +protected: + /// @brief perform y=beta*y+alpha*H*x without the log barrier term from H + void HessianTimesVec_noLogBarrierTerm(double beta, hiopVector& y, double alpha, const hiopVector& x) + { + hiopHessianLowRank* HessLowR = dynamic_cast(Hess_); + assert(nullptr != HessLowR); + if(HessLowR) { + HessLowR->timesVec_noLogBarrierTerm(beta, y, alpha, x); + } + } +#endif + +private: + hiopNlpDenseConstraints* nlpD; + hiopHessianLowRank* HessLowRank; + + /// The kxk reduced matrix + hiopMatrixDense* N; +#ifdef HIOP_DEEPCHECKS + /// A copy of the above to compute the residual + hiopMatrixDense* Nmat; +#endif + //internal buffers: k is usually 2 x quasi-Newton memory; n is the size of primal variable vector + hiopMatrixDense* kxn_mat_; + hiopVector* k_vec1_; +}; +}; //end namespace + +#endif // HIOP_KKTLINSYSY_LOWRANK diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index 973c6d47..2209fcb3 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -56,6 +56,7 @@ #include "hiopAlgFilterIPM.hpp" #include "hiopKKTLinSys.hpp" +#include "KktLinSysLowRank.hpp" #include "hiopKKTLinSysDense.hpp" #include "hiopKKTLinSysMDS.hpp" #include "hiopKKTLinSysSparse.hpp" @@ -1047,7 +1048,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() theta_max = theta_max_fact_*fmax(1.0,resid->get_theta()); theta_min = theta_min_fact_*fmax(1.0,resid->get_theta()); - hiopKKTLinSysLowRank* kkt=new hiopKKTLinSysLowRank(nlp); + KktLinSysLowRank* kkt = new KktLinSysLowRank(nlp); assert(kkt != nullptr); // assign an Null pd_perturb, i.e., all the deltas = 0.0 diff --git a/src/Optimization/hiopIterate.hpp b/src/Optimization/hiopIterate.hpp index 8c4ce5dc..4039c452 100644 --- a/src/Optimization/hiopIterate.hpp +++ b/src/Optimization/hiopIterate.hpp @@ -165,7 +165,7 @@ class hiopIterate friend class hiopKKTLinSysCompressedXDYcYd; friend class hiopKKTLinSysDenseXYcYd; friend class hiopKKTLinSysDenseXDYcYd; - friend class hiopKKTLinSysLowRank; + friend class KktLinSysLowRank; friend class hiopHessianLowRank; friend class hiopKKTLinSysCompressedMDSXYcYd; friend class hiopHessianInvLowRank_obsolette; diff --git a/src/Optimization/hiopKKTLinSys.cpp b/src/Optimization/hiopKKTLinSys.cpp index d91a3d64..a2704909 100644 --- a/src/Optimization/hiopKKTLinSys.cpp +++ b/src/Optimization/hiopKKTLinSys.cpp @@ -1022,499 +1022,6 @@ errorCompressedLinsys(const hiopVector& rx, const hiopVector& rd, #endif -//////////////////////////////////////////////////////////////////////// -//////////////////////////////////////////////////////////////////////// -// hiopKKTLinSysLowRank -//////////////////////////////////////////////////////////////////////// -//////////////////////////////////////////////////////////////////////// - -hiopKKTLinSysLowRank::hiopKKTLinSysLowRank(hiopNlpFormulation* nlp) - : hiopKKTLinSysCompressedXYcYd(nlp) -{ - nlpD = dynamic_cast(nlp_); - - _kxn_mat = nlpD->alloc_multivector_primal(nlpD->m()); //!opt - assert("DEFAULT" == toupper(nlpD->options->GetString("mem_space"))); - N = LinearAlgebraFactory::create_matrix_dense(nlpD->options->GetString("mem_space"), - nlpD->m(), - nlpD->m()); -#ifdef HIOP_DEEPCHECKS - Nmat=N->alloc_clone(); -#endif - _k_vec1 = dynamic_cast(nlpD->alloc_dual_vec()); -} - -hiopKKTLinSysLowRank::~hiopKKTLinSysLowRank() -{ - if(N) delete N; -#ifdef HIOP_DEEPCHECKS - if(Nmat) delete Nmat; -#endif - if(_kxn_mat) delete _kxn_mat; - if(_k_vec1) delete _k_vec1; -} - -bool hiopKKTLinSysLowRank:: -update(const hiopIterate* iter, - const hiopVector* grad_f, - const hiopMatrixDense* Jac_c, const hiopMatrixDense* Jac_d, - hiopHessianLowRank* Hess) -{ - nlp_->runStats.tmSolverInternal.start(); - - iter_=iter; - grad_f_ = dynamic_cast(grad_f); - Jac_c_ = Jac_c; Jac_d_ = Jac_d; - //Hess = dynamic_cast(Hess_); - Hess_=HessLowRank=Hess; - - //compute the diagonals - //Dx=(Sxl)^{-1}Zl + (Sxu)^{-1}Zu - Dx_->setToZero(); - Dx_->axdzpy_w_pattern(1.0, *iter_->zl, *iter_->sxl, nlp_->get_ixl()); - Dx_->axdzpy_w_pattern(1.0, *iter_->zu, *iter_->sxu, nlp_->get_ixu()); - nlp_->log->write("Dx in KKT", *Dx_, hovMatrices); - - HessLowRank->updateLogBarrierDiagonal(*Dx_); - - //Dd=(Sdl)^{-1}Vu + (Sdu)^{-1}Vu - Dd_inv_->setToZero(); - Dd_inv_->axdzpy_w_pattern(1.0, *iter_->vl, *iter_->sdl, nlp_->get_idl()); - Dd_inv_->axdzpy_w_pattern(1.0, *iter_->vu, *iter_->sdu, nlp_->get_idu()); -#ifdef HIOP_DEEPCHECKS - assert(true==Dd_inv_->allPositive()); -#endif - Dd_->copyFrom(*Dd_inv_); - Dd_inv_->invert(); - - nlp_->runStats.tmSolverInternal.stop(); - - nlp_->log->write("Dd_inv in KKT", *Dd_inv_, hovMatrices); - return true; -} - - -/* Solves the system corresponding to directions for x, yc, and yd, namely - * [ H_BFGS + Dx Jc^T Jd^T ] [ dx] [ rx ] - * [ Jc 0 0 ] [dyc] = [ ryc ] - * [ Jd 0 -Dd^{-1}] [dyd] [ ryd ] - * - * This is done by forming and solving - * [ Jc*(H+Dx)^{-1}*Jc^T Jc*(H+Dx)^{-1}*Jd^T ] [dyc] = [ Jc(H+Dx)^{-1} rx - ryc ] - * [ Jd*(H+Dx)^{-1}*Jc^T Jd*(H+Dx)^{-1}*Jd^T + Dd^{-1}] [dyd] [ Jd(H+dx)^{-1} rx - ryd ] - * and then solving for dx from - * dx = - (H+Dx)^{-1}*(Jc^T*dyc+Jd^T*dyd - rx) - * - * Note that ops H+Dx are provided by hiopHessianLowRank - */ -bool hiopKKTLinSysLowRank:: -solveCompressed(hiopVector& rx, hiopVector& ryc, hiopVector& ryd, - hiopVector& dx, hiopVector& dyc, hiopVector& dyd) -{ -#ifdef HIOP_DEEPCHECKS - //some outputing - nlp_->log->write("KKT Low rank: solve compressed RHS", hovIteration); - nlp_->log->write(" rx: ", rx, hovIteration); nlp_->log->write(" ryc: ", ryc, hovIteration); nlp_->log->write(" ryd: ", ryd, hovIteration); - nlp_->log->write(" Jc: ", *Jac_c_, hovMatrices); - nlp_->log->write(" Jd: ", *Jac_d_, hovMatrices); - nlp_->log->write(" Dd_inv: ", *Dd_inv_, hovMatrices); - assert(Dd_inv_->isfinite_local() && "Something bad happened: nan or inf value"); -#endif - - hiopMatrixDense& J = *_kxn_mat; - const hiopMatrixDense* Jac_c_de = dynamic_cast(Jac_c_); assert(Jac_c_de); - const hiopMatrixDense* Jac_d_de = dynamic_cast(Jac_d_); assert(Jac_d_de); - J.copyRowsFrom(*Jac_c_de, nlp_->m_eq(), 0); //!opt - J.copyRowsFrom(*Jac_d_de, nlp_->m_ineq(), nlp_->m_eq());//!opt - - //N = J*(Hess\J') - //Hess->symmetricTimesMat(0.0, *N, 1.0, J); - HessLowRank->symMatTimesInverseTimesMatTrans(0.0, *N, 1.0, J); - - //subdiag of N += 1., Dd_inv - N->addSubDiagonal(1., nlp_->m_eq(), *Dd_inv_); -#ifdef HIOP_DEEPCHECKS - assert(J.isfinite()); - nlp_->log->write("solveCompressed: N is", *N, hovMatrices); - nlp_->log->write("solveCompressed: rx is", rx, hovMatrices); - nlp_->log->printf(hovLinAlgScalars, "inf norm of Dd_inv is %g\n", Dd_inv_->infnorm()); - N->assertSymmetry(1e-10); -#endif - - //compute the rhs of the lin sys involving N - // 1. first compute (H+Dx)^{-1} rx_tilde and store it temporarily in dx - HessLowRank->solve(rx, dx); -#ifdef HIOP_DEEPCHECKS - assert(rx.isfinite_local() && "Something bad happened: nan or inf value"); - assert(dx.isfinite_local() && "Something bad happened: nan or inf value"); -#endif - - // 2 . then rhs = [ Jc(H+Dx)^{-1}*rx - ryc ] - // [ Jd(H+dx)^{-1}*rx - ryd ] - hiopVector& rhs=*_k_vec1; - rhs.copyFromStarting(0, ryc); - rhs.copyFromStarting(nlp_->m_eq(), ryd); - J.timesVec(-1.0, rhs, 1.0, dx); - -#ifdef HIOP_DEEPCHECKS - nlp_->log->write("solveCompressed: dx sol is", dx, hovMatrices); - nlp_->log->write("solveCompressed: rhs for N is", rhs, hovMatrices); - Nmat->copyFrom(*N); - hiopVector* r=rhs.new_copy(); //save the rhs to check the norm of the residual -#endif - - // - //solve N * dyc_dyd = rhs - // - int ierr = solveWithRefin(*N,rhs); - //int ierr = solve(*N,rhs); - - hiopVector& dyc_dyd= rhs; - dyc_dyd.copyToStarting(0, dyc); - dyc_dyd.copyToStarting(nlp_->m_eq(), dyd); - - //now solve for dx = - (H+Dx)^{-1}*(Jc^T*dyc+Jd^T*dyd - rx) - //first rx = -(Jc^T*dyc+Jd^T*dyd - rx) - J.transTimesVec(1.0, rx, -1.0, dyc_dyd); - //then dx = (H+Dx)^{-1} rx - HessLowRank->solve(rx, dx); - -#ifdef HIOP_DEEPCHECKS - //some outputing - nlp_->log->write("KKT Low rank: solve compressed SOL", hovIteration); - nlp_->log->write(" dx: ", dx, hovIteration); nlp_->log->write(" dyc: ", dyc, hovIteration); nlp_->log->write(" dyd: ", dyd, hovIteration); - delete r; -#endif - - return ierr==0; -} - -int hiopKKTLinSysLowRank::solveWithRefin(hiopMatrixDense& M, hiopVector& rhs) -{ - // 1. Solve dposvx (solve + equilibrating + iterative refinement + forward and backward error estimates) - // 2. Check the residual norm - // 3. If residual norm is not small enough, then perform iterative refinement. This is because dposvx - // does not always provide a small enough residual since it stops (possibly without refinement) based on - // the forward and backward estimates - - int N=M.n(); - if(N<=0) return 0; - - hiopMatrixDense* Aref = M.new_copy(); - hiopVector* rhsref = rhs.new_copy(); - - char FACT='E'; - char UPLO='L'; - - int NRHS=1; - double* A=M.local_data(); - int LDA=N; - double* AF=new double[N*N]; - int LDAF=N; - char EQUED='N'; //it is an output if FACT='E' - double* S = new double[N]; - double* B = rhs.local_data(); - int LDB=N; - double* X = new double[N]; - int LDX = N; - double RCOND, FERR, BERR; - double* WORK = new double[3*N]; - int* IWORK = new int[N]; - int INFO; - - // - // 1. solve - // - DPOSVX(&FACT, &UPLO, &N, &NRHS, - A, &LDA, - AF, &LDAF, - &EQUED, - S, - B, &LDB, - X, &LDX, - &RCOND, &FERR, &BERR, - WORK, IWORK, - &INFO); - //printf("INFO ===== %d RCOND=%g FERR=%g BERR=%g EQUED=%c\n", INFO, RCOND, FERR, BERR, EQUED); - // - // 2. check residual - // - hiopVector* x = rhs.alloc_clone(); - hiopVector* dx = rhs.alloc_clone(); - hiopVector* resid = rhs.alloc_clone(); - int nIterRefin=0;double nrmResid; - int info; - const int MAX_ITER_REFIN=3; - while(true) { - x->copyFrom(X); - resid->copyFrom(*rhsref); - Aref->timesVec(1.0, *resid, -1.0, *x); - - nlp_->log->write("resid", *resid, hovLinAlgScalars); - - nrmResid= resid->infnorm(); - nlp_->log->printf(hovScalars, "hiopKKTLinSysLowRank::solveWithRefin iterrefin=%d residual norm=%g\n", nIterRefin, nrmResid); - - if(nrmResid<1e-8) break; - - if(nIterRefin>=MAX_ITER_REFIN) { - nlp_->log->write("N", *Aref, hovMatrices); - nlp_->log->write("sol", *x, hovMatrices); - nlp_->log->write("rhs", *rhsref, hovMatrices); - - nlp_->log->printf(hovWarning, "hiopKKTLinSysLowRank::solveWithRefin reduced residual to ONLY (inf-norm) %g after %d iterative refinements\n", nrmResid, nIterRefin); - break; - //assert(false && "too many refinements"); - } - if(0) { //iter refin based on symmetric indefinite factorization+solve - - - int _V_ipiv_vec[1000]; double _V_work_vec[1000]; int lwork=1000; - M.copyFrom(*Aref); - DSYTRF(&UPLO, &N, M.local_data(), &LDA, _V_ipiv_vec, _V_work_vec, &lwork, &info); - assert(info==0); - DSYTRS(&UPLO, &N, &NRHS, M.local_data(), &LDA, _V_ipiv_vec, resid->local_data(), &LDB, &info); - assert(info==0); - } else { //iter refin based on symmetric positive definite factorization+solve - M.copyFrom(*Aref); - //for(int i=0; i<4; i++) M.local_data()[i][i] +=1e-8; - DPOTRF(&UPLO, &N, M.local_data(), &LDA, &info); - if(info>0) - nlp_->log->printf(hovError, "hiopKKTLinSysLowRank::factorizeMat: dpotrf (Chol fact) " - "detected %d minor being indefinite.\n", info); - else - if(info<0) - nlp_->log->printf(hovError, "hiopKKTLinSysLowRank::factorizeMat: dpotrf returned " - "error %d\n", info); - - DPOTRS(&UPLO,&N, &NRHS, M.local_data(), &LDA, resid->local_data(), &LDA, &info); - if(info<0) - nlp_->log->printf(hovError, "hiopKKTLinSysLowRank::solveWithFactors: dpotrs returned " - "error %d\n", info); - } - - // //FACT='F'; EQUED='Y'; //reuse the factorization and the equilibration - // M.copyFrom(*Aref); - // A = M.local_buffer(); - // dposvx_(&FACT, &UPLO, &N, &NRHS, - // A, &LDA, - // AF, &LDAF, - // &EQUED, - // S, - // resid.local_data(), &LDB, - // X, &LDX, - // &RCOND, &FERR, &BERR, - // WORK, IWORK, - // &INFO); - // printf("INFO ===== %d RCOND=%g FERR=%g BERR=%g EQUED=%c\n", INFO, RCOND, FERR, BERR, EQUED); - - dx->copyFrom(*resid); - x->axpy(1., *dx); - - nIterRefin++; - } - - rhs.copyFrom(*x); - delete[] AF; - delete[] S; - delete[] X; - delete[] WORK; - delete[] IWORK; - delete Aref; - delete rhsref; - delete x; - delete dx; - delete resid; - -// #ifdef HIOP_DEEPCHECKS -// hiopVectorPar sol(rhs.get_size()); -// hiopVectorPar rhss(rhs.get_size()); -// sol.copyFrom(rhs); rhss.copyFrom(*r); -// double relErr=solveError(*Nmat, rhs, *r); -// if(relErr>1e-5) { -// nlp_->log->printf(hovWarning, "large rel. error (%g) in linear solver occured the Cholesky solve (hiopKKTLinSys)\n", relErr); - -// nlp_->log->write("matrix N=", *Nmat, hovError); -// nlp_->log->write("rhs", rhss, hovError); -// nlp_->log->write("sol", sol, hovError); - -// assert(false && "large error (%g) in linear solve (hiopKKTLinSys), equilibrating the matrix and/or iterative refinement are needed (see dposvx/x)"); -// } else -// if(relErr>1e-16) -// nlp_->log->printf(hovWarning, "considerable rel. error (%g) in linear solver occured the Cholesky solve (hiopKKTLinSys)\n", relErr); - -// nlp_->log->printf(hovLinAlgScalars, "hiopKKTLinSysLowRank::solveCompressed: Cholesky solve: relative error %g\n", relErr); -// delete r; -// #endif - return 0; -} - -int hiopKKTLinSysLowRank::solve(hiopMatrixDense& M, hiopVector& rhs) -{ - char FACT='E'; - char UPLO='L'; - int N=M.n(); - int NRHS=1; - double* A=M.local_data(); - int LDA=N; - double* AF=new double[N*N]; - int LDAF=N; - char EQUED='N'; //it is an output if FACT='E' - double* S = new double[N]; - double* B = rhs.local_data(); - int LDB=N; - double* X = new double[N]; - int LDX = N; - double RCOND, FERR, BERR; - double* WORK = new double[3*N]; - int* IWORK = new int[N]; - int INFO; - - DPOSVX(&FACT, &UPLO, &N, &NRHS, - A, &LDA, - AF, &LDAF, - &EQUED, - S, - B, &LDB, - X, &LDX, - &RCOND, &FERR, &BERR, - WORK, IWORK, - &INFO); - - rhs.copyFrom(S); - nlp_->log->write("Scaling S", rhs, hovSummary); - - //printf("INFO ===== %d RCOND=%g FERR=%g BERR=%g EQUED=%c\n", INFO, RCOND, FERR, BERR, EQUED); - - rhs.copyFrom(X); - delete [] AF; - delete [] S; - delete [] X; - delete [] WORK; - delete [] IWORK; - return 0; -} - -/* this code works fine but requires xblas -int hiopKKTLinSysLowRank::solveWithRefin(hiopMatrixDense& M, hiopVectorPar& rhs) -{ - char FACT='E'; - char UPLO='L'; - int N=M.n(); - int NRHS=1; - double* A=M.local_buffer(); - int LDA=N; - double* AF=new double[N*N]; - int LDAF=N; - char EQUED='N'; //it is an output if FACT='E' - double* S = new double[N]; - double* B = rhs.local_data(); - int LDB=N; - double* X = new double[N]; - int LDX = N; - double RCOND, BERR; - double RPVGRW; //Reciprocal pivot growth - int N_ERR_BNDS=3; - double* ERR_BNDS_NORM = new double[NRHS*N_ERR_BNDS]; - double* ERR_BNDS_COMP = new double[NRHS*N_ERR_BNDS]; - int NPARAMS=3; - double PARAMS[NPARAMS]; - PARAMS[0]=1.0; //Use the extra-precise refinement algorithm - PARAMS[1]=3.0; //Maximum number of residual computations allowed for refinement - PARAMS[2]=1.0; //attempt to find a solution with small componentwise - double* WORK = new double[4*N]; - int* IWORK = new int[N]; - int INFO; - - dposvxx_(&FACT, &UPLO, &N, &NRHS, - A, &LDA, - AF, &LDAF, - &EQUED, - S, - B, &LDB, - X, &LDX, - &RCOND, &RPVGRW, &BERR, - &N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, - &NPARAMS, PARAMS, - WORK, IWORK, - &INFO); - - //rhs.copyFrom(S); - //nlp_->log->write("Scaling S", rhs, hovSummary); - - //M.copyFrom(AF); - //nlp_->log->write("Factoriz ", M, hovSummary); - - printf("INFO ===== %d RCOND=%g RPVGRW=%g BERR=%g EQUED=%c\n", INFO, RCOND, RPVGRW, BERR, EQUED); - printf(" ERR_BNDS_NORM=%g %g %g ERR_BNDS_COMP=%g %g %g \n", ERR_BNDS_NORM[0], ERR_BNDS_NORM[1], ERR_BNDS_NORM[2], ERR_BNDS_COMP[0], ERR_BNDS_COMP[1], ERR_BNDS_COMP[2]); - printf(" PARAMS=%g %g %g \n", PARAMS[0], PARAMS[1], PARAMS[2]); - - - rhs.copyFrom(X); - delete [] AF; - delete [] S; - delete [] X; - delete [] ERR_BNDS_NORM; - delete [] ERR_BNDS_COMP; - delete [] WORK; - delete [] IWORK; - return 0; -} -*/ - -#ifdef HIOP_DEEPCHECKS - -double hiopKKTLinSysLowRank:: -errorCompressedLinsys(const hiopVector& rx, const hiopVector& ryc, const hiopVector& ryd, - const hiopVector& dx, const hiopVector& dyc, const hiopVector& dyd) -{ - nlp_->log->printf(hovLinAlgScalars, "hiopKKTLinSysLowRank::errorCompressedLinsys residuals norm:\n"); - - double derr=-1., aux; - hiopVector *RX=rx.new_copy(); - //RX=rx-H*dx-J'c*dyc-J'*dyd - HessLowRank->timesVec(1.0, *RX, -1.0, dx); - //RX->axzpy(-1.0,*Dx,dx); - Jac_c_->transTimesVec(1.0, *RX, -1.0, dyc); - Jac_d_->transTimesVec(1.0, *RX, -1.0, dyd); - aux=RX->twonorm(); - derr=fmax(derr,aux); - nlp_->log->printf(hovLinAlgScalars, " >>> rx=%g\n", aux); - // if(aux>1e-8) { - // nlp_->log->write("Low rank Hessian is:", *Hess, hovLinAlgScalars); - // } - delete RX; RX=NULL; - - hiopVector* RC=ryc.new_copy(); - Jac_c_->timesVec(1.0,*RC, -1.0,dx); - aux = RC->twonorm(); - derr=fmax(derr,aux); - nlp_->log->printf(hovLinAlgScalars, " >>> ryc=%g\n", aux); - delete RC; RC=NULL; - - hiopVector* RD=ryd.new_copy(); - Jac_d_->timesVec(1.0,*RD, -1.0, dx); - RD->axzpy(1.0, *Dd_inv_, dyd); - aux = RD->twonorm(); - derr=fmax(derr,aux); - nlp_->log->printf(hovLinAlgScalars, " >>> ryd=%g\n", aux); - delete RD; RD=NULL; - - return derr; -} - -double hiopKKTLinSysLowRank::solveError(const hiopMatrixDense& M, const hiopVector& x, hiopVector& rhs) -{ - double relError; - M.timesVec(1.0,rhs,-1.0,x); - double resnorm=rhs.infnorm(); - - relError=resnorm;// / (1+rhsnorm); - return relError; -} -#endif - - //////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////// // hiopKKTLinSysFull diff --git a/src/Optimization/hiopKKTLinSys.hpp b/src/Optimization/hiopKKTLinSys.hpp index 83334ba6..817106a8 100644 --- a/src/Optimization/hiopKKTLinSys.hpp +++ b/src/Optimization/hiopKKTLinSys.hpp @@ -382,82 +382,7 @@ class hiopKKTLinSysCompressedXDYcYd : public hiopKKTLinSysCompressed #endif }; -class hiopKKTLinSysLowRank : public hiopKKTLinSysCompressedXYcYd -{ -public: - hiopKKTLinSysLowRank(hiopNlpFormulation* nlp); - virtual ~hiopKKTLinSysLowRank(); - - bool update(const hiopIterate* iter, - const hiopVector* grad_f, - const hiopMatrix* Jac_c, const hiopMatrix* Jac_d, - hiopMatrix* Hess) - { - const hiopMatrixDense* Jac_c_ = dynamic_cast(Jac_c); - const hiopMatrixDense* Jac_d_ = dynamic_cast(Jac_d); - hiopHessianLowRank* Hess_ = dynamic_cast(Hess); - if(Jac_c_==NULL || Jac_d_==NULL || Hess_==NULL) { - assert(false); - return false; - } - return update(iter, grad_f_, Jac_c_, Jac_d_, Hess_); - } - virtual bool update(const hiopIterate* iter, - const hiopVector* grad_f, - const hiopMatrixDense* Jac_c, const hiopMatrixDense* Jac_d, - hiopHessianLowRank* Hess); - - virtual bool build_kkt_matrix(const hiopPDPerturbation& pdreg) - { - assert(false && "not yet implemented"); - return false; - } - - /* Solves the system corresponding to directions for x, yc, and yd, namely - * [ H_BFGS + Dx Jc^T Jd^T ] [ dx] [ rx_tilde ] - * [ Jc 0 0 ] [dyc] = [ ryc ] - * [ Jd 0 -Dd^{-1}] [dyd] [ ryd_tilde] - * - * This is done by forming and solving - * [ Jc*(H+Dx)^{-1}*Jc^T Jc*(H+Dx)^{-1}*Jd^T ] [dyc] = [ Jc(H+Dx)^{-1} rx - ryc ] - * [ Jd*(H+Dx)^{-1}*Jc^T Jd*(H+Dx)^{-1}*Jd^T + Dd^{-1}] [dyd] [ Jd(H+dx)^{-1} rx - ryd ] - * and then solving for dx from - * dx = - (H+Dx)^{-1}*(Jc^T*dyc+Jd^T*dyd - rx) - * - */ - virtual bool solveCompressed(hiopVector& rx, hiopVector& ryc, hiopVector& ryd, - hiopVector& dx, hiopVector& dyc, hiopVector& dyd); - - //LAPACK wrappers - int solve(hiopMatrixDense& M, hiopVector& rhs); - int solveWithRefin(hiopMatrixDense& M, hiopVector& rhs); -#ifdef HIOP_DEEPCHECKS - static double solveError(const hiopMatrixDense& M, const hiopVector& x, hiopVector& rhs); - double errorCompressedLinsys(const hiopVector& rx, const hiopVector& ryc, const hiopVector& ryd, - const hiopVector& dx, const hiopVector& dyc, const hiopVector& dyd); -protected: - //y=beta*y+alpha*H*x - void HessianTimesVec_noLogBarrierTerm(double beta, hiopVector& y, double alpha, const hiopVector& x) - { - hiopHessianLowRank* HessLowR = dynamic_cast(Hess_); - assert(NULL != HessLowR); - if(HessLowR) HessLowR->timesVec_noLogBarrierTerm(beta, y, alpha, x); - } -#endif - -private: - hiopNlpDenseConstraints* nlpD; - hiopHessianLowRank* HessLowRank; - - hiopMatrixDense* N; //the kxk reduced matrix -#ifdef HIOP_DEEPCHECKS - hiopMatrixDense* Nmat; //a copy of the above to compute the residual -#endif - //internal buffers - hiopMatrixDense* _kxn_mat; //!opt (work directly with the Jacobian) - hiopVector* _k_vec1; -}; /*