diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp index 195ab4dc0..062042c92 100644 --- a/src/Algorithm/IpPDFullSpaceSolver.cpp +++ b/src/Algorithm/IpPDFullSpaceSolver.cpp @@ -77,6 +77,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.", @@ -111,6 +119,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); @@ -138,6 +147,11 @@ bool PDFullSpaceSolver::Solve( DBG_ASSERT(!allow_inexact || !improve_solution); DBG_ASSERT(!improve_solution || beta == 0.); + if( gmres_refinement_ ) + { + return SolveGMRES( alpha, beta, rhs, res, allow_inexact, improve_solution ); + } + // Timing of PDSystem solver starts here IpData().TimingStats().PDSystemSolverTotal().Start(); @@ -233,7 +247,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; } @@ -241,11 +255,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, 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); + 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; @@ -263,9 +280,9 @@ 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); + residual_ratio = ComputeResidualRatio(rhs, res, *resid, A_nrm_inf); Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "residual_ratio = %e\n", residual_ratio); @@ -374,6 +391,250 @@ 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::SolveGMRES", 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 ) + { + + // 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); + + 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, A_nrm_inf); + 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, A_nrm_inf); + 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 +942,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 +969,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 +998,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 +1006,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 +1014,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 +1022,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); @@ -795,7 +1056,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); @@ -812,11 +1074,303 @@ 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); + /* 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 ); + } +} + +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. ); + 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()); + Pd_U.ComputeRowA1(*tmp.s_NonConst(), false); + 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() + Number(1.) + delta_d); + + 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(); + 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()); + tmp.z_U_NonConst()->ElementWiseMultiply(z_U); + tmp.z_U_NonConst()->ElementWiseAbs(); + 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()); + + 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(); + 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()); + tmp.v_U_NonConst()->ElementWiseMultiply(v_U); + tmp.v_U_NonConst()->ElementWiseAbs(); + 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()); + + if( norm_underestimated ) + { + Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, + "Cannot accurately compute matrix norm!\n"); + } + + return A_inf; +} + +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); + + // 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_ ) + { + restart = max_refinement_steps_; + } + Index totit = 0, ldH = restart + 1; + bool done = false; + 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); + 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 ) + { + 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, 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, res, *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_ ) || + NRBE_inf < std::numeric_limits::epsilon() || + rho < std::numeric_limits::epsilon() ) + { + resid.Copy( *V[0] ); + done = true; + break; + } + V[0]->Scal( 1./rho ); + b_[0] = rho; + 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., 0., res, *Z[it], *V[it+1]); + 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] ); + } + for( Index k=0; k<=it; k++ ) + { + Number 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(); + 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++ ) + { + 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; + } + 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]; + b_[it+1] = -givens_s[it]*b_[it]; + b_[it] = givens_c[it]*b_[it]; + + 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]; + } + 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, res, resid); + NRBE_inf = ComputeResidualRatio( rhs, res, 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_) || + NRBE_inf < std::numeric_limits::epsilon() || + (totit >= max_refinement_steps_) ) + { + done = true; + break; + } + } } + + return ( totit < max_refinement_steps_ ); } } // namespace Ipopt diff --git a/src/Algorithm/IpPDFullSpaceSolver.hpp b/src/Algorithm/IpPDFullSpaceSolver.hpp index 6e117821a..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 */ @@ -224,6 +228,62 @@ 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 + ); + + 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..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 { @@ -70,6 +72,42 @@ void SymScaledMatrix::ComputeRowAMaxImpl( THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SymScaledMatrix::ComputeRowAMaxImpl not implemented"); } +void SymScaledMatrix::ComputeRowA1Impl( + Vector& rows_norms, + bool init +) const +{ + 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( 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,