From 5e298d3083f3158d5d0ddc54bfbb69eb6b6a68fd Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Wed, 10 Apr 2024 21:36:55 -0700 Subject: [PATCH 01/20] Add FGMRES refinement in Primal Dual Full System solver --- src/Algorithm/IpPDFullSpaceSolver.cpp | 432 +++++++++++++++++++++++++- src/Algorithm/IpPDFullSpaceSolver.hpp | 36 +++ 2 files changed, 453 insertions(+), 15 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 195ab4dc0..ba1217c02 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -8,6 +8,7 @@ #include "IpDebug.hpp" #include +#include namespace Ipopt { @@ -138,6 +139,9 @@ bool PDFullSpaceSolver::Solve( DBG_ASSERT(!allow_inexact || !improve_solution); DBG_ASSERT(!improve_solution || beta == 0.); + return SolveGMRES( alpha, beta, rhs, res, allow_inexact, improve_solution ); + + // Timing of PDSystem solver starts here IpData().TimingStats().PDSystemSolverTotal().Start(); @@ -233,7 +237,7 @@ bool PDFullSpaceSolver::Solve( { SmartPtr resid = res.MakeNewIteratesVector(true); ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, - *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, alpha, beta, rhs, res, *resid); + *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, /*alpha*/1., /*beta*/-1., rhs, res, *resid); } break; } @@ -243,7 +247,7 @@ bool PDFullSpaceSolver::Solve( // ToDo don't to that after max refinement? ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, - *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, alpha, beta, rhs, res, *resid); + *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, /*alpha*/1., /*beta*/-1., rhs, res, *resid); Number residual_ratio = ComputeResidualRatio(rhs, res, *resid); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, @@ -263,7 +267,7 @@ bool PDFullSpaceSolver::Solve( ASSERT_EXCEPTION(solve_retval, INTERNAL_ABORT, "SolveOnce returns false during iterative refinement."); ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, - *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, alpha, beta, rhs, res, *resid); + *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, /*alpha*/1., /*beta*/-1., rhs, res, *resid); residual_ratio = ComputeResidualRatio(rhs, res, *resid); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, @@ -374,6 +378,252 @@ bool PDFullSpaceSolver::Solve( return true; } +bool PDFullSpaceSolver::SolveGMRES( + Number alpha, + Number beta, + const IteratesVector& rhs, + IteratesVector& res, + bool allow_inexact, + bool improve_solution /* = false */ +) +{ + DBG_START_METH("PDFullSpaceSolver::Solve", dbg_verbosity); + DBG_ASSERT(!allow_inexact || !improve_solution); + DBG_ASSERT(!improve_solution || beta == 0.); + + // Timing of PDSystem solver starts here + IpData().TimingStats().PDSystemSolverTotal().Start(); + + DBG_PRINT_VECTOR(2, "rhs_x", *rhs.x()); + DBG_PRINT_VECTOR(2, "rhs_s", *rhs.s()); + DBG_PRINT_VECTOR(2, "rhs_c", *rhs.y_c()); + DBG_PRINT_VECTOR(2, "rhs_d", *rhs.y_d()); + DBG_PRINT_VECTOR(2, "rhs_zL", *rhs.z_L()); + DBG_PRINT_VECTOR(2, "rhs_zU", *rhs.z_U()); + DBG_PRINT_VECTOR(2, "rhs_vL", *rhs.v_L()); + DBG_PRINT_VECTOR(2, "rhs_vU", *rhs.v_U()); + DBG_PRINT_VECTOR(2, "res_x in", *res.x()); + DBG_PRINT_VECTOR(2, "res_s in", *res.s()); + DBG_PRINT_VECTOR(2, "res_c in", *res.y_c()); + DBG_PRINT_VECTOR(2, "res_d in", *res.y_d()); + DBG_PRINT_VECTOR(2, "res_zL in", *res.z_L()); + DBG_PRINT_VECTOR(2, "res_zU in", *res.z_U()); + DBG_PRINT_VECTOR(2, "res_vL in", *res.v_L()); + DBG_PRINT_VECTOR(2, "res_vU in", *res.v_U()); + + // if beta is nonzero, keep a copy of the incoming values in res_ */ + SmartPtr copy_res; + if( beta != 0. ) + { + copy_res = res.MakeNewIteratesVectorCopy(); + } + + // Receive data about matrix +// SmartPtr x = IpData().curr()->x(); +// SmartPtr s = IpData().curr()->s(); + SmartPtr W = IpData().W(); + SmartPtr J_c = IpCq().curr_jac_c(); + SmartPtr J_d = IpCq().curr_jac_d(); + SmartPtr Px_L = IpNLP().Px_L(); + SmartPtr Px_U = IpNLP().Px_U(); + SmartPtr Pd_L = IpNLP().Pd_L(); + SmartPtr Pd_U = IpNLP().Pd_U(); + SmartPtr z_L = IpData().curr()->z_L(); + SmartPtr z_U = IpData().curr()->z_U(); + SmartPtr v_L = IpData().curr()->v_L(); + SmartPtr v_U = IpData().curr()->v_U(); + SmartPtr slack_x_L = IpCq().curr_slack_x_L(); + SmartPtr slack_x_U = IpCq().curr_slack_x_U(); + SmartPtr slack_s_L = IpCq().curr_slack_s_L(); + SmartPtr slack_s_U = IpCq().curr_slack_s_U(); + SmartPtr sigma_x = IpCq().curr_sigma_x(); + SmartPtr sigma_s = IpCq().curr_sigma_s(); + DBG_PRINT_VECTOR(2, "Sigma_x", *sigma_x); + DBG_PRINT_VECTOR(2, "Sigma_s", *sigma_s); + + bool done = false; + // The following flag is set to true, if we asked the linear + // solver to improve the quality of the solution in + // the next solve + bool resolve_with_better_quality = false; + // the following flag is set to true, if iterative refinement + // failed and we want to try if a modified system is able to + // remedy that problem by pretending the matrix is singular + bool pretend_singular = false; + bool pretend_singular_last_time = false; + + // Beginning of loop for solving the system (including all + // modifications for the linear system to ensure good solution + // quality) + while( !done ) + { + + // allow_inexact = true; + // std::cout << "allow_inexact= " << allow_inexact << std::endl; + // std::cout << "resolve_with_better_quality= " << resolve_with_better_quality << std::endl; + // std::cout << "pretend_singular= " << pretend_singular << std::endl; + + // if improve_solution is true, we are given already a solution + // from the calling function, so we can skip the first solve + bool solve_retval = true; + if( !improve_solution ) + { + solve_retval = SolveOnce(resolve_with_better_quality, pretend_singular, *W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, + *Pd_U, *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, 1., 0., + rhs, res); + resolve_with_better_quality = false; + pretend_singular = false; + } + improve_solution = false; + + if( !solve_retval ) + { + // If system seems not to be solvable, we return with false + // and let the calling routine deal with it. + IpData().TimingStats().PDSystemSolverTotal().End(); + return false; + } + + if( allow_inexact ) + { + // no safety checks required + if( Jnlst().ProduceOutput(J_MOREDETAILED, J_LINEAR_ALGEBRA) ) + { + SmartPtr resid = res.MakeNewIteratesVector(true); + ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, + *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, /*alpha*/1., /*beta*/-1., rhs, res, *resid); + } + break; + } + + // Get space for the residual + SmartPtr resid = res.MakeNewIteratesVector(true); + + // ToDo don't to that after max refinement? + ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, + *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, /*alpha*/1., /*beta*/-1., rhs, res, *resid); + + Number residual_ratio = ComputeResidualRatio(rhs, res, *resid); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "BEFORE GMRES residual_ratio = %e\n", residual_ratio); + + + // bool solve_retval = true; + bool gmres_conv = GMRES(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, + *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, + rhs, res, *resid, /*improve_solution*/true, resolve_with_better_quality, + pretend_singular, solve_retval); + if( !solve_retval ) + { + // If system seems not to be solvable, we return with false + // and let the calling routine deal with it. + IpData().TimingStats().PDSystemSolverTotal().End(); + return false; + } + improve_solution = false; + + residual_ratio = ComputeResidualRatio(rhs, res, *resid); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "residual_ratio = %e\n", residual_ratio); + + // Check if we have to give up on iterative refinement + if( !gmres_conv ) + { + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "GMRES refinement failed with residual_ratio = %e\n", residual_ratio); + // quit_refinement = true; + + // Pretend singularity only once - if it didn't help, we + // have to live with what we got so far + resolve_with_better_quality = false; + DBG_PRINT((1, "pretend_singular = %d\n", pretend_singular)); + if( !pretend_singular_last_time ) + { + // First try if we can ask the augmented system solver to + // improve the quality of the solution (only if that hasn't + // been done before for this linear system) + if( !augsys_improved_ ) + { + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "Asking augmented system solver to improve quality of its solutions.\n"); + augsys_improved_ = augSysSolver_->IncreaseQuality(); + if( augsys_improved_ ) + { + IpData().Append_info_string("q"); + resolve_with_better_quality = true; + } + else + { + // solver said it cannot improve quality, so let + // possibly conclude that the current modification is + // singular + pretend_singular = true; + } + } + else + { + // we had already asked the solver before to improve the + // quality of the solution, so let's now pretend that the + // modification is possibly singular + pretend_singular = true; + } + pretend_singular_last_time = pretend_singular; + if( pretend_singular ) + { + // let's only conclude that the current linear system + // including modifications is singular, if the residual is + // quite bad + if( residual_ratio < residual_ratio_singular_ ) + { + pretend_singular = false; + IpData().Append_info_string("S"); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "Just accept current solution.\n"); + } + else + { + IpData().Append_info_string("s"); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "Pretend that the current system (including modifications) is singular.\n"); + } + } + } + else + { + pretend_singular = false; + DBG_PRINT((1, "Resetting pretend_singular to false.\n")); + } + } + + done = !(resolve_with_better_quality) && !(pretend_singular); + + } // End of loop for solving the linear system (incl. modifications) + + // Finally let's assemble the res result vectors + if( alpha != 0. ) + { + res.Scal(alpha); + } + + if( beta != 0. ) + { + res.Axpy(beta, *copy_res); + } + + DBG_PRINT_VECTOR(2, "res_x", *res.x()); + DBG_PRINT_VECTOR(2, "res_s", *res.s()); + DBG_PRINT_VECTOR(2, "res_c", *res.y_c()); + DBG_PRINT_VECTOR(2, "res_d", *res.y_d()); + DBG_PRINT_VECTOR(2, "res_zL", *res.z_L()); + DBG_PRINT_VECTOR(2, "res_zU", *res.z_U()); + DBG_PRINT_VECTOR(2, "res_vL", *res.v_L()); + DBG_PRINT_VECTOR(2, "res_vU", *res.v_U()); + + IpData().TimingStats().PDSystemSolverTotal().End(); + + return true; +} + bool PDFullSpaceSolver::SolveOnce( bool resolve_with_better_quality, bool pretend_singular, @@ -681,8 +931,8 @@ void PDFullSpaceSolver::ComputeResiduals( const Vector& slack_s_U, const Vector& /*sigma_x*/, const Vector& /*sigma_s*/, - Number /*alpha*/, - Number /*beta*/, + Number alpha, + Number beta, const IteratesVector& rhs, const IteratesVector& res, IteratesVector& resid @@ -708,27 +958,27 @@ void PDFullSpaceSolver::ComputeResiduals( J_d.TransMultVector(1., *res.y_d(), 1., *resid.x_NonConst()); Px_L.MultVector(-1., *res.z_L(), 1., *resid.x_NonConst()); Px_U.MultVector(1., *res.z_U(), 1., *resid.x_NonConst()); - resid.x_NonConst()->AddTwoVectors(delta_x, *res.x(), -1., *rhs.x(), 1.); + resid.x_NonConst()->AddTwoVectors(alpha*delta_x, *res.x(), beta, *rhs.x(), alpha); // s Pd_U.MultVector(1., *res.v_U(), 0., *resid.s_NonConst()); Pd_L.MultVector(-1., *res.v_L(), 1., *resid.s_NonConst()); - resid.s_NonConst()->AddTwoVectors(-1., *res.y_d(), -1., *rhs.s(), 1.); + resid.s_NonConst()->AddTwoVectors(-alpha, *res.y_d(), beta, *rhs.s(), alpha); if( delta_s != 0. ) { - resid.s_NonConst()->Axpy(delta_s, *res.s()); + resid.s_NonConst()->Axpy(alpha*delta_s, *res.s()); } // c J_c.MultVector(1., *res.x(), 0., *resid.y_c_NonConst()); - resid.y_c_NonConst()->AddTwoVectors(-delta_c, *res.y_c(), -1., *rhs.y_c(), 1.); + resid.y_c_NonConst()->AddTwoVectors(-alpha*delta_c, *res.y_c(), beta, *rhs.y_c(), alpha); // d J_d.MultVector(1., *res.x(), 0., *resid.y_d_NonConst()); - resid.y_d_NonConst()->AddTwoVectors(-1., *res.s(), -1., *rhs.y_d(), 1.); + resid.y_d_NonConst()->AddTwoVectors(alpha, *res.s(), beta, *rhs.y_d(), alpha); if( delta_d != 0. ) { - resid.y_d_NonConst()->Axpy(-delta_d, *res.y_d()); + resid.y_d_NonConst()->Axpy(-alpha*delta_d, *res.y_d()); } // zL @@ -737,7 +987,7 @@ void PDFullSpaceSolver::ComputeResiduals( tmp = z_L.MakeNew(); Px_L.TransMultVector(1., *res.x(), 0., *tmp); tmp->ElementWiseMultiply(z_L); - resid.z_L_NonConst()->AddTwoVectors(1., *tmp, -1., *rhs.z_L(), 1.); + resid.z_L_NonConst()->AddTwoVectors(alpha, *tmp, beta, *rhs.z_L(), alpha); // zU resid.z_U_NonConst()->Copy(*res.z_U()); @@ -745,7 +995,7 @@ void PDFullSpaceSolver::ComputeResiduals( tmp = z_U.MakeNew(); Px_U.TransMultVector(1., *res.x(), 0., *tmp); tmp->ElementWiseMultiply(z_U); - resid.z_U_NonConst()->AddTwoVectors(-1., *tmp, -1., *rhs.z_U(), 1.); + resid.z_U_NonConst()->AddTwoVectors(-alpha, *tmp, beta, *rhs.z_U(), alpha); // vL resid.v_L_NonConst()->Copy(*res.v_L()); @@ -753,7 +1003,7 @@ void PDFullSpaceSolver::ComputeResiduals( tmp = v_L.MakeNew(); Pd_L.TransMultVector(1., *res.s(), 0., *tmp); tmp->ElementWiseMultiply(v_L); - resid.v_L_NonConst()->AddTwoVectors(1., *tmp, -1., *rhs.v_L(), 1.); + resid.v_L_NonConst()->AddTwoVectors(alpha, *tmp, beta, *rhs.v_L(), alpha); // vU resid.v_U_NonConst()->Copy(*res.v_U()); @@ -761,7 +1011,7 @@ void PDFullSpaceSolver::ComputeResiduals( tmp = v_U.MakeNew(); Pd_U.TransMultVector(1., *res.s(), 0., *tmp); tmp->ElementWiseMultiply(v_U); - resid.v_U_NonConst()->AddTwoVectors(-1., *tmp, -1., *rhs.v_U(), 1.); + resid.v_U_NonConst()->AddTwoVectors(-alpha, *tmp, beta, *rhs.v_U(), alpha); DBG_PRINT_VECTOR(2, "resid", resid); @@ -819,4 +1069,156 @@ Number PDFullSpaceSolver::ComputeResidualRatio( } } +bool PDFullSpaceSolver::GMRES( + const SymMatrix& W, + const Matrix& J_c, + const Matrix& J_d, + const Matrix& Px_L, + const Matrix& Px_U, + const Matrix& Pd_L, + const Matrix& Pd_U, + const Vector& z_L, + const Vector& z_U, + const Vector& v_L, + const Vector& v_U, + const Vector& slack_x_L, + const Vector& slack_x_U, + const Vector& slack_s_L, + const Vector& slack_s_U, + const Vector& sigma_x, + const Vector& sigma_s, + const IteratesVector& rhs, + IteratesVector& res, + IteratesVector& resid, + bool improve_solution, + bool resolve_with_better_quality, + bool pretend_singular, + bool& solve_retval +) +{ + DBG_START_METH("PDFullSpaceSolver::GMRES", dbg_verbosity); + + auto x = res.MakeNewIteratesVectorCopy(); + auto NullVec = res.MakeNewIteratesVectorCopy(); + NullVec->Set( 0. ); + + int restart = 500, totit = 0, ldH = restart + 1; + if (restart > max_refinement_steps_) + restart = max_refinement_steps_; + bool done = false; + Number tol = residual_ratio_max_, rho0 = 0., rho; + while ( !done ) + { + std::vector givens_c(restart), givens_s(restart), + b_(restart+1), H(restart*(restart+1)); + std::vector> V, Z; + V.emplace_back(rhs.MakeNewIteratesVector()); + V[0]->Set( 0. ); + if ( !improve_solution ) + { + x->Set(0.); + } + ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, + slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, *x, *V[0]); + improve_solution = true; // after restart, improve x from previous cycle + rho = V[0]->Nrm2(); + if( !totit ) + { + rho0 = rho; + } + V[0]->Scal( 1./rho ); + b_[0] = rho; + int nrit = restart - 1; + std::cout << "GMRES it. " << totit << "\tres = " << std::setw(12) + << rho << "\trel.res = " << std::setw(12) + << rho/rho0 << "\t restart!" + << " ||b||2= " << rhs.Nrm2() << " ||b||Inf= " << rhs.Amax() + << std::endl; + if( rho/rho0 < tol || rho < tol ) + { + done = true; + break; + } + for( Index it = 0; it < restart; it++ ) + { + totit++; + Z.emplace_back( rhs.MakeNewIteratesVector() ); + Z[it]->Set( 0. ); + solve_retval = SolveOnce(resolve_with_better_quality, pretend_singular, W, J_c, J_d, Px_L, Px_U, Pd_L, + Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U, sigma_x, sigma_s, 1., 0., + *V[it], *Z[it]); + if( !solve_retval ) + { + return false; + } + V.emplace_back( rhs.MakeNewIteratesVector() ); + V[it+1]->Set( 0. ); + ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, + slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., *NullVec, *Z[it], *V[it+1]); + // for( Index reps=0; reps<2; reps++ ) + { + for( Index k=0; k<=it; k++ ) + { + H[k+it*ldH] = V[it+1]->Dot( *V[k] ); + V[it+1]->Axpy( -H[k+it*ldH], *V[k] ); + } + } + H[it+1+it*ldH] = V[it+1]->Nrm2(); + V[it+1]->Scal( 1. / H[it+1+it*ldH] ); + for( Index k = 1; k < it+1; k++ ) + { + auto gamma = givens_c[k-1]*H[k-1+it*ldH] + givens_s[k-1]*H[k+it*ldH]; + H[k+it*ldH] = -givens_s[k-1]*H[k-1+it*ldH] + givens_c[k-1]*H[k+it*ldH]; + H[k-1+it*ldH] = gamma; + } + auto delta = std::sqrt( H[it+it*ldH]*H[it+it*ldH] + H[it+1+it*ldH]*H[it+1+it*ldH] ); + givens_c[it] = H[it+it*ldH] / delta; + givens_s[it] = H[it+1+it*ldH] / delta; + H[it+it*ldH] = givens_c[it]*H[it+it*ldH] + givens_s[it]*H[it+1+it*ldH]; + b_[it+1] = -givens_s[it]*b_[it]; + b_[it] = givens_c[it]*b_[it]; + rho = std::abs( b_[it+1] ); + std::cout << "GMRES it. " << totit << "\tres = " << std::setw(12) + << rho << "\trel.res = " << std::setw(12) + << rho/rho0 << std::endl; + if( (rho < tol) || (rho/rho0 < tol) || (totit >= max_refinement_steps_) ) + { + done = true; + nrit = it; + break; + } + } + for( Index k = nrit; k >= 0; k-- ) + { + for ( Index i = k+1; i <= nrit; i++ ) + b_[k] -= H[k+i*ldH]*b_[i]; + b_[k] /= H[k+k*ldH]; + x->Axpy( b_[k], *Z[k] ); + } + } + ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, + slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, *x, resid); + std::cout << "GMRES it. " << totit + << " rho= " << rho << " rho/rho0= " << rho/rho0 + << " ||r||2= " << resid.Nrm2() + << " ||b||2= " << rhs.Nrm2() + << " ||r||2/||b||2= " << resid.Nrm2() / rhs.Nrm2() + << " ||r||Inf= " << resid.Amax() + << " ||b||Inf= " << rhs.Amax() + << " ||r||Inf/||b||Inf= " << resid.Amax() / rhs.Amax() << std::endl; + + if ( totit < max_refinement_steps_ && + ( resid.Nrm2() / rhs.Nrm2() < tol ) ) + { + std::cout << "Accepting GMRES solution" << std::endl; + res.Copy( *x ); + return true; + } + else + { + return false; + } + //return totit < max_refinement_steps_; +} + } // namespace Ipopt diff --git a/src/Algorithm/IpPDFullSpaceSolver.hpp b/src/Algorithm/IpPDFullSpaceSolver.hpp index 6e117821a..ec1112de1 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.hpp +++ b/src/Algorithm/IpPDFullSpaceSolver.hpp @@ -224,6 +224,42 @@ class PDFullSpaceSolver: public PDSystemSolver Vector& X ); ///@} + + bool SolveGMRES( + Number alpha, + Number beta, + const IteratesVector& rhs, + IteratesVector& res, + bool allow_inexact = false, + bool improve_solution = false + ); + + bool GMRES( + const SymMatrix& W, + const Matrix& J_c, + const Matrix& J_d, + const Matrix& Px_L, + const Matrix& Px_U, + const Matrix& Pd_L, + const Matrix& Pd_U, + const Vector& z_L, + const Vector& z_U, + const Vector& v_L, + const Vector& v_U, + const Vector& slack_x_L, + const Vector& slack_x_U, + const Vector& slack_s_L, + const Vector& slack_s_U, + const Vector& sigma_x, + const Vector& sigma_s, + const IteratesVector& rhs, + IteratesVector& res, + IteratesVector& resid, + bool improve_solution, + bool resolve_with_better_quality, + bool pretend_singular, + bool& solve_retval + ); }; } // namespace Ipopt From ad8b0e51ff997fffb2adba3bacebeddba8bf43e9 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Wed, 10 Apr 2024 21:48:12 -0700 Subject: [PATCH 02/20] small sign fix --- src/Algorithm/IpPDFullSpaceSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index ba1217c02..0635cfca6 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -975,7 +975,7 @@ void PDFullSpaceSolver::ComputeResiduals( // d J_d.MultVector(1., *res.x(), 0., *resid.y_d_NonConst()); - resid.y_d_NonConst()->AddTwoVectors(alpha, *res.s(), beta, *rhs.y_d(), alpha); + resid.y_d_NonConst()->AddTwoVectors(-alpha, *res.s(), beta, *rhs.y_d(), alpha); if( delta_d != 0. ) { resid.y_d_NonConst()->Axpy(-alpha*delta_d, *res.y_d()); From e6c2c0fd4e9d4bf25dcf66c5d2821aa84eff17b9 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Wed, 10 Apr 2024 22:13:13 -0700 Subject: [PATCH 03/20] Few fixes in FGMRES for PDFull --- src/Algorithm/IpPDFullSpaceSolver.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 0635cfca6..72acd5188 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1102,9 +1102,10 @@ bool PDFullSpaceSolver::GMRES( auto NullVec = res.MakeNewIteratesVectorCopy(); NullVec->Set( 0. ); - int restart = 500, totit = 0, ldH = restart + 1; + Index restart = 10; if (restart > max_refinement_steps_) restart = max_refinement_steps_; + Index totit = 0, ldH = restart + 1; bool done = false; Number tol = residual_ratio_max_, rho0 = 0., rho; while ( !done ) @@ -1155,7 +1156,7 @@ bool PDFullSpaceSolver::GMRES( V[it+1]->Set( 0. ); ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., *NullVec, *Z[it], *V[it+1]); - // for( Index reps=0; reps<2; reps++ ) + for( Index reps=0; reps<2; reps++ ) { for( Index k=0; k<=it; k++ ) { @@ -1210,7 +1211,6 @@ bool PDFullSpaceSolver::GMRES( if ( totit < max_refinement_steps_ && ( resid.Nrm2() / rhs.Nrm2() < tol ) ) { - std::cout << "Accepting GMRES solution" << std::endl; res.Copy( *x ); return true; } @@ -1218,7 +1218,6 @@ bool PDFullSpaceSolver::GMRES( { return false; } - //return totit < max_refinement_steps_; } } // namespace Ipopt From a947da1c96fe824d54c21115490502e62dc3ae3d Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Wed, 10 Apr 2024 22:27:00 -0700 Subject: [PATCH 04/20] Do MGS twice --- src/Algorithm/IpPDFullSpaceSolver.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 72acd5188..bb8a6b1bb 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1156,7 +1156,6 @@ bool PDFullSpaceSolver::GMRES( V[it+1]->Set( 0. ); ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., *NullVec, *Z[it], *V[it+1]); - for( Index reps=0; reps<2; reps++ ) { for( Index k=0; k<=it; k++ ) { @@ -1164,6 +1163,14 @@ bool PDFullSpaceSolver::GMRES( V[it+1]->Axpy( -H[k+it*ldH], *V[k] ); } } + { // second time for numerical stability, ToDo optional? + for( Index k=0; k<=it; k++ ) + { + auto tmp = V[it+1]->Dot( *V[k] ); + H[k+it*ldH] += tmp; + V[it+1]->Axpy( -tmp, *V[k] ); + } + } H[it+1+it*ldH] = V[it+1]->Nrm2(); V[it+1]->Scal( 1. / H[it+1+it*ldH] ); for( Index k = 1; k < it+1; k++ ) From d2e34b86d7d9cc3bffdb5fe7b407de957a5617c1 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Tue, 23 Apr 2024 15:04:47 -0700 Subject: [PATCH 05/20] Update FGMRES for primal-dual system. Add code for computing matrix Inf norm --- src/Algorithm/IpPDFullSpaceSolver.cpp | 245 ++++++++++++++++----- src/Algorithm/IpPDFullSpaceSolver.hpp | 26 ++- src/LinAlg/IpCompoundMatrix.cpp | 112 ++++++++++ src/LinAlg/IpCompoundMatrix.hpp | 10 + src/LinAlg/IpCompoundSymMatrix.cpp | 51 +++++ src/LinAlg/IpCompoundSymMatrix.hpp | 5 + src/LinAlg/IpDenseGenMatrix.cpp | 16 ++ src/LinAlg/IpDenseGenMatrix.hpp | 10 + src/LinAlg/IpDenseSymMatrix.cpp | 9 + src/LinAlg/IpDenseSymMatrix.hpp | 5 + src/LinAlg/IpDiagMatrix.cpp | 19 ++ src/LinAlg/IpDiagMatrix.hpp | 5 + src/LinAlg/IpExpandedMultiVectorMatrix.cpp | 16 ++ src/LinAlg/IpExpandedMultiVectorMatrix.hpp | 12 +- src/LinAlg/IpExpansionMatrix.cpp | 32 +++ src/LinAlg/IpExpansionMatrix.hpp | 10 + src/LinAlg/IpIdentityMatrix.cpp | 15 ++ src/LinAlg/IpIdentityMatrix.hpp | 5 + src/LinAlg/IpLowRankUpdateSymMatrix.cpp | 16 ++ src/LinAlg/IpLowRankUpdateSymMatrix.hpp | 12 +- src/LinAlg/IpMatrix.hpp | 56 +++++ src/LinAlg/IpMultiVectorMatrix.cpp | 16 ++ src/LinAlg/IpMultiVectorMatrix.hpp | 10 + src/LinAlg/IpScaledMatrix.cpp | 79 +++++++ src/LinAlg/IpScaledMatrix.hpp | 10 + src/LinAlg/IpSumMatrix.cpp | 16 ++ src/LinAlg/IpSumMatrix.hpp | 10 + src/LinAlg/IpSumSymMatrix.cpp | 24 ++ src/LinAlg/IpSumSymMatrix.hpp | 10 + src/LinAlg/IpSymMatrix.hpp | 13 ++ src/LinAlg/IpSymScaledMatrix.cpp | 8 + src/LinAlg/IpSymScaledMatrix.hpp | 5 + src/LinAlg/IpTransposeMatrix.hpp | 16 ++ src/LinAlg/IpZeroMatrix.hpp | 12 + src/LinAlg/IpZeroSymMatrix.hpp | 12 + src/LinAlg/TMatrices/IpGenTMatrix.cpp | 54 +++++ src/LinAlg/TMatrices/IpGenTMatrix.hpp | 10 + src/LinAlg/TMatrices/IpSymTMatrix.cpp | 36 +++ src/LinAlg/TMatrices/IpSymTMatrix.hpp | 5 + 39 files changed, 974 insertions(+), 59 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index bb8a6b1bb..42df1a830 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -78,6 +78,14 @@ void PDFullSpaceSolver::RegisterOptions( "If the improvement of the residual test ratio made by one iterative refinement step is not better than this factor, " "iterative refinement is aborted.", true); + roptions->AddStringOption2( + "gmres_refinement", + "Whether to use flexible GMRES for refinement instead of iterative refinement." + "The min_refinement_steps, max_refinement_steps and residual_ratio_max" + "options are also used with GMRES.", + "yes", + "yes", "use FGMRES instead of iterative refinement", + "no", "use original IPOPT approach, with iterative refinement"); roptions->AddLowerBoundedNumberOption( "neg_curv_test_tol", "Tolerance for heuristic to ignore wrong inertia.", @@ -112,6 +120,7 @@ bool PDFullSpaceSolver::InitializeImpl( ASSERT_EXCEPTION(residual_ratio_singular_ >= residual_ratio_max_, OPTION_INVALID, "Option \"residual_ratio_singular\": This value must be not smaller than residual_ratio_max."); options.GetNumericValue("residual_improvement_factor", residual_improvement_factor_, prefix); + options.GetBoolValue("gmres_refinement", gmres_refinement_, prefix); options.GetNumericValue("neg_curv_test_tol", neg_curv_test_tol_, prefix); options.GetBoolValue("neg_curv_test_reg", neg_curv_test_reg_, prefix); @@ -139,8 +148,10 @@ bool PDFullSpaceSolver::Solve( DBG_ASSERT(!allow_inexact || !improve_solution); DBG_ASSERT(!improve_solution || beta == 0.); - return SolveGMRES( alpha, beta, rhs, res, allow_inexact, improve_solution ); - + if( gmres_refinement_ ) + { + return SolveGMRES( alpha, beta, rhs, res, allow_inexact, improve_solution ); + } // Timing of PDSystem solver starts here IpData().TimingStats().PDSystemSolverTotal().Start(); @@ -245,11 +256,14 @@ bool PDFullSpaceSolver::Solve( // Get space for the residual SmartPtr resid = res.MakeNewIteratesVector(true); + Number A_nrm_inf = NrmInf(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, + *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U, *resid); + // ToDo don't to that after max refinement? ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, /*alpha*/1., /*beta*/-1., rhs, res, *resid); - Number residual_ratio = ComputeResidualRatio(rhs, res, *resid); + Number residual_ratio = ComputeResidualRatio(rhs, res, *resid, A_nrm_inf); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "residual_ratio = %e\n", residual_ratio); Number residual_ratio_old = residual_ratio; @@ -269,7 +283,7 @@ bool PDFullSpaceSolver::Solve( ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, /*alpha*/1., /*beta*/-1., rhs, res, *resid); - residual_ratio = ComputeResidualRatio(rhs, res, *resid); + residual_ratio = ComputeResidualRatio(rhs, res, *resid, A_nrm_inf); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "residual_ratio = %e\n", residual_ratio); @@ -458,11 +472,6 @@ bool PDFullSpaceSolver::SolveGMRES( while( !done ) { - // allow_inexact = true; - // std::cout << "allow_inexact= " << allow_inexact << std::endl; - // std::cout << "resolve_with_better_quality= " << resolve_with_better_quality << std::endl; - // std::cout << "pretend_singular= " << pretend_singular << std::endl; - // if improve_solution is true, we are given already a solution // from the calling function, so we can skip the first solve bool solve_retval = true; @@ -499,11 +508,14 @@ bool PDFullSpaceSolver::SolveGMRES( // Get space for the residual SmartPtr resid = res.MakeNewIteratesVector(true); + Number A_nrm_inf = NrmInf(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, + *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U, *resid); + // ToDo don't to that after max refinement? ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, /*alpha*/1., /*beta*/-1., rhs, res, *resid); - Number residual_ratio = ComputeResidualRatio(rhs, res, *resid); + Number residual_ratio = ComputeResidualRatio(rhs, res, *resid, A_nrm_inf); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "BEFORE GMRES residual_ratio = %e\n", residual_ratio); @@ -522,7 +534,7 @@ bool PDFullSpaceSolver::SolveGMRES( } improve_solution = false; - residual_ratio = ComputeResidualRatio(rhs, res, *resid); + residual_ratio = ComputeResidualRatio(rhs, res, *resid, A_nrm_inf); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "residual_ratio = %e\n", residual_ratio); @@ -1045,7 +1057,8 @@ void PDFullSpaceSolver::ComputeResiduals( Number PDFullSpaceSolver::ComputeResidualRatio( const IteratesVector& rhs, const IteratesVector& res, - const IteratesVector& resid + const IteratesVector& resid, + Number AInfNrm ) { DBG_START_METH("PDFullSpaceSolver::ComputeResidualRatio", dbg_verbosity); @@ -1066,9 +1079,91 @@ Number PDFullSpaceSolver::ComputeResidualRatio( // safeguard to use against incredibly large solution vectors Number max_cond = 1e6; return nrm_resid / (Min(nrm_res, max_cond * nrm_rhs) + nrm_rhs); + // return nrm_resid / (AInfNrm * nrm_res /*+ nrm_rhs*/ ); } } +Number PDFullSpaceSolver::NrmInf( + const SymMatrix& W, + const Matrix& J_c, + const Matrix& J_d, + const Matrix& Px_L, + const Matrix& Px_U, + const Matrix& Pd_L, + const Matrix& Pd_U, + const Vector& z_L, + const Vector& z_U, + const Vector& v_L, + const Vector& v_U, + const Vector& slack_x_L, + const Vector& slack_x_U, + const Vector& slack_s_L, + const Vector& slack_s_U, + IteratesVector& tmp +) +{ + DBG_START_METH("PDFullSpaceSolver::NrmInf", dbg_verbosity); + + // Get the current sizes of the perturbation factors + Number delta_x; + Number delta_s; + Number delta_c; + Number delta_d; + perturbHandler_->CurrentPerturbation(delta_x, delta_s, delta_c, delta_d); + + tmp.Set( 0. ); + W.ComputeRowA1(*tmp.x_NonConst(), false); + J_c.ComputeColA1(*tmp.x_NonConst(), true); + J_d.ComputeColA1(*tmp.x_NonConst(), true); + Px_L.ComputeRowA1(*tmp.x_NonConst(), true); + Px_U.ComputeRowA1(*tmp.x_NonConst(), true); + Number A_inf = tmp.x_NonConst()->Amax() + std::abs(delta_x); + + Pd_L.ComputeRowA1(*tmp.s_NonConst(), false); + Pd_U.ComputeRowA1(*tmp.s_NonConst(), true); + A_inf = std::max(A_inf, tmp.s_NonConst()->Amax() + 1. + delta_s); + + J_c.ComputeRowA1(*tmp.y_c_NonConst(), false); + A_inf = std::max(A_inf, tmp.y_c_NonConst()->Amax() + delta_c); + + J_d.ComputeRowA1(*tmp.y_d_NonConst(), false); + A_inf = std::max(A_inf, tmp.y_d_NonConst()->Amax() + 1. + delta_d); + + Px_L.ComputeColA1(*tmp.z_L_NonConst(), false); + tmp.z_L_NonConst()->ElementWiseMultiply(z_L); + tmp.z_L_NonConst()->ElementWiseAbs(); + auto tmp_slack_x_L = slack_x_L.MakeNewCopy(); + tmp_slack_x_L->ElementWiseAbs(); + tmp.z_L_NonConst()->Axpy(1., *tmp_slack_x_L); + A_inf = std::max(A_inf, tmp.z_L_NonConst()->Amax()); + + Px_U.ComputeColA1(*tmp.z_U_NonConst(), false); + tmp.z_U_NonConst()->ElementWiseMultiply(z_U); + tmp.z_U_NonConst()->ElementWiseAbs(); + auto tmp_slack_x_U = slack_x_U.MakeNewCopy(); + tmp_slack_x_U->ElementWiseAbs(); + tmp.z_U_NonConst()->Axpy(1., *tmp_slack_x_U); + A_inf = std::max(A_inf, tmp.z_U_NonConst()->Amax()); + + Pd_L.ComputeColA1(*tmp.v_L_NonConst(), false); + tmp.v_L_NonConst()->ElementWiseMultiply(v_L); + tmp.v_L_NonConst()->ElementWiseAbs(); + auto tmp_slack_s_L = slack_s_L.MakeNewCopy(); + tmp_slack_s_L->ElementWiseAbs(); + tmp.v_L_NonConst()->Axpy(1., *tmp_slack_s_L); + A_inf = std::max(A_inf, tmp.v_L_NonConst()->Amax()); + + Pd_U.ComputeColA1(*tmp.v_U_NonConst(), false); + tmp.v_U_NonConst()->ElementWiseMultiply(v_U); + tmp.v_U_NonConst()->ElementWiseAbs(); + auto tmp_slack_s_U = slack_s_U.MakeNewCopy(); + tmp_slack_s_U->ElementWiseAbs(); + tmp.v_U_NonConst()->Axpy(1., *tmp_slack_s_U); + A_inf = std::max(A_inf, tmp.v_U_NonConst()->Amax()); + + return A_inf; +} + bool PDFullSpaceSolver::GMRES( const SymMatrix& W, const Matrix& J_c, @@ -1098,16 +1193,23 @@ bool PDFullSpaceSolver::GMRES( { DBG_START_METH("PDFullSpaceSolver::GMRES", dbg_verbosity); + const auto default_precision{std::cout.precision()}; + std::cout << std::setprecision(3); + auto x = res.MakeNewIteratesVectorCopy(); auto NullVec = res.MakeNewIteratesVectorCopy(); NullVec->Set( 0. ); - Index restart = 10; - if (restart > max_refinement_steps_) - restart = max_refinement_steps_; + // Index restart = 10; + // if (restart > max_refinement_steps_) + // restart = max_refinement_steps_; + Index restart = max_refinement_steps_; Index totit = 0, ldH = restart + 1; bool done = false; Number tol = residual_ratio_max_, rho0 = 0., rho; + Number b_nrm_2 = rhs.Nrm2(), b_nrm_max = rhs.Amax(); + Number A_nrm_inf = NrmInf(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, + slack_x_L, slack_x_U, slack_s_L, slack_s_U, resid); while ( !done ) { std::vector givens_c(restart), givens_s(restart), @@ -1127,19 +1229,41 @@ bool PDFullSpaceSolver::GMRES( { rho0 = rho; } - V[0]->Scal( 1./rho ); - b_[0] = rho; - int nrit = restart - 1; - std::cout << "GMRES it. " << totit << "\tres = " << std::setw(12) - << rho << "\trel.res = " << std::setw(12) - << rho/rho0 << "\t restart!" - << " ||b||2= " << rhs.Nrm2() << " ||b||Inf= " << rhs.Amax() + Number resid_nrm_max = resid.Amax(); + // norm-wise relative backward error, in 2 or Inf norm + // Number NRBE_2 = rho / ( A_nrm_inf * x->Nrm2() + b_nrm_2 ); + Number NRBE_inf = resid_nrm_max / ( A_nrm_inf * x->Amax() /*+ b_nrm_max*/ ); + // Number NRBE_inf = ComputeResidualRatio(rhs, *x, *V[0], A_nrm_inf); + // auto abs_resid = V[0]->MakeNewIteratesVectorCopy(); + // auto abs_rhs = rhs.MakeNewIteratesVectorCopy(); + // auto abs_x = x->MakeNewIteratesVectorCopy(); + // abs_resid->ElementWiseAbs(); + // abs_rhs->ElementWiseAbs(); + // abs_x->ElementWiseAbs(); + // abs_rhs->Axpy(A_nrm_inf, *abs_x); + // abs_rhs->AddScalar(Number(1.e-16)); + // abs_resid->ElementWiseDivide(*abs_rhs); + // Number CRBE = abs_resid->Amax(); + std::cout << "GMRES it. " << totit + // << " NRBE_2= " << std::setw(8) << NRBE_2 + << " NRBE_inf= " << std::setw(8) << NRBE_inf + // << " CRBE= " << std::setw(8) << CRBE + << " ||r||2= " << std::setw(8) << rho + << " ||r||2/||b||2= " << std::setw(8) << rho / b_nrm_2 + << " ||r||Inf= " << std::setw(8) << resid_nrm_max + << " ||r||Inf/||b||Inf= " << std::setw(8) << resid_nrm_max / b_nrm_max + << " ||A||Inf= " << std::setw(8) << A_nrm_inf + // << " ||b||Inf= " << std::setw(8) << b_nrm_max << std::endl; - if( rho/rho0 < tol || rho < tol ) + if( NRBE_inf < tol && totit >= min_refinement_steps_) { + resid.Copy( *V[0] ); done = true; break; } + V[0]->Scal( 1./rho ); + b_[0] = rho; + // int nrit = restart - 1; for( Index it = 0; it < restart; it++ ) { totit++; @@ -1186,45 +1310,54 @@ bool PDFullSpaceSolver::GMRES( b_[it+1] = -givens_s[it]*b_[it]; b_[it] = givens_c[it]*b_[it]; rho = std::abs( b_[it+1] ); - std::cout << "GMRES it. " << totit << "\tres = " << std::setw(12) - << rho << "\trel.res = " << std::setw(12) - << rho/rho0 << std::endl; - if( (rho < tol) || (rho/rho0 < tol) || (totit >= max_refinement_steps_) ) + + std::vector y = b_; + for( Index k = it; k >= 0; k-- ) + { + for ( Index i = k+1; i <= it; i++ ) + y[k] -= H[k+i*ldH] * y[i]; + y[k] /= H[k+k*ldH]; + } + x->Axpy( y[it], *Z[it] ); + + ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, + slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, *x, resid); + Number resid_nrm_2 = resid.Nrm2(); + resid_nrm_max = resid.Amax(); + // NRBE_2 = resid_nrm_2 / ( A_nrm_inf * x->Nrm2() + b_nrm_2 ); + NRBE_inf = resid_nrm_max / ( A_nrm_inf * x->Amax() /*+ b_nrm_max*/ ); + // NRBE_inf = ComputeResidualRatio(rhs, *x, resid, A_nrm_inf); + // abs_resid->Copy( resid ); + // abs_x->Copy( *x ); + // abs_resid->ElementWiseAbs(); + // abs_x->ElementWiseAbs(); + // abs_rhs->Axpy(A_nrm_inf, *abs_x); + // abs_rhs->AddScalar(Number(1.e-16)); + // abs_resid->ElementWiseDivide(*abs_rhs); + // Number CRBE = abs_resid->Amax(); + std::cout << "GMRES it. " << totit + // << " NRBE_2= " << std::setw(8) << NRBE_2 + << " NRBE_inf= " << std::setw(8) << NRBE_inf + // << " CRBE= " << std::setw(8) << CRBE + << " ||r||2= " << std::setw(8) << resid_nrm_2 + << " ||r||2/||b||2= " << std::setw(8) << resid_nrm_2 / b_nrm_2 + << " ||r||Inf= " << std::setw(8) << resid_nrm_max + << " ||r||Inf/||b||Inf= " << std::setw(8) << resid_nrm_max / b_nrm_max + // << " rho= " << std::setw(8) << rho + << std::endl; + if( (NRBE_inf < tol && totit >= min_refinement_steps_) || + (totit >= max_refinement_steps_) ) { done = true; - nrit = it; break; } } - for( Index k = nrit; k >= 0; k-- ) - { - for ( Index i = k+1; i <= nrit; i++ ) - b_[k] -= H[k+i*ldH]*b_[i]; - b_[k] /= H[k+k*ldH]; - x->Axpy( b_[k], *Z[k] ); - } - } - ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, - slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, *x, resid); - std::cout << "GMRES it. " << totit - << " rho= " << rho << " rho/rho0= " << rho/rho0 - << " ||r||2= " << resid.Nrm2() - << " ||b||2= " << rhs.Nrm2() - << " ||r||2/||b||2= " << resid.Nrm2() / rhs.Nrm2() - << " ||r||Inf= " << resid.Amax() - << " ||b||Inf= " << rhs.Amax() - << " ||r||Inf/||b||Inf= " << resid.Amax() / rhs.Amax() << std::endl; - - if ( totit < max_refinement_steps_ && - ( resid.Nrm2() / rhs.Nrm2() < tol ) ) - { - res.Copy( *x ); - return true; - } - else - { - return false; } + + std::cout << std::setprecision(default_precision); + + res.Copy( *x ); + return ( totit < max_refinement_steps_ ); } } // namespace Ipopt diff --git a/src/Algorithm/IpPDFullSpaceSolver.hpp b/src/Algorithm/IpPDFullSpaceSolver.hpp index ec1112de1..a94a73e62 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.hpp +++ b/src/Algorithm/IpPDFullSpaceSolver.hpp @@ -132,6 +132,9 @@ class PDFullSpaceSolver: public PDSystemSolver */ Number residual_improvement_factor_; + /** Use (flexible) GMRES refinement instead of iterative refinement */ + bool gmres_refinement_; + /** Tolerance for heuristic to ignore wrong inertia */ Number neg_curv_test_tol_; @@ -208,7 +211,8 @@ class PDFullSpaceSolver: public PDSystemSolver Number ComputeResidualRatio( const IteratesVector& rhs, const IteratesVector& res, - const IteratesVector& resid + const IteratesVector& resid, + Number AInfNrm ); /** @name Auxiliary functions */ @@ -260,6 +264,26 @@ class PDFullSpaceSolver: public PDSystemSolver bool pretend_singular, bool& solve_retval ); + + Number NrmInf( + const SymMatrix& W, + const Matrix& J_c, + const Matrix& J_d, + const Matrix& Px_L, + const Matrix& Px_U, + const Matrix& Pd_L, + const Matrix& Pd_U, + const Vector& z_L, + const Vector& z_U, + const Vector& v_L, + const Vector& v_U, + const Vector& slack_x_L, + const Vector& slack_x_U, + const Vector& slack_s_L, + const Vector& slack_s_U, + IteratesVector& tmp + ); + }; } // namespace Ipopt diff --git a/src/LinAlg/IpCompoundMatrix.cpp b/src/LinAlg/IpCompoundMatrix.cpp index c994c78b8..d47f55570 100644 --- a/src/LinAlg/IpCompoundMatrix.cpp +++ b/src/LinAlg/IpCompoundMatrix.cpp @@ -715,6 +715,118 @@ void CompoundMatrix::ComputeColAMaxImpl( } } +void CompoundMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool /*init*/ +) const +{ + if( !matrices_valid_ ) + { + matrices_valid_ = MatricesValid(); + } + DBG_ASSERT(matrices_valid_); + + // The vector is assumed to be compound Vectors as well except if + // there is only one component + CompoundVector* comp_vec = dynamic_cast(&rows_norms); + +#ifndef ALLOW_NESTED + // A few sanity checks + if (comp_vec) + { + DBG_ASSERT(NComps_Rows() == comp_vec->NComps()); + } + else + { + DBG_ASSERT(NComps_Rows() == 1); + } +#endif + if( comp_vec ) + { + if( NComps_Rows() != comp_vec->NComps() ) + { + comp_vec = NULL; + } + } + + for( Index jcol = 0; jcol < NComps_Cols(); jcol++ ) + { + for( Index irow = 0; irow < NComps_Rows(); irow++ ) + { + if( ConstComp(irow, jcol) ) + { + SmartPtr vec_i; + if( comp_vec ) + { + vec_i = comp_vec->GetCompNonConst(irow); + } + else + { + vec_i = &rows_norms; + } + DBG_ASSERT(IsValid(vec_i)); + ConstComp(irow, jcol)->ComputeRowA1(*vec_i, false); + } + } + } +} + +void CompoundMatrix::ComputeColA1Impl( + Vector& cols_norms, + bool /*init*/ +) const +{ + if( !matrices_valid_ ) + { + matrices_valid_ = MatricesValid(); + } + DBG_ASSERT(matrices_valid_); + + // The vector is assumed to be compound Vectors as well except if + // there is only one component + CompoundVector* comp_vec = dynamic_cast(&cols_norms); + +#ifndef ALLOW_NESTED + // A few sanity checks + if (comp_vec) + { + DBG_ASSERT(NComps_Cols() == comp_vec->NComps()); + } + else + { + DBG_ASSERT(NComps_Cols() == 1); + } +#endif + if( comp_vec ) + { + if( NComps_Cols() != comp_vec->NComps() ) + { + comp_vec = NULL; + } + } + + for( Index irow = 0; irow < NComps_Rows(); irow++ ) + { + for( Index jcol = 0; jcol < NComps_Cols(); jcol++ ) + { + if( ConstComp(irow, jcol) ) + { + SmartPtr vec_j; + if( comp_vec ) + { + vec_j = comp_vec->GetCompNonConst(jcol); + } + else + { + vec_j = &cols_norms; + } + DBG_ASSERT(IsValid(vec_j)); + ConstComp(irow, jcol)->ComputeColA1(*vec_j, false); + } + } + } +} + void CompoundMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpCompoundMatrix.hpp b/src/LinAlg/IpCompoundMatrix.hpp index 821942559..67fe820e5 100644 --- a/src/LinAlg/IpCompoundMatrix.hpp +++ b/src/LinAlg/IpCompoundMatrix.hpp @@ -155,6 +155,16 @@ class IPOPTLIB_EXPORT CompoundMatrix: public Matrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpCompoundSymMatrix.cpp b/src/LinAlg/IpCompoundSymMatrix.cpp index 0c190a43c..9de8f5e56 100644 --- a/src/LinAlg/IpCompoundSymMatrix.cpp +++ b/src/LinAlg/IpCompoundSymMatrix.cpp @@ -237,6 +237,57 @@ void CompoundSymMatrix::ComputeRowAMaxImpl( } } +void CompoundSymMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool /*init*/ +) const +{ + if( !matrices_valid_ ) + { + matrices_valid_ = MatricesValid(); + } + DBG_ASSERT(matrices_valid_); + + // The vector is assumed to be compound Vectors as well except if + // there is only one component + CompoundVector* comp_vec = dynamic_cast(&rows_norms); + + // A few sanity checks + if( comp_vec ) + { + DBG_ASSERT(NComps_Dim() == comp_vec->NComps()); + } + else + { + DBG_ASSERT(NComps_Dim() == 1); + } + + for( Index jcol = 0; jcol < NComps_Dim(); jcol++ ) + { + for( Index irow = 0; irow < NComps_Dim(); irow++ ) + { + SmartPtr vec_i; + if( comp_vec ) + { + vec_i = comp_vec->GetCompNonConst(irow); + } + else + { + vec_i = &rows_norms; + } + DBG_ASSERT(IsValid(vec_i)); + if( jcol <= irow && ConstComp(irow, jcol) ) + { + ConstComp(irow, jcol)->ComputeRowA1(*vec_i, false); + } + else if( jcol > irow && ConstComp(jcol, irow) ) + { + ConstComp(jcol, irow)->ComputeRowA1(*vec_i, false); + } + } + } +} + void CompoundSymMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpCompoundSymMatrix.hpp b/src/LinAlg/IpCompoundSymMatrix.hpp index c1571b03b..b5183b39c 100644 --- a/src/LinAlg/IpCompoundSymMatrix.hpp +++ b/src/LinAlg/IpCompoundSymMatrix.hpp @@ -117,6 +117,11 @@ class IPOPTLIB_EXPORT CompoundSymMatrix: public SymMatrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpDenseGenMatrix.cpp b/src/LinAlg/IpDenseGenMatrix.cpp index 481f65bf5..d4d3d0522 100644 --- a/src/LinAlg/IpDenseGenMatrix.cpp +++ b/src/LinAlg/IpDenseGenMatrix.cpp @@ -415,6 +415,22 @@ void DenseGenMatrix::ComputeColAMaxImpl( } } +void DenseGenMatrix::ComputeRowA1Impl( + Vector& /*rows_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "DenseGenMatrix::ComputeRowA1Impl not implemented"); +} + +void DenseGenMatrix::ComputeColA1Impl( + Vector& /*cols_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "DenseGenMatrix::ComputeColA1Impl not implemented"); +} + bool DenseGenMatrix::HasValidNumbersImpl() const { DBG_ASSERT(initialized_); diff --git a/src/LinAlg/IpDenseGenMatrix.hpp b/src/LinAlg/IpDenseGenMatrix.hpp index b3bebb10e..47cd8190a 100644 --- a/src/LinAlg/IpDenseGenMatrix.hpp +++ b/src/LinAlg/IpDenseGenMatrix.hpp @@ -217,6 +217,16 @@ class IPOPTLIB_EXPORT DenseGenMatrix: public Matrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpDenseSymMatrix.cpp b/src/LinAlg/IpDenseSymMatrix.cpp index 0761d7845..3ad15284a 100644 --- a/src/LinAlg/IpDenseSymMatrix.cpp +++ b/src/LinAlg/IpDenseSymMatrix.cpp @@ -254,6 +254,15 @@ void DenseSymMatrix::ComputeRowAMaxImpl( } } +void DenseSymMatrix::ComputeRowA1Impl( + Vector& /*rows_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "DenseSymMatrix::ComputeRowA1Impl not implemented"); +} + + void DenseSymMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpDenseSymMatrix.hpp b/src/LinAlg/IpDenseSymMatrix.hpp index 506fcd1de..1f3a8931a 100644 --- a/src/LinAlg/IpDenseSymMatrix.hpp +++ b/src/LinAlg/IpDenseSymMatrix.hpp @@ -133,6 +133,11 @@ class IPOPTLIB_EXPORT DenseSymMatrix: public SymMatrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpDiagMatrix.cpp b/src/LinAlg/IpDiagMatrix.cpp index db7d29fee..51a3efea8 100644 --- a/src/LinAlg/IpDiagMatrix.cpp +++ b/src/LinAlg/IpDiagMatrix.cpp @@ -71,6 +71,25 @@ void DiagMatrix::ComputeRowAMaxImpl( } } +void DiagMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool init +) const +{ + DBG_ASSERT(IsValid(diag_)); + if( init ) + { + rows_norms.Copy(*diag_); + rows_norms.ElementWiseAbs(); + } + else + { + SmartPtr v = diag_->MakeNewCopy(); + v->ElementWiseAbs(); + rows_norms.Axpy(Number(1.), *v); + } +} + void DiagMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpDiagMatrix.hpp b/src/LinAlg/IpDiagMatrix.hpp index f12f84fd5..672b45713 100644 --- a/src/LinAlg/IpDiagMatrix.hpp +++ b/src/LinAlg/IpDiagMatrix.hpp @@ -63,6 +63,11 @@ class IPOPTLIB_EXPORT DiagMatrix: public SymMatrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpExpandedMultiVectorMatrix.cpp b/src/LinAlg/IpExpandedMultiVectorMatrix.cpp index 9b45e88c7..5bec4cfa6 100644 --- a/src/LinAlg/IpExpandedMultiVectorMatrix.cpp +++ b/src/LinAlg/IpExpandedMultiVectorMatrix.cpp @@ -192,6 +192,22 @@ void ExpandedMultiVectorMatrix::ComputeColAMaxImpl( THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ExpandedMultiVectorMatrix::ComputeColAMaxImpl not implemented"); } +void ExpandedMultiVectorMatrix::ComputeRowA1Impl( + Vector& /*rows_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ExpandedMultiVectorMatrix::ComputeRowA1Impl not implemented"); +} + +void ExpandedMultiVectorMatrix::ComputeColA1Impl( + Vector& /*cols_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ExpandedMultiVectorMatrix::ComputeColA1Impl not implemented"); +} + void ExpandedMultiVectorMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpExpandedMultiVectorMatrix.hpp b/src/LinAlg/IpExpandedMultiVectorMatrix.hpp index 1bf709249..a35e7d8a7 100644 --- a/src/LinAlg/IpExpandedMultiVectorMatrix.hpp +++ b/src/LinAlg/IpExpandedMultiVectorMatrix.hpp @@ -101,7 +101,17 @@ class ExpandedMultiVectorMatrix: public Matrix bool init ) const; - virtual void PrintImpl( + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, EJournalCategory category, diff --git a/src/LinAlg/IpExpansionMatrix.cpp b/src/LinAlg/IpExpansionMatrix.cpp index 47edcb460..27b286d77 100644 --- a/src/LinAlg/IpExpansionMatrix.cpp +++ b/src/LinAlg/IpExpansionMatrix.cpp @@ -405,6 +405,38 @@ void ExpansionMatrix::ComputeColAMaxImpl( } } +void ExpansionMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool /*init*/ +) const +{ + DenseVector* dense_vec = static_cast(&rows_norms); + DBG_ASSERT(dynamic_cast(&rows_norms)); + Number* vec_vals = dense_vec->Values(); + + const Index* exp_pos = ExpandedPosIndices(); + + for( Index i = 0; i < NCols(); i++ ) + { + vec_vals[exp_pos[i]] += Number(1.); + } +} + +void ExpansionMatrix::ComputeColA1Impl( + Vector& cols_norms, + bool init +) const +{ + if( init ) + { + cols_norms.Set(1.); + } + else + { + cols_norms.AddScalar(Number(1.)); + } +} + void ExpansionMatrix::PrintImplOffset( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpExpansionMatrix.hpp b/src/LinAlg/IpExpansionMatrix.hpp index e23355db1..1208d2365 100644 --- a/src/LinAlg/IpExpansionMatrix.hpp +++ b/src/LinAlg/IpExpansionMatrix.hpp @@ -103,6 +103,16 @@ class IPOPTLIB_EXPORT ExpansionMatrix: public Matrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpIdentityMatrix.cpp b/src/LinAlg/IpIdentityMatrix.cpp index 10ff48e33..e84c533a1 100644 --- a/src/LinAlg/IpIdentityMatrix.cpp +++ b/src/LinAlg/IpIdentityMatrix.cpp @@ -62,6 +62,21 @@ void IdentityMatrix::ComputeRowAMaxImpl( } } +void IdentityMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool init +) const +{ + if( init ) + { + rows_norms.Set(1.); + } + else + { + rows_norms.AddScalar(Number(1.)); + } +} + void IdentityMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpIdentityMatrix.hpp b/src/LinAlg/IpIdentityMatrix.hpp index 6f0e4296c..efbfad228 100644 --- a/src/LinAlg/IpIdentityMatrix.hpp +++ b/src/LinAlg/IpIdentityMatrix.hpp @@ -71,6 +71,11 @@ class IPOPTLIB_EXPORT IdentityMatrix: public SymMatrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpLowRankUpdateSymMatrix.cpp b/src/LinAlg/IpLowRankUpdateSymMatrix.cpp index a4128da16..5e97e0e62 100644 --- a/src/LinAlg/IpLowRankUpdateSymMatrix.cpp +++ b/src/LinAlg/IpLowRankUpdateSymMatrix.cpp @@ -170,6 +170,22 @@ void LowRankUpdateSymMatrix::ComputeColAMaxImpl( THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "LowRankUpdateSymMatrix::ComputeColAMaxImpl not implemented"); } +void LowRankUpdateSymMatrix::ComputeRowA1Impl( + Vector& /*rows_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "LowRankUpdateSymMatrix::ComputeRowA1Impl not implemented"); +} + +void LowRankUpdateSymMatrix::ComputeColA1Impl( + Vector& /*cols_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "LowRankUpdateSymMatrix::ComputeColA1Impl not implemented"); +} + void LowRankUpdateSymMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpLowRankUpdateSymMatrix.hpp b/src/LinAlg/IpLowRankUpdateSymMatrix.hpp index c88c8cd66..74f75182c 100644 --- a/src/LinAlg/IpLowRankUpdateSymMatrix.hpp +++ b/src/LinAlg/IpLowRankUpdateSymMatrix.hpp @@ -122,7 +122,17 @@ class LowRankUpdateSymMatrix: public SymMatrix bool init ) const; - virtual void PrintImpl( + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, EJournalCategory category, diff --git a/src/LinAlg/IpMatrix.hpp b/src/LinAlg/IpMatrix.hpp index 6b2555004..2e12735bf 100644 --- a/src/LinAlg/IpMatrix.hpp +++ b/src/LinAlg/IpMatrix.hpp @@ -163,6 +163,44 @@ class IPOPTLIB_EXPORT Matrix: public TaggedObject } ///@} + /** Compute the 1-norm of the rows in the matrix. + * + * The result is stored in rows_norms. + * The vector is assumed to be initialized if init is false. + */ + void ComputeRowA1( + Vector& rows_norms, + bool init = true + ) const + { + DBG_ASSERT(NRows() == rows_norms.Dim()); + if( init ) + { + rows_norms.Set(0.); + } + ComputeRowA1Impl(rows_norms, init); + } + ///@} + + /** Compute the 1-norm of the columns in the matrix. + * + * The result is stored in cols_norms. + * The vector is assumed to be initialized if init is false. + */ + void ComputeColA1( + Vector& cols_norms, + bool init = true + ) const + { + DBG_ASSERT(NCols() == cols_norms.Dim()); + if( init ) + { + cols_norms.Set(0.); + } + ComputeColA1Impl(cols_norms, init); + } + ///@} + /** Print detailed information about the matrix. * * @attention Do not overload. Overload PrintImpl instead. @@ -270,6 +308,24 @@ class IPOPTLIB_EXPORT Matrix: public TaggedObject bool init ) const = 0; + /** Compute the 1-norm of the rows in the matrix. + * + * The result is stored in rows_norms. + * The vector is assumed to be initialized if init is false. */ + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const = 0; + + /** Compute the 1-norm of the cols in the matrix. + * + * The result is stored in cols_norms. + * The vector is assumed to be initialized if init is false. */ + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const = 0; + /** Print detailed information about the matrix. */ virtual void PrintImpl( const Journalist& jnlst, diff --git a/src/LinAlg/IpMultiVectorMatrix.cpp b/src/LinAlg/IpMultiVectorMatrix.cpp index 71c764e2c..473964a02 100644 --- a/src/LinAlg/IpMultiVectorMatrix.cpp +++ b/src/LinAlg/IpMultiVectorMatrix.cpp @@ -301,6 +301,22 @@ void MultiVectorMatrix::ComputeColAMaxImpl( THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "MultiVectorMatrix::ComputeColAMaxImpl not implemented"); } +void MultiVectorMatrix::ComputeRowA1Impl( + Vector& /*rows_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "MultiVectorMatrix::ComputeRowA1Impl not implemented"); +} + +void MultiVectorMatrix::ComputeColA1Impl( + Vector& /*cols_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "MultiVectorMatrix::ComputeColA1Impl not implemented"); +} + void MultiVectorMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpMultiVectorMatrix.hpp b/src/LinAlg/IpMultiVectorMatrix.hpp index 319ac66d9..fd37cfd37 100644 --- a/src/LinAlg/IpMultiVectorMatrix.hpp +++ b/src/LinAlg/IpMultiVectorMatrix.hpp @@ -165,6 +165,16 @@ class IPOPTLIB_EXPORT MultiVectorMatrix: public Matrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpScaledMatrix.cpp b/src/LinAlg/IpScaledMatrix.cpp index c75330184..9e4c24adb 100644 --- a/src/LinAlg/IpScaledMatrix.cpp +++ b/src/LinAlg/IpScaledMatrix.cpp @@ -5,6 +5,11 @@ // Authors: Carl Laird, Andreas Waechter IBM 2004-08-13 #include "IpScaledMatrix.hpp" +#include "IpGenTMatrix.hpp" +#include "IpDenseVector.hpp" + +#include +#include namespace Ipopt { @@ -117,6 +122,80 @@ void ScaledMatrix::ComputeColAMaxImpl( THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ScaledMatrix::ComputeColAMaxImpl not implemented"); } +void ScaledMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool init +) const +{ + DBG_ASSERT(IsValid(matrix_)); + + if( IsValid(owner_space_->ColumnScaling()) ) + { + const GenTMatrix* mt = dynamic_cast(GetRawPtr(matrix_)); + if ( mt ) + { + const DenseVector* Dc = dynamic_cast(GetRawPtr(owner_space_->ColumnScaling())); + DenseVector* rn = dynamic_cast(&rows_norms); + DBG_ASSERT(Dc); + DBG_ASSERT(rn); + for ( Index i = 0; i < mt->Nonzeros(); i++ ) + { + rn->Values()[mt->Irows()[i] - 1] += std::abs(mt->Values()[i] * Dc->Values()[mt->Jcols()[i] - 1]); + } + } + else + { + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ScaledMatrix::ComputeRowA1Impl not implemented"); + } + } + else + { + matrix_->ComputeRowA1(rows_norms, init); + } + if( IsValid(owner_space_->RowScaling()) ) + { + rows_norms.ElementWiseMultiply(*owner_space_->RowScaling()); + rows_norms.ElementWiseAbs(); + } +} + +void ScaledMatrix::ComputeColA1Impl( + Vector& cols_norms, + bool init +) const +{ + DBG_ASSERT(IsValid(matrix_)); + + if( IsValid(owner_space_->RowScaling()) ) + { + const GenTMatrix* mt = dynamic_cast(GetRawPtr(matrix_)); + if ( mt ) + { + const DenseVector* Dr = dynamic_cast(GetRawPtr(owner_space_->RowScaling())); + DenseVector* cn = dynamic_cast(&cols_norms); + DBG_ASSERT(Dr); + DBG_ASSERT(cn); + for ( Index i = 0; i < mt->Nonzeros(); i++ ) + { + cn->Values()[mt->Jcols()[i] - 1] += std::abs(mt->Values()[i] * Dr->Values()[mt->Irows()[i] - 1]); + } + } + else + { + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ScaledMatrix::ComputeColA1Impl not implemented"); + } + } + else + { + matrix_->ComputeColA1(cols_norms, init); + } + if( IsValid(owner_space_->ColumnScaling()) ) + { + cols_norms.ElementWiseMultiply(*owner_space_->ColumnScaling()); + cols_norms.ElementWiseAbs(); + } +} + void ScaledMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpScaledMatrix.hpp b/src/LinAlg/IpScaledMatrix.hpp index f75b3a6d7..f9c695b0e 100644 --- a/src/LinAlg/IpScaledMatrix.hpp +++ b/src/LinAlg/IpScaledMatrix.hpp @@ -89,6 +89,16 @@ class IPOPTLIB_EXPORT ScaledMatrix: public Matrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpSumMatrix.cpp b/src/LinAlg/IpSumMatrix.cpp index c50746f89..f42cd5e0e 100644 --- a/src/LinAlg/IpSumMatrix.cpp +++ b/src/LinAlg/IpSumMatrix.cpp @@ -136,6 +136,22 @@ void SumMatrix::ComputeColAMaxImpl( THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SumMatrix::ComputeColAMaxImpl not implemented"); } +void SumMatrix::ComputeRowA1Impl( + Vector& /*rows_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SumMatrix::ComputeRowA1Impl not implemented"); +} + +void SumMatrix::ComputeColA1Impl( + Vector& /*cols_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SumMatrix::ComputeColA1Impl not implemented"); +} + void SumMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpSumMatrix.hpp b/src/LinAlg/IpSumMatrix.hpp index 48a16d879..0e602d091 100644 --- a/src/LinAlg/IpSumMatrix.hpp +++ b/src/LinAlg/IpSumMatrix.hpp @@ -83,6 +83,16 @@ class SumMatrix: public Matrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpSumSymMatrix.cpp b/src/LinAlg/IpSumSymMatrix.cpp index dfb2fe6d4..51ba0fd51 100644 --- a/src/LinAlg/IpSumSymMatrix.cpp +++ b/src/LinAlg/IpSumSymMatrix.cpp @@ -108,6 +108,30 @@ void SumSymMatrix::ComputeColAMaxImpl( THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SumSymMatrix::ComputeColAMaxImpl not implemented"); } +void SumSymMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool /*init*/ +) const +{ + for( Index iterm = 0; iterm < NTerms(); iterm++ ) + { + DBG_ASSERT(IsValid(matrices_[iterm])); + matrices_[iterm]->ComputeRowA1(rows_norms, false); + } +} + +void SumSymMatrix::ComputeColA1Impl( + Vector& cols_norms, + bool /*init*/ +) const +{ + for( Index iterm = 0; iterm < NTerms(); iterm++ ) + { + DBG_ASSERT(IsValid(matrices_[iterm])); + matrices_[iterm]->ComputeColA1(cols_norms, false); + } +} + void SumSymMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpSumSymMatrix.hpp b/src/LinAlg/IpSumSymMatrix.hpp index 4f1108308..a021cd163 100644 --- a/src/LinAlg/IpSumSymMatrix.hpp +++ b/src/LinAlg/IpSumSymMatrix.hpp @@ -81,6 +81,16 @@ class IPOPTLIB_EXPORT SumSymMatrix: public SymMatrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpSymMatrix.hpp b/src/LinAlg/IpSymMatrix.hpp index 294b6c1db..800c30e8e 100644 --- a/src/LinAlg/IpSymMatrix.hpp +++ b/src/LinAlg/IpSymMatrix.hpp @@ -74,6 +74,19 @@ class IPOPTLIB_EXPORT SymMatrix: public Matrix } ///@} + /** Implementation of ComputeColA1Impl, which calls ComputeRowA1Impl. + * + * Since the matrix is symmetric, the row and column 1 norms are identical. + */ + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const + { + ComputeRowA1Impl(cols_norms, init); + } + ///@} + private: /** Copy of the owner space ptr as a SymMatrixSpace instead * of a MatrixSpace diff --git a/src/LinAlg/IpSymScaledMatrix.cpp b/src/LinAlg/IpSymScaledMatrix.cpp index 6f967612a..ce32f1b7a 100644 --- a/src/LinAlg/IpSymScaledMatrix.cpp +++ b/src/LinAlg/IpSymScaledMatrix.cpp @@ -70,6 +70,14 @@ void SymScaledMatrix::ComputeRowAMaxImpl( THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SymScaledMatrix::ComputeRowAMaxImpl not implemented"); } +void SymScaledMatrix::ComputeRowA1Impl( + Vector& /*rows_norms*/, + bool /*init*/ +) const +{ + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SymScaledMatrix::ComputeRowA1Impl not implemented"); +} + void SymScaledMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpSymScaledMatrix.hpp b/src/LinAlg/IpSymScaledMatrix.hpp index c63c14b26..fa259ae3b 100644 --- a/src/LinAlg/IpSymScaledMatrix.hpp +++ b/src/LinAlg/IpSymScaledMatrix.hpp @@ -74,6 +74,11 @@ class IPOPTLIB_EXPORT SymScaledMatrix: public SymMatrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpTransposeMatrix.hpp b/src/LinAlg/IpTransposeMatrix.hpp index b2efaed72..6bbdda16e 100644 --- a/src/LinAlg/IpTransposeMatrix.hpp +++ b/src/LinAlg/IpTransposeMatrix.hpp @@ -85,6 +85,22 @@ class TransposeMatrix: public Matrix orig_matrix_->ComputeRowAMax(rows_norms, init); } + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const { + DBG_ASSERT(IsValid(orig_matrix_)); + orig_matrix_->ComputeColA1(rows_norms, init); + } + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const { + DBG_ASSERT(IsValid(orig_matrix_)); + orig_matrix_->ComputeRowA1(cols_norms, init); + } + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpZeroMatrix.hpp b/src/LinAlg/IpZeroMatrix.hpp index 48e9b727e..38513d5ae 100644 --- a/src/LinAlg/IpZeroMatrix.hpp +++ b/src/LinAlg/IpZeroMatrix.hpp @@ -58,6 +58,18 @@ class ZeroMatrix: public Matrix ) const { } + virtual void ComputeRowA1Impl( + Vector& /*rows_norms*/, + bool /*init*/ + ) const + { } + + virtual void ComputeColA1Impl( + Vector& /*cols_norms*/, + bool /*init*/ + ) const + { } + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/IpZeroSymMatrix.hpp b/src/LinAlg/IpZeroSymMatrix.hpp index 4d1e67f00..f9fdb598e 100644 --- a/src/LinAlg/IpZeroSymMatrix.hpp +++ b/src/LinAlg/IpZeroSymMatrix.hpp @@ -57,6 +57,18 @@ class IPOPTLIB_EXPORT ZeroSymMatrix: public SymMatrix ) const { } + virtual void ComputeRowA1Impl( + Vector& /*rows_norms*/, + bool /*init*/ + ) const + { } + + virtual void ComputeColA1Impl( + Vector& /*cols_norms*/, + bool /*init*/ + ) const + { } + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/TMatrices/IpGenTMatrix.cpp b/src/LinAlg/TMatrices/IpGenTMatrix.cpp index bcfca90e0..cb2b8f4d5 100644 --- a/src/LinAlg/TMatrices/IpGenTMatrix.cpp +++ b/src/LinAlg/TMatrices/IpGenTMatrix.cpp @@ -238,6 +238,60 @@ void GenTMatrix::ComputeColAMaxImpl( } } +void GenTMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool /*init*/ +) const +{ + DBG_ASSERT(initialized_); + + if( NRows() == 0 ) + { + return; + } + + DenseVector* dense_vec = static_cast(&rows_norms); + DBG_ASSERT(dynamic_cast(&rows_norms)); + + const Index* irows = Irows(); + const Number* val = values_; + Number* vec_vals = dense_vec->Values(); + DBG_ASSERT(vec_vals != NULL); + + vec_vals--; + for( Index i = 0; i < Nonzeros(); i++ ) + { + vec_vals[irows[i]] += std::abs(val[i]); + } +} + +void GenTMatrix::ComputeColA1Impl( + Vector& cols_norms, + bool /*init*/ +) const +{ + DBG_ASSERT(initialized_); + + if( NCols() == 0 ) + { + return; + } + + DenseVector* dense_vec = static_cast(&cols_norms); + DBG_ASSERT(dynamic_cast(&cols_norms)); + + const Index* jcols = Jcols(); + const Number* val = values_; + Number* vec_vals = dense_vec->Values(); + DBG_ASSERT(vec_vals != NULL); + + vec_vals--; // to deal with 1-based indexing in jcols, I believe + for( Index i = 0; i < Nonzeros(); i++ ) + { + vec_vals[jcols[i]] += std::abs(val[i]); + } +} + void GenTMatrix::PrintImplOffset( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/TMatrices/IpGenTMatrix.hpp b/src/LinAlg/TMatrices/IpGenTMatrix.hpp index 5a04572bd..5b37076d8 100644 --- a/src/LinAlg/TMatrices/IpGenTMatrix.hpp +++ b/src/LinAlg/TMatrices/IpGenTMatrix.hpp @@ -119,6 +119,16 @@ class IPOPTLIB_EXPORT GenTMatrix: public Matrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + + virtual void ComputeColA1Impl( + Vector& cols_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/TMatrices/IpSymTMatrix.cpp b/src/LinAlg/TMatrices/IpSymTMatrix.cpp index e85681f35..67b4e1e74 100644 --- a/src/LinAlg/TMatrices/IpSymTMatrix.cpp +++ b/src/LinAlg/TMatrices/IpSymTMatrix.cpp @@ -197,6 +197,42 @@ void SymTMatrix::ComputeRowAMaxImpl( } } +void SymTMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool /*init*/ +) const +{ + DBG_ASSERT(initialized_); + + if( NRows() == 0 ) + { + return; + } + + DenseVector* dense_vec = static_cast(&rows_norms); + DBG_ASSERT(dynamic_cast(&rows_norms)); + + const Index* irn = Irows(); + const Index* jcn = Jcols(); + const Number* val = values_; + Number* vec_vals = dense_vec->Values(); + DBG_ASSERT(vec_vals != NULL); + + // const Number zero = 0.; + // IpBlasCopy(NRows(), &zero, 0, vec_vals, 1); + + vec_vals--; // to deal with 1-based indexing in irn and jcn (I believe) + for( Index i = 0; i < Nonzeros(); i++ ) + { + const Number f = std::abs(*val); + vec_vals[*irn] += f; + vec_vals[*jcn] += f; + val++; + irn++; + jcn++; + } +} + void SymTMatrix::PrintImpl( const Journalist& jnlst, EJournalLevel level, diff --git a/src/LinAlg/TMatrices/IpSymTMatrix.hpp b/src/LinAlg/TMatrices/IpSymTMatrix.hpp index 98e9330dc..e29780306 100644 --- a/src/LinAlg/TMatrices/IpSymTMatrix.hpp +++ b/src/LinAlg/TMatrices/IpSymTMatrix.hpp @@ -129,6 +129,11 @@ class IPOPTLIB_EXPORT SymTMatrix: public SymMatrix bool init ) const; + virtual void ComputeRowA1Impl( + Vector& rows_norms, + bool init + ) const; + virtual void PrintImpl( const Journalist& jnlst, EJournalLevel level, From 72dea2f39d7f4535f1853face8aa4a7ac19e6a65 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Tue, 23 Apr 2024 15:07:04 -0700 Subject: [PATCH 06/20] Remove some commented out code --- src/Algorithm/IpPDFullSpaceSolver.cpp | 31 +-------------------------- 1 file changed, 1 insertion(+), 30 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 42df1a830..3303ad295 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1230,30 +1230,15 @@ bool PDFullSpaceSolver::GMRES( rho0 = rho; } Number resid_nrm_max = resid.Amax(); - // norm-wise relative backward error, in 2 or Inf norm - // Number NRBE_2 = rho / ( A_nrm_inf * x->Nrm2() + b_nrm_2 ); + // norm-wise relative backward error Number NRBE_inf = resid_nrm_max / ( A_nrm_inf * x->Amax() /*+ b_nrm_max*/ ); - // Number NRBE_inf = ComputeResidualRatio(rhs, *x, *V[0], A_nrm_inf); - // auto abs_resid = V[0]->MakeNewIteratesVectorCopy(); - // auto abs_rhs = rhs.MakeNewIteratesVectorCopy(); - // auto abs_x = x->MakeNewIteratesVectorCopy(); - // abs_resid->ElementWiseAbs(); - // abs_rhs->ElementWiseAbs(); - // abs_x->ElementWiseAbs(); - // abs_rhs->Axpy(A_nrm_inf, *abs_x); - // abs_rhs->AddScalar(Number(1.e-16)); - // abs_resid->ElementWiseDivide(*abs_rhs); - // Number CRBE = abs_resid->Amax(); std::cout << "GMRES it. " << totit - // << " NRBE_2= " << std::setw(8) << NRBE_2 << " NRBE_inf= " << std::setw(8) << NRBE_inf - // << " CRBE= " << std::setw(8) << CRBE << " ||r||2= " << std::setw(8) << rho << " ||r||2/||b||2= " << std::setw(8) << rho / b_nrm_2 << " ||r||Inf= " << std::setw(8) << resid_nrm_max << " ||r||Inf/||b||Inf= " << std::setw(8) << resid_nrm_max / b_nrm_max << " ||A||Inf= " << std::setw(8) << A_nrm_inf - // << " ||b||Inf= " << std::setw(8) << b_nrm_max << std::endl; if( NRBE_inf < tol && totit >= min_refinement_steps_) { @@ -1263,7 +1248,6 @@ bool PDFullSpaceSolver::GMRES( } V[0]->Scal( 1./rho ); b_[0] = rho; - // int nrit = restart - 1; for( Index it = 0; it < restart; it++ ) { totit++; @@ -1324,26 +1308,13 @@ bool PDFullSpaceSolver::GMRES( slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, *x, resid); Number resid_nrm_2 = resid.Nrm2(); resid_nrm_max = resid.Amax(); - // NRBE_2 = resid_nrm_2 / ( A_nrm_inf * x->Nrm2() + b_nrm_2 ); NRBE_inf = resid_nrm_max / ( A_nrm_inf * x->Amax() /*+ b_nrm_max*/ ); - // NRBE_inf = ComputeResidualRatio(rhs, *x, resid, A_nrm_inf); - // abs_resid->Copy( resid ); - // abs_x->Copy( *x ); - // abs_resid->ElementWiseAbs(); - // abs_x->ElementWiseAbs(); - // abs_rhs->Axpy(A_nrm_inf, *abs_x); - // abs_rhs->AddScalar(Number(1.e-16)); - // abs_resid->ElementWiseDivide(*abs_rhs); - // Number CRBE = abs_resid->Amax(); std::cout << "GMRES it. " << totit - // << " NRBE_2= " << std::setw(8) << NRBE_2 << " NRBE_inf= " << std::setw(8) << NRBE_inf - // << " CRBE= " << std::setw(8) << CRBE << " ||r||2= " << std::setw(8) << resid_nrm_2 << " ||r||2/||b||2= " << std::setw(8) << resid_nrm_2 / b_nrm_2 << " ||r||Inf= " << std::setw(8) << resid_nrm_max << " ||r||Inf/||b||Inf= " << std::setw(8) << resid_nrm_max / b_nrm_max - // << " rho= " << std::setw(8) << rho << std::endl; if( (NRBE_inf < tol && totit >= min_refinement_steps_) || (totit >= max_refinement_steps_) ) From cecb67029ad679ae2964b4673397c801ef351787 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Tue, 23 Apr 2024 23:27:17 -0700 Subject: [PATCH 07/20] Cleanup in FGMRES code --- src/Algorithm/IpPDFullSpaceSolver.cpp | 93 ++++++++++++++++----------- 1 file changed, 57 insertions(+), 36 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 3303ad295..806892eea 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -401,7 +401,7 @@ bool PDFullSpaceSolver::SolveGMRES( bool improve_solution /* = false */ ) { - DBG_START_METH("PDFullSpaceSolver::Solve", dbg_verbosity); + DBG_START_METH("PDFullSpaceSolver::SolveGMRES", dbg_verbosity); DBG_ASSERT(!allow_inexact || !improve_solution); DBG_ASSERT(!improve_solution || beta == 0.); @@ -1075,11 +1075,13 @@ Number PDFullSpaceSolver::ComputeResidualRatio( } else { - // ToDo: determine how to include norm of matrix, and what - // safeguard to use against incredibly large solution vectors - Number max_cond = 1e6; - return nrm_resid / (Min(nrm_res, max_cond * nrm_rhs) + nrm_rhs); - // return nrm_resid / (AInfNrm * nrm_res /*+ nrm_rhs*/ ); + /* Compute the Inf-norm-wise relative backward error + See: + Nicholas J. Higham. Accuracy and Stability of Numerical + Algorithms. Society for Industrial and Applied Mathematics, + Philadelphia, PA, USA, second edition, 2002. + */ + return nrm_resid / (AInfNrm * nrm_res + nrm_rhs ); } } @@ -1193,12 +1195,8 @@ bool PDFullSpaceSolver::GMRES( { DBG_START_METH("PDFullSpaceSolver::GMRES", dbg_verbosity); - const auto default_precision{std::cout.precision()}; - std::cout << std::setprecision(3); - + // ToDO work on res directly ? auto x = res.MakeNewIteratesVectorCopy(); - auto NullVec = res.MakeNewIteratesVectorCopy(); - NullVec->Set( 0. ); // Index restart = 10; // if (restart > max_refinement_steps_) @@ -1229,17 +1227,29 @@ bool PDFullSpaceSolver::GMRES( { rho0 = rho; } - Number resid_nrm_max = resid.Amax(); - // norm-wise relative backward error - Number NRBE_inf = resid_nrm_max / ( A_nrm_inf * x->Amax() /*+ b_nrm_max*/ ); - std::cout << "GMRES it. " << totit - << " NRBE_inf= " << std::setw(8) << NRBE_inf - << " ||r||2= " << std::setw(8) << rho - << " ||r||2/||b||2= " << std::setw(8) << rho / b_nrm_2 - << " ||r||Inf= " << std::setw(8) << resid_nrm_max - << " ||r||Inf/||b||Inf= " << std::setw(8) << resid_nrm_max / b_nrm_max - << " ||A||Inf= " << std::setw(8) << A_nrm_inf - << std::endl; + + // norm-wise relative backward error, Inf norm + Number NRBE_inf = ComputeResidualRatio( rhs, *x, *V[0], A_nrm_inf ); + + if( Jnlst().ProduceOutput(J_DETAILED, J_LINEAR_ALGEBRA) ) + { + Number resid_nrm_max = resid.Amax(); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "GMRES it = %d\n", totit); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " norm-wise backward error = %e\n", NRBE_inf); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||r||2 = %e\n", rho); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||r||2/||b||2 = %e\n", rho / b_nrm_2); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||r||Inf = %e\n", resid_nrm_max); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||r||Inf/||b||Inf = %e\n", resid_nrm_max / b_nrm_max); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||A||Inf = %e\n", A_nrm_inf); + } + if( NRBE_inf < tol && totit >= min_refinement_steps_) { resid.Copy( *V[0] ); @@ -1263,7 +1273,7 @@ bool PDFullSpaceSolver::GMRES( V.emplace_back( rhs.MakeNewIteratesVector() ); V[it+1]->Set( 0. ); ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, - slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., *NullVec, *Z[it], *V[it+1]); + slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 0., *x, *Z[it], *V[it+1]); { for( Index k=0; k<=it; k++ ) { @@ -1271,7 +1281,9 @@ bool PDFullSpaceSolver::GMRES( V[it+1]->Axpy( -H[k+it*ldH], *V[k] ); } } - { // second time for numerical stability, ToDo optional? + { + // second time for numerical stability, + // ToDo Brown-Hindmarsh condition for( Index k=0; k<=it; k++ ) { auto tmp = V[it+1]->Dot( *V[k] ); @@ -1306,16 +1318,27 @@ bool PDFullSpaceSolver::GMRES( ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, *x, resid); - Number resid_nrm_2 = resid.Nrm2(); - resid_nrm_max = resid.Amax(); - NRBE_inf = resid_nrm_max / ( A_nrm_inf * x->Amax() /*+ b_nrm_max*/ ); - std::cout << "GMRES it. " << totit - << " NRBE_inf= " << std::setw(8) << NRBE_inf - << " ||r||2= " << std::setw(8) << resid_nrm_2 - << " ||r||2/||b||2= " << std::setw(8) << resid_nrm_2 / b_nrm_2 - << " ||r||Inf= " << std::setw(8) << resid_nrm_max - << " ||r||Inf/||b||Inf= " << std::setw(8) << resid_nrm_max / b_nrm_max - << std::endl; + NRBE_inf = ComputeResidualRatio( rhs, *x, resid, A_nrm_inf ); + + if( Jnlst().ProduceOutput(J_DETAILED, J_LINEAR_ALGEBRA) ) + { + Number resid_nrm_2 = resid.Nrm2(); + Number resid_nrm_max = resid.Amax(); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "GMRES it = %d\n", totit); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " norm-wise backward error = %e\n", NRBE_inf); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||r||2 = %e\n", resid_nrm_2); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||r||2/||b||2 = %e\n", resid_nrm_2 / b_nrm_2); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||r||Inf = %e\n", resid_nrm_max); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||r||Inf/||b||Inf = %e\n", resid_nrm_max / b_nrm_max); + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + " ||A||Inf = %e\n", A_nrm_inf); + } if( (NRBE_inf < tol && totit >= min_refinement_steps_) || (totit >= max_refinement_steps_) ) { @@ -1325,8 +1348,6 @@ bool PDFullSpaceSolver::GMRES( } } - std::cout << std::setprecision(default_precision); - res.Copy( *x ); return ( totit < max_refinement_steps_ ); } From 223e08234eceb22e6a8bf312d1b36f5f7afc074e Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Thu, 9 May 2024 09:15:37 -0700 Subject: [PATCH 08/20] Add restart to GMRES (can give better accuracy). Do not use Brown-Hindmarsh condition, can worsen results. --- src/Algorithm/IpPDFullSpaceSolver.cpp | 64 +++++++++++---------------- 1 file changed, 26 insertions(+), 38 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 806892eea..8991d614a 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1195,16 +1195,14 @@ bool PDFullSpaceSolver::GMRES( { DBG_START_METH("PDFullSpaceSolver::GMRES", dbg_verbosity); - // ToDO work on res directly ? - auto x = res.MakeNewIteratesVectorCopy(); - - // Index restart = 10; - // if (restart > max_refinement_steps_) - // restart = max_refinement_steps_; - Index restart = max_refinement_steps_; + Index restart = 3; + if( restart > max_refinement_steps_ ) + { + restart = max_refinement_steps_; + } Index totit = 0, ldH = restart + 1; bool done = false; - Number tol = residual_ratio_max_, rho0 = 0., rho; + Number tol = residual_ratio_max_; Number b_nrm_2 = rhs.Nrm2(), b_nrm_max = rhs.Amax(); Number A_nrm_inf = NrmInf(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U, resid); @@ -1217,19 +1215,14 @@ bool PDFullSpaceSolver::GMRES( V[0]->Set( 0. ); if ( !improve_solution ) { - x->Set(0.); + res.Set(0.); } ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, - slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, *x, *V[0]); - improve_solution = true; // after restart, improve x from previous cycle - rho = V[0]->Nrm2(); - if( !totit ) - { - rho0 = rho; - } - + slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, res, *V[0]); + improve_solution = true; // after restart, improve res from previous cycle + Number rho = V[0]->Nrm2(); // norm-wise relative backward error, Inf norm - Number NRBE_inf = ComputeResidualRatio( rhs, *x, *V[0], A_nrm_inf ); + Number NRBE_inf = ComputeResidualRatio( rhs, res, *V[0], A_nrm_inf ); if( Jnlst().ProduceOutput(J_DETAILED, J_LINEAR_ALGEBRA) ) { @@ -1273,26 +1266,23 @@ bool PDFullSpaceSolver::GMRES( V.emplace_back( rhs.MakeNewIteratesVector() ); V[it+1]->Set( 0. ); ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, - slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 0., *x, *Z[it], *V[it+1]); + slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 0., res, *Z[it], *V[it+1]); + for( Index k=0; k<=it; k++ ) { - for( Index k=0; k<=it; k++ ) - { - H[k+it*ldH] = V[it+1]->Dot( *V[k] ); - V[it+1]->Axpy( -H[k+it*ldH], *V[k] ); - } + H[k+it*ldH] = V[it+1]->Dot( *V[k] ); + V[it+1]->Axpy( -H[k+it*ldH], *V[k] ); } + for( Index k=0; k<=it; k++ ) { - // second time for numerical stability, - // ToDo Brown-Hindmarsh condition - for( Index k=0; k<=it; k++ ) - { - auto tmp = V[it+1]->Dot( *V[k] ); - H[k+it*ldH] += tmp; - V[it+1]->Axpy( -tmp, *V[k] ); - } + auto tmp = V[it+1]->Dot( *V[k] ); + H[k+it*ldH] += tmp; + V[it+1]->Axpy( -tmp, *V[k] ); } H[it+1+it*ldH] = V[it+1]->Nrm2(); - V[it+1]->Scal( 1. / H[it+1+it*ldH] ); + if( H[it+1+it*ldH] != 0. ) + { + V[it+1]->Scal( 1. / H[it+1+it*ldH] ); + } for( Index k = 1; k < it+1; k++ ) { auto gamma = givens_c[k-1]*H[k-1+it*ldH] + givens_s[k-1]*H[k+it*ldH]; @@ -1305,7 +1295,6 @@ bool PDFullSpaceSolver::GMRES( H[it+it*ldH] = givens_c[it]*H[it+it*ldH] + givens_s[it]*H[it+1+it*ldH]; b_[it+1] = -givens_s[it]*b_[it]; b_[it] = givens_c[it]*b_[it]; - rho = std::abs( b_[it+1] ); std::vector y = b_; for( Index k = it; k >= 0; k-- ) @@ -1314,11 +1303,11 @@ bool PDFullSpaceSolver::GMRES( y[k] -= H[k+i*ldH] * y[i]; y[k] /= H[k+k*ldH]; } - x->Axpy( y[it], *Z[it] ); + res.Axpy( y[it], *Z[it] ); ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, - slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, *x, resid); - NRBE_inf = ComputeResidualRatio( rhs, *x, resid, A_nrm_inf ); + slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, res, resid); + NRBE_inf = ComputeResidualRatio( rhs, res, resid, A_nrm_inf ); if( Jnlst().ProduceOutput(J_DETAILED, J_LINEAR_ALGEBRA) ) { @@ -1348,7 +1337,6 @@ bool PDFullSpaceSolver::GMRES( } } - res.Copy( *x ); return ( totit < max_refinement_steps_ ); } From c6a8be22ecc6575657eb72b2040f0f131d5eea68 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Thu, 20 Jun 2024 18:29:47 -0700 Subject: [PATCH 09/20] Comment on the GMRES restart length --- src/Algorithm/IpPDFullSpaceSolver.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 8991d614a..f94eec15b 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1195,6 +1195,12 @@ bool PDFullSpaceSolver::GMRES( { DBG_START_METH("PDFullSpaceSolver::GMRES", dbg_verbosity); + // A larger GMRES restart is often better for convergence, but + // convergence should be fast since we use a direct solver as + // preconditioner. And restarting can improve attainable accuracy, + // see: Buttari, Alfredo, Nicholas J. Higham, Theo Mary, and + // Bastien Vieuble. "A modular framework for the backward error + // analysis of GMRES." (2024). Index restart = 3; if( restart > max_refinement_steps_ ) { From 8f94f83db136e2138925e7b57c78ba2ea3ab2117 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Fri, 21 Jun 2024 08:46:44 -0700 Subject: [PATCH 10/20] Remove unused header include --- src/Algorithm/IpPDFullSpaceSolver.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index f94eec15b..2f03d29ad 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -8,7 +8,6 @@ #include "IpDebug.hpp" #include -#include namespace Ipopt { From 326f2f3b83cb207beaa2361c992a57e705793783 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Fri, 21 Jun 2024 08:48:07 -0700 Subject: [PATCH 11/20] Implement missing matrix row sum computation for SymScaledMatrix --- src/LinAlg/IpSymScaledMatrix.cpp | 36 +++++++++++++++++++++++++++++--- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/src/LinAlg/IpSymScaledMatrix.cpp b/src/LinAlg/IpSymScaledMatrix.cpp index ce32f1b7a..d8f66dd93 100644 --- a/src/LinAlg/IpSymScaledMatrix.cpp +++ b/src/LinAlg/IpSymScaledMatrix.cpp @@ -5,6 +5,8 @@ // Authors: Carl Laird, Andreas Waechter IBM 2004-08-13 #include "IpSymScaledMatrix.hpp" +#include "IpSymTMatrix.hpp" +#include "IpDenseVector.hpp" namespace Ipopt { @@ -71,11 +73,39 @@ void SymScaledMatrix::ComputeRowAMaxImpl( } void SymScaledMatrix::ComputeRowA1Impl( - Vector& /*rows_norms*/, - bool /*init*/ + Vector& rows_norms, + bool init ) const { - THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SymScaledMatrix::ComputeRowA1Impl not implemented"); + DBG_ASSERT(IsValid(matrix_)); + + if( IsValid(owner_space_->RowColScaling()) ) + { + const SymTMatrix* mt = dynamic_cast(GetRawPtr(matrix_)); + if( mt ) + { + const DenseVector* D = dynamic_cast(GetRawPtr(owner_space_->RowColScaling())); + DenseVector* rn = dynamic_cast(&rows_norms); + DBG_ASSERT(D); + DBG_ASSERT(rn); + for ( Index i = 0; i < mt->Nonzeros(); i++ ) + { + rn->Values()[mt->Irows()[i] - 1] += std::abs(mt->Values()[i] * D->Values()[mt->Jcols()[i] - 1]); + if( mt->Irows()[i] != mt->Jcols()[i] ) + { + rn->Values()[mt->Jcols()[i] - 1] += std::abs(mt->Values()[i] * D->Values()[mt->Irows()[i] - 1]); + } + } + } + else + { + THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ScaledMatrix::ComputeRowA1Impl not implemented"); + } + } + else + { + matrix_->ComputeRowA1(rows_norms, init); + } } void SymScaledMatrix::PrintImpl( From 6e0f2310eeaf35170e770ce74561546df1e792b2 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Fri, 21 Jun 2024 10:02:48 -0700 Subject: [PATCH 12/20] Continue with matrix norm (under-)estimate if the norm cannot be computed accurately --- src/Algorithm/IpPDFullSpaceSolver.cpp | 51 ++++++++++++++++++++------- 1 file changed, 38 insertions(+), 13 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 2f03d29ad..da2d52d35 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1113,24 +1113,43 @@ Number PDFullSpaceSolver::NrmInf( perturbHandler_->CurrentPerturbation(delta_x, delta_s, delta_c, delta_d); tmp.Set( 0. ); - W.ComputeRowA1(*tmp.x_NonConst(), false); - J_c.ComputeColA1(*tmp.x_NonConst(), true); - J_d.ComputeColA1(*tmp.x_NonConst(), true); - Px_L.ComputeRowA1(*tmp.x_NonConst(), true); - Px_U.ComputeRowA1(*tmp.x_NonConst(), true); + bool norm_underestimated = false; + try + { + W.ComputeRowA1(*tmp.x_NonConst()); + } catch (UNIMPLEMENTED_LINALG_METHOD_CALLED) + { + norm_underestimated = true; + } + try + { + J_c.ComputeColA1(*tmp.x_NonConst(), false); + } catch (UNIMPLEMENTED_LINALG_METHOD_CALLED) + { + norm_underestimated = true; + } + try + { + J_d.ComputeColA1(*tmp.x_NonConst(), false); + } catch (UNIMPLEMENTED_LINALG_METHOD_CALLED) + { + norm_underestimated = true; + } + Px_L.ComputeRowA1(*tmp.x_NonConst(), false); + Px_U.ComputeRowA1(*tmp.x_NonConst(), false); Number A_inf = tmp.x_NonConst()->Amax() + std::abs(delta_x); - Pd_L.ComputeRowA1(*tmp.s_NonConst(), false); - Pd_U.ComputeRowA1(*tmp.s_NonConst(), true); + Pd_L.ComputeRowA1(*tmp.s_NonConst()); + Pd_U.ComputeRowA1(*tmp.s_NonConst(), false); A_inf = std::max(A_inf, tmp.s_NonConst()->Amax() + 1. + delta_s); - J_c.ComputeRowA1(*tmp.y_c_NonConst(), false); + J_c.ComputeRowA1(*tmp.y_c_NonConst()); A_inf = std::max(A_inf, tmp.y_c_NonConst()->Amax() + delta_c); - J_d.ComputeRowA1(*tmp.y_d_NonConst(), false); + J_d.ComputeRowA1(*tmp.y_d_NonConst()); A_inf = std::max(A_inf, tmp.y_d_NonConst()->Amax() + 1. + delta_d); - Px_L.ComputeColA1(*tmp.z_L_NonConst(), false); + Px_L.ComputeColA1(*tmp.z_L_NonConst()); tmp.z_L_NonConst()->ElementWiseMultiply(z_L); tmp.z_L_NonConst()->ElementWiseAbs(); auto tmp_slack_x_L = slack_x_L.MakeNewCopy(); @@ -1138,7 +1157,7 @@ Number PDFullSpaceSolver::NrmInf( tmp.z_L_NonConst()->Axpy(1., *tmp_slack_x_L); A_inf = std::max(A_inf, tmp.z_L_NonConst()->Amax()); - Px_U.ComputeColA1(*tmp.z_U_NonConst(), false); + Px_U.ComputeColA1(*tmp.z_U_NonConst()); tmp.z_U_NonConst()->ElementWiseMultiply(z_U); tmp.z_U_NonConst()->ElementWiseAbs(); auto tmp_slack_x_U = slack_x_U.MakeNewCopy(); @@ -1146,7 +1165,7 @@ Number PDFullSpaceSolver::NrmInf( tmp.z_U_NonConst()->Axpy(1., *tmp_slack_x_U); A_inf = std::max(A_inf, tmp.z_U_NonConst()->Amax()); - Pd_L.ComputeColA1(*tmp.v_L_NonConst(), false); + Pd_L.ComputeColA1(*tmp.v_L_NonConst()); tmp.v_L_NonConst()->ElementWiseMultiply(v_L); tmp.v_L_NonConst()->ElementWiseAbs(); auto tmp_slack_s_L = slack_s_L.MakeNewCopy(); @@ -1154,7 +1173,7 @@ Number PDFullSpaceSolver::NrmInf( tmp.v_L_NonConst()->Axpy(1., *tmp_slack_s_L); A_inf = std::max(A_inf, tmp.v_L_NonConst()->Amax()); - Pd_U.ComputeColA1(*tmp.v_U_NonConst(), false); + Pd_U.ComputeColA1(*tmp.v_U_NonConst()); tmp.v_U_NonConst()->ElementWiseMultiply(v_U); tmp.v_U_NonConst()->ElementWiseAbs(); auto tmp_slack_s_U = slack_s_U.MakeNewCopy(); @@ -1162,6 +1181,12 @@ Number PDFullSpaceSolver::NrmInf( tmp.v_U_NonConst()->Axpy(1., *tmp_slack_s_U); A_inf = std::max(A_inf, tmp.v_U_NonConst()->Amax()); + if( norm_underestimated ) + { + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "Cannot accurately compute matrix norm!\n"); + } + return A_inf; } From 533c7175f313b284db3a91604f76b1a2bf84f251 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Fri, 21 Jun 2024 13:53:10 -0700 Subject: [PATCH 13/20] Cast to Number --- src/Algorithm/IpPDFullSpaceSolver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index da2d52d35..0be085dd5 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1141,13 +1141,13 @@ Number PDFullSpaceSolver::NrmInf( Pd_L.ComputeRowA1(*tmp.s_NonConst()); Pd_U.ComputeRowA1(*tmp.s_NonConst(), false); - A_inf = std::max(A_inf, tmp.s_NonConst()->Amax() + 1. + delta_s); + A_inf = std::max(A_inf, tmp.s_NonConst()->Amax() + Number(1.) + delta_s); J_c.ComputeRowA1(*tmp.y_c_NonConst()); A_inf = std::max(A_inf, tmp.y_c_NonConst()->Amax() + delta_c); J_d.ComputeRowA1(*tmp.y_d_NonConst()); - A_inf = std::max(A_inf, tmp.y_d_NonConst()->Amax() + 1. + delta_d); + A_inf = std::max(A_inf, tmp.y_d_NonConst()->Amax() + Number(1.) + delta_d); Px_L.ComputeColA1(*tmp.z_L_NonConst()); tmp.z_L_NonConst()->ElementWiseMultiply(z_L); From 94ee41775d417fc30214c9c8008450d1205ffe71 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Thu, 27 Jun 2024 09:48:27 -0700 Subject: [PATCH 14/20] Avoid divide by zero, happy breakdown, in GMRES --- src/Algorithm/IpPDFullSpaceSolver.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 0be085dd5..25cb38e9e 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1273,7 +1273,8 @@ bool PDFullSpaceSolver::GMRES( " ||A||Inf = %e\n", A_nrm_inf); } - if( NRBE_inf < tol && totit >= min_refinement_steps_) + if( ( NRBE_inf < tol && totit >= min_refinement_steps_ ) || + rho < std::numeric_limits::epsilon() ) { resid.Copy( *V[0] ); done = true; From 530ae0a189c7550d8d822a2abea7e41da74d11e4 Mon Sep 17 00:00:00 2001 From: Pieter Ghysels Date: Wed, 3 Jul 2024 13:48:00 -0700 Subject: [PATCH 15/20] Add some printing to debug AppVeyor test --- src/Algorithm/IpPDFullSpaceSolver.cpp | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 25cb38e9e..dfcad6ae8 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1254,9 +1254,18 @@ bool PDFullSpaceSolver::GMRES( // norm-wise relative backward error, Inf norm Number NRBE_inf = ComputeResidualRatio( rhs, res, *V[0], A_nrm_inf ); + Number resid_nrm_max = resid.Amax(); + std::cout << "GMRES it = " << totit << std::endl + << " norm-wise backward error = " << NRBE_inf << std::endl + << " ||r||2 = resid_nrm_2" << std::endl + << " ||r||2/||b||2 = " << rho / b_nrm_2 << std::endl + << " ||r||Inf = " << resid_nrm_max << std::endl + << " ||r||Inf/||b||Inf = " << resid_nrm_max / b_nrm_max << std::endl + << " ||A||Inf = " << A_nrm_inf << std::endl; + if( Jnlst().ProduceOutput(J_DETAILED, J_LINEAR_ALGEBRA) ) { - Number resid_nrm_max = resid.Amax(); + // Number resid_nrm_max = resid.Amax(); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "GMRES it = %d\n", totit); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, @@ -1340,10 +1349,20 @@ bool PDFullSpaceSolver::GMRES( slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, res, resid); NRBE_inf = ComputeResidualRatio( rhs, res, resid, A_nrm_inf ); + Number resid_nrm_2 = resid.Nrm2(); + Number resid_nrm_max = resid.Amax(); + std::cout << "GMRES it = " << totit << std::endl + << " norm-wise backward error = " << NRBE_inf << std::endl + << " ||r||2 = resid_nrm_2" << std::endl + << " ||r||2/||b||2 = " << resid_nrm_2 / b_nrm_2 << std::endl + << " ||r||Inf = " << resid_nrm_max << std::endl + << " ||r||Inf/||b||Inf = " << resid_nrm_max / b_nrm_max << std::endl + << " ||A||Inf = " << A_nrm_inf << std::endl; + if( Jnlst().ProduceOutput(J_DETAILED, J_LINEAR_ALGEBRA) ) { - Number resid_nrm_2 = resid.Nrm2(); - Number resid_nrm_max = resid.Amax(); + // Number resid_nrm_2 = resid.Nrm2(); + // Number resid_nrm_max = resid.Amax(); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "GMRES it = %d\n", totit); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, From cafb7e8668a11c46a8fe33b1deee2a9bc453f84f Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Wed, 17 Jul 2024 15:39:56 -0600 Subject: [PATCH 16/20] adding residual_improvement_factor_; remove use of auto --- src/Algorithm/IpPDFullSpaceSolver.cpp | 29 +++++++++++++++++---------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index dfcad6ae8..fc70d4b4b 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1152,7 +1152,7 @@ Number PDFullSpaceSolver::NrmInf( Px_L.ComputeColA1(*tmp.z_L_NonConst()); tmp.z_L_NonConst()->ElementWiseMultiply(z_L); tmp.z_L_NonConst()->ElementWiseAbs(); - auto tmp_slack_x_L = slack_x_L.MakeNewCopy(); + SmartPtr tmp_slack_x_L = slack_x_L.MakeNewCopy(); tmp_slack_x_L->ElementWiseAbs(); tmp.z_L_NonConst()->Axpy(1., *tmp_slack_x_L); A_inf = std::max(A_inf, tmp.z_L_NonConst()->Amax()); @@ -1160,7 +1160,7 @@ Number PDFullSpaceSolver::NrmInf( Px_U.ComputeColA1(*tmp.z_U_NonConst()); tmp.z_U_NonConst()->ElementWiseMultiply(z_U); tmp.z_U_NonConst()->ElementWiseAbs(); - auto tmp_slack_x_U = slack_x_U.MakeNewCopy(); + SmartPtr tmp_slack_x_U = slack_x_U.MakeNewCopy(); tmp_slack_x_U->ElementWiseAbs(); tmp.z_U_NonConst()->Axpy(1., *tmp_slack_x_U); A_inf = std::max(A_inf, tmp.z_U_NonConst()->Amax()); @@ -1168,7 +1168,7 @@ Number PDFullSpaceSolver::NrmInf( Pd_L.ComputeColA1(*tmp.v_L_NonConst()); tmp.v_L_NonConst()->ElementWiseMultiply(v_L); tmp.v_L_NonConst()->ElementWiseAbs(); - auto tmp_slack_s_L = slack_s_L.MakeNewCopy(); + SmartPtr tmp_slack_s_L = slack_s_L.MakeNewCopy(); tmp_slack_s_L->ElementWiseAbs(); tmp.v_L_NonConst()->Axpy(1., *tmp_slack_s_L); A_inf = std::max(A_inf, tmp.v_L_NonConst()->Amax()); @@ -1176,7 +1176,7 @@ Number PDFullSpaceSolver::NrmInf( Pd_U.ComputeColA1(*tmp.v_U_NonConst()); tmp.v_U_NonConst()->ElementWiseMultiply(v_U); tmp.v_U_NonConst()->ElementWiseAbs(); - auto tmp_slack_s_U = slack_s_U.MakeNewCopy(); + SmartPtr tmp_slack_s_U = slack_s_U.MakeNewCopy(); tmp_slack_s_U->ElementWiseAbs(); tmp.v_U_NonConst()->Axpy(1., *tmp_slack_s_U); A_inf = std::max(A_inf, tmp.v_U_NonConst()->Amax()); @@ -1236,11 +1236,14 @@ bool PDFullSpaceSolver::GMRES( Number b_nrm_2 = rhs.Nrm2(), b_nrm_max = rhs.Amax(); Number A_nrm_inf = NrmInf(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U, resid); + // for first loop + Number residual_ratio_old = std::numeric_limits::max(); + Number NRBE_inf = std::numeric_limits::max(); while ( !done ) { std::vector givens_c(restart), givens_s(restart), b_(restart+1), H(restart*(restart+1)); - std::vector> V, Z; + std::vector > V, Z; V.emplace_back(rhs.MakeNewIteratesVector()); V[0]->Set( 0. ); if ( !improve_solution ) @@ -1251,8 +1254,9 @@ bool PDFullSpaceSolver::GMRES( slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, res, *V[0]); improve_solution = true; // after restart, improve res from previous cycle Number rho = V[0]->Nrm2(); + residual_ratio_old = NRBE_inf; // norm-wise relative backward error, Inf norm - Number NRBE_inf = ComputeResidualRatio( rhs, res, *V[0], A_nrm_inf ); + NRBE_inf = ComputeResidualRatio( rhs, res, *V[0], A_nrm_inf ); Number resid_nrm_max = resid.Amax(); std::cout << "GMRES it = " << totit << std::endl @@ -1283,7 +1287,8 @@ bool PDFullSpaceSolver::GMRES( } if( ( NRBE_inf < tol && totit >= min_refinement_steps_ ) || - rho < std::numeric_limits::epsilon() ) + rho < std::numeric_limits::epsilon() || + NRBE_inf > residual_improvement_factor_ * residual_ratio_old ) { resid.Copy( *V[0] ); done = true; @@ -1314,7 +1319,7 @@ bool PDFullSpaceSolver::GMRES( } for( Index k=0; k<=it; k++ ) { - auto tmp = V[it+1]->Dot( *V[k] ); + Number tmp = V[it+1]->Dot( *V[k] ); H[k+it*ldH] += tmp; V[it+1]->Axpy( -tmp, *V[k] ); } @@ -1325,11 +1330,11 @@ bool PDFullSpaceSolver::GMRES( } for( Index k = 1; k < it+1; k++ ) { - auto gamma = givens_c[k-1]*H[k-1+it*ldH] + givens_s[k-1]*H[k+it*ldH]; + Number gamma = givens_c[k-1]*H[k-1+it*ldH] + givens_s[k-1]*H[k+it*ldH]; H[k+it*ldH] = -givens_s[k-1]*H[k-1+it*ldH] + givens_c[k-1]*H[k+it*ldH]; H[k-1+it*ldH] = gamma; } - auto delta = std::sqrt( H[it+it*ldH]*H[it+it*ldH] + H[it+1+it*ldH]*H[it+1+it*ldH] ); + Number delta = std::sqrt( H[it+it*ldH]*H[it+it*ldH] + H[it+1+it*ldH]*H[it+1+it*ldH] ); givens_c[it] = H[it+it*ldH] / delta; givens_s[it] = H[it+1+it*ldH] / delta; H[it+it*ldH] = givens_c[it]*H[it+it*ldH] + givens_s[it]*H[it+1+it*ldH]; @@ -1347,6 +1352,7 @@ bool PDFullSpaceSolver::GMRES( ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, res, resid); + residual_ratio_old = NRBE_inf; NRBE_inf = ComputeResidualRatio( rhs, res, resid, A_nrm_inf ); Number resid_nrm_2 = resid.Nrm2(); @@ -1379,7 +1385,8 @@ bool PDFullSpaceSolver::GMRES( " ||A||Inf = %e\n", A_nrm_inf); } if( (NRBE_inf < tol && totit >= min_refinement_steps_) || - (totit >= max_refinement_steps_) ) + (totit >= max_refinement_steps_) || + (NRBE_inf > residual_improvement_factor_ * residual_ratio_old ) ) { done = true; break; From 0bf732c9f72dfd5915cd977a0c6370eccacb1f58 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Thu, 18 Jul 2024 12:39:29 -0600 Subject: [PATCH 17/20] Revert "adding residual_improvement_factor_; remove use of auto" This reverts commit cafb7e8668a11c46a8fe33b1deee2a9bc453f84f. --- src/Algorithm/IpPDFullSpaceSolver.cpp | 29 ++++++++++----------------- 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index fc70d4b4b..dfcad6ae8 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1152,7 +1152,7 @@ Number PDFullSpaceSolver::NrmInf( Px_L.ComputeColA1(*tmp.z_L_NonConst()); tmp.z_L_NonConst()->ElementWiseMultiply(z_L); tmp.z_L_NonConst()->ElementWiseAbs(); - SmartPtr tmp_slack_x_L = slack_x_L.MakeNewCopy(); + auto tmp_slack_x_L = slack_x_L.MakeNewCopy(); tmp_slack_x_L->ElementWiseAbs(); tmp.z_L_NonConst()->Axpy(1., *tmp_slack_x_L); A_inf = std::max(A_inf, tmp.z_L_NonConst()->Amax()); @@ -1160,7 +1160,7 @@ Number PDFullSpaceSolver::NrmInf( Px_U.ComputeColA1(*tmp.z_U_NonConst()); tmp.z_U_NonConst()->ElementWiseMultiply(z_U); tmp.z_U_NonConst()->ElementWiseAbs(); - SmartPtr tmp_slack_x_U = slack_x_U.MakeNewCopy(); + auto tmp_slack_x_U = slack_x_U.MakeNewCopy(); tmp_slack_x_U->ElementWiseAbs(); tmp.z_U_NonConst()->Axpy(1., *tmp_slack_x_U); A_inf = std::max(A_inf, tmp.z_U_NonConst()->Amax()); @@ -1168,7 +1168,7 @@ Number PDFullSpaceSolver::NrmInf( Pd_L.ComputeColA1(*tmp.v_L_NonConst()); tmp.v_L_NonConst()->ElementWiseMultiply(v_L); tmp.v_L_NonConst()->ElementWiseAbs(); - SmartPtr tmp_slack_s_L = slack_s_L.MakeNewCopy(); + auto tmp_slack_s_L = slack_s_L.MakeNewCopy(); tmp_slack_s_L->ElementWiseAbs(); tmp.v_L_NonConst()->Axpy(1., *tmp_slack_s_L); A_inf = std::max(A_inf, tmp.v_L_NonConst()->Amax()); @@ -1176,7 +1176,7 @@ Number PDFullSpaceSolver::NrmInf( Pd_U.ComputeColA1(*tmp.v_U_NonConst()); tmp.v_U_NonConst()->ElementWiseMultiply(v_U); tmp.v_U_NonConst()->ElementWiseAbs(); - SmartPtr tmp_slack_s_U = slack_s_U.MakeNewCopy(); + auto tmp_slack_s_U = slack_s_U.MakeNewCopy(); tmp_slack_s_U->ElementWiseAbs(); tmp.v_U_NonConst()->Axpy(1., *tmp_slack_s_U); A_inf = std::max(A_inf, tmp.v_U_NonConst()->Amax()); @@ -1236,14 +1236,11 @@ bool PDFullSpaceSolver::GMRES( Number b_nrm_2 = rhs.Nrm2(), b_nrm_max = rhs.Amax(); Number A_nrm_inf = NrmInf(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U, resid); - // for first loop - Number residual_ratio_old = std::numeric_limits::max(); - Number NRBE_inf = std::numeric_limits::max(); while ( !done ) { std::vector givens_c(restart), givens_s(restart), b_(restart+1), H(restart*(restart+1)); - std::vector > V, Z; + std::vector> V, Z; V.emplace_back(rhs.MakeNewIteratesVector()); V[0]->Set( 0. ); if ( !improve_solution ) @@ -1254,9 +1251,8 @@ bool PDFullSpaceSolver::GMRES( slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, res, *V[0]); improve_solution = true; // after restart, improve res from previous cycle Number rho = V[0]->Nrm2(); - residual_ratio_old = NRBE_inf; // norm-wise relative backward error, Inf norm - NRBE_inf = ComputeResidualRatio( rhs, res, *V[0], A_nrm_inf ); + Number NRBE_inf = ComputeResidualRatio( rhs, res, *V[0], A_nrm_inf ); Number resid_nrm_max = resid.Amax(); std::cout << "GMRES it = " << totit << std::endl @@ -1287,8 +1283,7 @@ bool PDFullSpaceSolver::GMRES( } if( ( NRBE_inf < tol && totit >= min_refinement_steps_ ) || - rho < std::numeric_limits::epsilon() || - NRBE_inf > residual_improvement_factor_ * residual_ratio_old ) + rho < std::numeric_limits::epsilon() ) { resid.Copy( *V[0] ); done = true; @@ -1319,7 +1314,7 @@ bool PDFullSpaceSolver::GMRES( } for( Index k=0; k<=it; k++ ) { - Number tmp = V[it+1]->Dot( *V[k] ); + auto tmp = V[it+1]->Dot( *V[k] ); H[k+it*ldH] += tmp; V[it+1]->Axpy( -tmp, *V[k] ); } @@ -1330,11 +1325,11 @@ bool PDFullSpaceSolver::GMRES( } for( Index k = 1; k < it+1; k++ ) { - Number gamma = givens_c[k-1]*H[k-1+it*ldH] + givens_s[k-1]*H[k+it*ldH]; + auto gamma = givens_c[k-1]*H[k-1+it*ldH] + givens_s[k-1]*H[k+it*ldH]; H[k+it*ldH] = -givens_s[k-1]*H[k-1+it*ldH] + givens_c[k-1]*H[k+it*ldH]; H[k-1+it*ldH] = gamma; } - Number delta = std::sqrt( H[it+it*ldH]*H[it+it*ldH] + H[it+1+it*ldH]*H[it+1+it*ldH] ); + auto delta = std::sqrt( H[it+it*ldH]*H[it+it*ldH] + H[it+1+it*ldH]*H[it+1+it*ldH] ); givens_c[it] = H[it+it*ldH] / delta; givens_s[it] = H[it+1+it*ldH] / delta; H[it+it*ldH] = givens_c[it]*H[it+it*ldH] + givens_s[it]*H[it+1+it*ldH]; @@ -1352,7 +1347,6 @@ bool PDFullSpaceSolver::GMRES( ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, z_U, v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, res, resid); - residual_ratio_old = NRBE_inf; NRBE_inf = ComputeResidualRatio( rhs, res, resid, A_nrm_inf ); Number resid_nrm_2 = resid.Nrm2(); @@ -1385,8 +1379,7 @@ bool PDFullSpaceSolver::GMRES( " ||A||Inf = %e\n", A_nrm_inf); } if( (NRBE_inf < tol && totit >= min_refinement_steps_) || - (totit >= max_refinement_steps_) || - (NRBE_inf > residual_improvement_factor_ * residual_ratio_old ) ) + (totit >= max_refinement_steps_) ) { done = true; break; From 610b9896c6122d6009a78c86f8ae6cd08b0e92c1 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Thu, 18 Jul 2024 12:38:14 -0600 Subject: [PATCH 18/20] stop if NRBE_inf is less than machine epsilon --- src/Algorithm/IpPDFullSpaceSolver.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index dfcad6ae8..1bd99fd38 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1257,7 +1257,7 @@ bool PDFullSpaceSolver::GMRES( Number resid_nrm_max = resid.Amax(); std::cout << "GMRES it = " << totit << std::endl << " norm-wise backward error = " << NRBE_inf << std::endl - << " ||r||2 = resid_nrm_2" << std::endl + << " ||r||2 = " << rho << std::endl << " ||r||2/||b||2 = " << rho / b_nrm_2 << std::endl << " ||r||Inf = " << resid_nrm_max << std::endl << " ||r||Inf/||b||Inf = " << resid_nrm_max / b_nrm_max << std::endl @@ -1283,6 +1283,7 @@ bool PDFullSpaceSolver::GMRES( } if( ( NRBE_inf < tol && totit >= min_refinement_steps_ ) || + NRBE_inf < std::numeric_limits::epsilon() || rho < std::numeric_limits::epsilon() ) { resid.Copy( *V[0] ); @@ -1353,7 +1354,7 @@ bool PDFullSpaceSolver::GMRES( Number resid_nrm_max = resid.Amax(); std::cout << "GMRES it = " << totit << std::endl << " norm-wise backward error = " << NRBE_inf << std::endl - << " ||r||2 = resid_nrm_2" << std::endl + << " ||r||2 = " << resid_nrm_2 << std::endl << " ||r||2/||b||2 = " << resid_nrm_2 / b_nrm_2 << std::endl << " ||r||Inf = " << resid_nrm_max << std::endl << " ||r||Inf/||b||Inf = " << resid_nrm_max / b_nrm_max << std::endl @@ -1379,6 +1380,7 @@ bool PDFullSpaceSolver::GMRES( " ||A||Inf = %e\n", A_nrm_inf); } if( (NRBE_inf < tol && totit >= min_refinement_steps_) || + NRBE_inf < std::numeric_limits::epsilon() || (totit >= max_refinement_steps_) ) { done = true; From 5639fd238e1b14b27265b686270d5587cec681be Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Thu, 18 Jul 2024 13:56:19 -0600 Subject: [PATCH 19/20] Revert "Add some printing to debug AppVeyor test" This reverts commit 530ae0a189c7550d8d822a2abea7e41da74d11e4. --- src/Algorithm/IpPDFullSpaceSolver.cpp | 25 +++---------------------- 1 file changed, 3 insertions(+), 22 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 1bd99fd38..d25107ece 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1254,18 +1254,9 @@ bool PDFullSpaceSolver::GMRES( // norm-wise relative backward error, Inf norm Number NRBE_inf = ComputeResidualRatio( rhs, res, *V[0], A_nrm_inf ); - Number resid_nrm_max = resid.Amax(); - std::cout << "GMRES it = " << totit << std::endl - << " norm-wise backward error = " << NRBE_inf << std::endl - << " ||r||2 = " << rho << std::endl - << " ||r||2/||b||2 = " << rho / b_nrm_2 << std::endl - << " ||r||Inf = " << resid_nrm_max << std::endl - << " ||r||Inf/||b||Inf = " << resid_nrm_max / b_nrm_max << std::endl - << " ||A||Inf = " << A_nrm_inf << std::endl; - if( Jnlst().ProduceOutput(J_DETAILED, J_LINEAR_ALGEBRA) ) { - // Number resid_nrm_max = resid.Amax(); + Number resid_nrm_max = resid.Amax(); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "GMRES it = %d\n", totit); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, @@ -1350,20 +1341,10 @@ bool PDFullSpaceSolver::GMRES( slack_s_L, slack_s_U, sigma_x, sigma_s, -1., 1., rhs, res, resid); NRBE_inf = ComputeResidualRatio( rhs, res, resid, A_nrm_inf ); - Number resid_nrm_2 = resid.Nrm2(); - Number resid_nrm_max = resid.Amax(); - std::cout << "GMRES it = " << totit << std::endl - << " norm-wise backward error = " << NRBE_inf << std::endl - << " ||r||2 = " << resid_nrm_2 << std::endl - << " ||r||2/||b||2 = " << resid_nrm_2 / b_nrm_2 << std::endl - << " ||r||Inf = " << resid_nrm_max << std::endl - << " ||r||Inf/||b||Inf = " << resid_nrm_max / b_nrm_max << std::endl - << " ||A||Inf = " << A_nrm_inf << std::endl; - if( Jnlst().ProduceOutput(J_DETAILED, J_LINEAR_ALGEBRA) ) { - // Number resid_nrm_2 = resid.Nrm2(); - // Number resid_nrm_max = resid.Amax(); + Number resid_nrm_2 = resid.Nrm2(); + Number resid_nrm_max = resid.Amax(); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "GMRES it = %d\n", totit); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, From 2aa4d4447eb343139bd617fe7c23a6d82effee67 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Thu, 18 Jul 2024 13:58:42 -0600 Subject: [PATCH 20/20] remove auto --- src/Algorithm/IpPDFullSpaceSolver.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index d25107ece..062042c92 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -1152,7 +1152,7 @@ Number PDFullSpaceSolver::NrmInf( Px_L.ComputeColA1(*tmp.z_L_NonConst()); tmp.z_L_NonConst()->ElementWiseMultiply(z_L); tmp.z_L_NonConst()->ElementWiseAbs(); - auto tmp_slack_x_L = slack_x_L.MakeNewCopy(); + SmartPtr tmp_slack_x_L = slack_x_L.MakeNewCopy(); tmp_slack_x_L->ElementWiseAbs(); tmp.z_L_NonConst()->Axpy(1., *tmp_slack_x_L); A_inf = std::max(A_inf, tmp.z_L_NonConst()->Amax()); @@ -1160,7 +1160,7 @@ Number PDFullSpaceSolver::NrmInf( Px_U.ComputeColA1(*tmp.z_U_NonConst()); tmp.z_U_NonConst()->ElementWiseMultiply(z_U); tmp.z_U_NonConst()->ElementWiseAbs(); - auto tmp_slack_x_U = slack_x_U.MakeNewCopy(); + SmartPtr tmp_slack_x_U = slack_x_U.MakeNewCopy(); tmp_slack_x_U->ElementWiseAbs(); tmp.z_U_NonConst()->Axpy(1., *tmp_slack_x_U); A_inf = std::max(A_inf, tmp.z_U_NonConst()->Amax()); @@ -1168,7 +1168,7 @@ Number PDFullSpaceSolver::NrmInf( Pd_L.ComputeColA1(*tmp.v_L_NonConst()); tmp.v_L_NonConst()->ElementWiseMultiply(v_L); tmp.v_L_NonConst()->ElementWiseAbs(); - auto tmp_slack_s_L = slack_s_L.MakeNewCopy(); + SmartPtr tmp_slack_s_L = slack_s_L.MakeNewCopy(); tmp_slack_s_L->ElementWiseAbs(); tmp.v_L_NonConst()->Axpy(1., *tmp_slack_s_L); A_inf = std::max(A_inf, tmp.v_L_NonConst()->Amax()); @@ -1176,7 +1176,7 @@ Number PDFullSpaceSolver::NrmInf( Pd_U.ComputeColA1(*tmp.v_U_NonConst()); tmp.v_U_NonConst()->ElementWiseMultiply(v_U); tmp.v_U_NonConst()->ElementWiseAbs(); - auto tmp_slack_s_U = slack_s_U.MakeNewCopy(); + SmartPtr tmp_slack_s_U = slack_s_U.MakeNewCopy(); tmp_slack_s_U->ElementWiseAbs(); tmp.v_U_NonConst()->Axpy(1., *tmp_slack_s_U); A_inf = std::max(A_inf, tmp.v_U_NonConst()->Amax()); @@ -1240,7 +1240,7 @@ bool PDFullSpaceSolver::GMRES( { std::vector givens_c(restart), givens_s(restart), b_(restart+1), H(restart*(restart+1)); - std::vector> V, Z; + std::vector > V, Z; V.emplace_back(rhs.MakeNewIteratesVector()); V[0]->Set( 0. ); if ( !improve_solution ) @@ -1306,7 +1306,7 @@ bool PDFullSpaceSolver::GMRES( } for( Index k=0; k<=it; k++ ) { - auto tmp = V[it+1]->Dot( *V[k] ); + Number tmp = V[it+1]->Dot( *V[k] ); H[k+it*ldH] += tmp; V[it+1]->Axpy( -tmp, *V[k] ); } @@ -1317,11 +1317,11 @@ bool PDFullSpaceSolver::GMRES( } for( Index k = 1; k < it+1; k++ ) { - auto gamma = givens_c[k-1]*H[k-1+it*ldH] + givens_s[k-1]*H[k+it*ldH]; + Number gamma = givens_c[k-1]*H[k-1+it*ldH] + givens_s[k-1]*H[k+it*ldH]; H[k+it*ldH] = -givens_s[k-1]*H[k-1+it*ldH] + givens_c[k-1]*H[k+it*ldH]; H[k-1+it*ldH] = gamma; } - auto delta = std::sqrt( H[it+it*ldH]*H[it+it*ldH] + H[it+1+it*ldH]*H[it+1+it*ldH] ); + Number delta = std::sqrt( H[it+it*ldH]*H[it+it*ldH] + H[it+1+it*ldH]*H[it+1+it*ldH] ); givens_c[it] = H[it+it*ldH] / delta; givens_s[it] = H[it+1+it*ldH] / delta; H[it+it*ldH] = givens_c[it]*H[it+it*ldH] + givens_s[it]*H[it+1+it*ldH];