<commit_before>#ifndef _MPRGPSOLVER_H_
#define _MPRGPSOLVER_H_

#include "MPRGPPrecondition.h"
#include "MPRGPProjection.h"

namespace MATH{

  /**
   * @class MPRGP a framework for the MPRGP method, with
   * preconditioning and projecting method can be changed through the
   * template parameters.
   * 
   * Solve such problem:
   * min_{x} 1/2*x^t*A*x-x^t*B s.t. n_i*x_j+p_i>= 0
   * 
   */
  template <typename T, typename MAT, typename PROJECTOIN,typename PRECONDITION>
  class MPRGP{

	typedef Eigen::Matrix<T,-1,1> Vec;
	
  public:
	MPRGP(const MAT &A,const Vec &B,
		  PRECONDITION &precond, PROJECTOIN &projector,
		  const int max_it = 1000, const T tol=1e-3):
	  _A(A), _B(B), _precond(precond), _projector(projector){

	  setParameters(tol,max_it);
	  _Gamma=1.0f;
	  _alphaBar=2.0f/specRad(_A,NULL,tol);
	  _iterationsOut = 0;
	  _residualOut = 0.0f;
	}
	
	int solve(Vec &result){

	  FUNC_TIMER();
	  
	  //declaration
	  T alphaCG,alphaF,beta;
	  Vec& AP=_gp;
	  Vec& AD=_gp;
	  Vec& y=_beta;	
	  Vec& xTmp=_beta;
	  Vec& D=_phi;

	  //initialization
	  _A.multiply(result,_g);
	  _g -= _B;
	  _projector.DECIDE_FACE(result);

	  assert_eq(_g,_g);
	  _projector.PHI(_g, _phi);
	  _precond.solve(_phi,_z);
	  _p = _z;
	  assert_eq(_p,_p);
	  int result_code = -1;
	  
	  int num_cg=0, num_exp=0, num_prop=0;

	  //MPRGP iteration
	  size_t iteration=0;
	  for(; iteration<_maxIterations; iteration++){

		DEBUG_LOG(setprecision(10)<<"func: "<<_A.funcValue(result,_B));

		//test termination
		assert_eq(_g,_g);
		_projector.PHI(_g,_phi);
		_projector.BETA(_g,_beta,_phi);
		assert_eq(_phi.size(), _beta.size());
		_gp = _phi+_beta;
		_residualOut=_gp.norm();

		assert_le(_phi.dot(_beta),ScalarUtil<T>::scalar_eps*_residualOut);
		DEBUG_LOG(setprecision(10)<<"||g||: "<<_g.norm());
		DEBUG_LOG(setprecision(10)<<"||beta||: "<<_beta.norm());
		DEBUG_LOG(setprecision(10)<<"||phi||: "<<_phi.norm());
		DEBUG_LOG(setprecision(10)<<"residual: "<<_residualOut);
		DEBUG_LOG("phi: "<<_phi.transpose());
		DEBUG_LOG("beta: "<<_beta.transpose());

		// debug
		DEBUG_FUN(assert(writeVTK(result, "beta_phi_g.vtk")));
		DEBUG_FUN(assert(printFace()));

		if(_residualOut <= _toleranceFactor){
		  _iterationsOut = iteration;
		  result_code = 0;
		  break;
		}

		//test proportional x: beta*beta <= gamma*gamma*phi*phiTilde
		const T beta_norm = _beta.norm();
		assert_eq(beta_norm, beta_norm);
		assert_eq(_g,_g);
		if(beta_norm*beta_norm <= _Gamma*_Gamma*_projector.PHITPHI(result,_alphaBar,_phi)){

		  //prepare conjugate gradient
		  _A.multiply(_p,AP);
		  const T pd = _p.dot(AP);
		  assert_eq(pd,pd);
		  assert_gt(pd,0); // pd = p^t*A*p > 0
		  alphaCG = (_z.dot(_g)) / pd;
		  assert_eq(alphaCG, alphaCG);
		  assert_ge(alphaCG,0.0f);
		  y = result-alphaCG*_p;
		  alphaF = _projector.stepLimit(result,_p,alphaCG);
		  assert_eq(alphaF, alphaF);
		  assert_ge(alphaF,0.0f);

		  DEBUG_LOG("alphaCG: "<<alphaCG);
		  DEBUG_LOG("alphaF: "<<alphaF);
		  DEBUG_LOG("alphaF-alphaCG: "<<alphaF-alphaCG);

		  if(alphaCG <= alphaF){
			//conjugate gradient step
			DEBUG_LOG("cg step");
			num_cg ++;
			assert_ge(alphaCG,0.0f);
			result = y;
			_g -= alphaCG*AP;

			assert_eq(_g,_g);
			_projector.PHI(_g, _phi);
			assert_ge(_g.dot(_phi),-ScalarUtil<T>::scalar_eps);
			_precond.solve(_phi,_z);
			assert_ge(_g.dot(_z),-ScalarUtil<T>::scalar_eps);
			beta = (_z.dot(AP)) / (_p.dot(AP));
			_p = _z-beta*_p;
			assert_eq(_p,_p);<commit_msg>removed a redundant variable y from the MPRGPSolver::solve() function.<commit_after>#ifndef _MPRGPSOLVER_H_
#define _MPRGPSOLVER_H_

#include "MPRGPPrecondition.h"
#include "MPRGPProjection.h"

namespace MATH{

  /**
   * @class MPRGP a framework for the MPRGP method, with
   * preconditioning and projecting method can be changed through the
   * template parameters.
   * 
   * Solve such problem:
   * min_{x} 1/2*x^t*A*x-x^t*B s.t. n_i*x_j+p_i>= 0
   * 
   */
  template <typename T, typename MAT, typename PROJECTOIN,typename PRECONDITION>
  class MPRGP{

	typedef Eigen::Matrix<T,-1,1> Vec;
	
  public:
	MPRGP(const MAT &A,const Vec &B,
		  PRECONDITION &precond, PROJECTOIN &projector,
		  const int max_it = 1000, const T tol=1e-3):
	  _A(A), _B(B), _precond(precond), _projector(projector){

	  setParameters(tol,max_it);
	  _Gamma=1.0f;
	  _alphaBar=2.0f/specRad(_A,NULL,tol);
	  _iterationsOut = 0;
	  _residualOut = 0.0f;
	}
	
	int solve(Vec &result){

	  FUNC_TIMER();
	  
	  //declaration
	  T alphaCG,alphaF,beta;
	  Vec& AP=_gp;
	  Vec& AD=_gp;
	  Vec& xTmp=_beta;
	  Vec& D=_phi;

	  //initialization
	  _A.multiply(result,_g);
	  _g -= _B;
	  _projector.DECIDE_FACE(result);

	  assert_eq(_g,_g);
	  _projector.PHI(_g, _phi);
	  _precond.solve(_phi,_z);
	  _p = _z;
	  assert_eq(_p,_p);
	  int result_code = -1;
	  
	  int num_cg=0, num_exp=0, num_prop=0;

	  //MPRGP iteration
	  size_t iteration=0;
	  for(; iteration<_maxIterations; iteration++){

		DEBUG_LOG(setprecision(10)<<"func: "<<_A.funcValue(result,_B));

		//test termination
		assert_eq(_g,_g);
		_projector.PHI(_g,_phi);
		_projector.BETA(_g,_beta,_phi);
		assert_eq(_phi.size(), _beta.size());
		_gp = _phi+_beta;
		_residualOut=_gp.norm();

		assert_le(_phi.dot(_beta),ScalarUtil<T>::scalar_eps*_residualOut);
		DEBUG_LOG(setprecision(10)<<"||g||: "<<_g.norm());
		DEBUG_LOG(setprecision(10)<<"||beta||: "<<_beta.norm());
		DEBUG_LOG(setprecision(10)<<"||phi||: "<<_phi.norm());
		DEBUG_LOG(setprecision(10)<<"residual: "<<_residualOut);
		DEBUG_LOG("phi: "<<_phi.transpose());
		DEBUG_LOG("beta: "<<_beta.transpose());

		// debug
		// DEBUG_FUN(assert(writeVTK(result, "beta_phi_g.vtk")));
		// DEBUG_FUN(assert(printFace()));

		if(_residualOut <= _toleranceFactor){
		  _iterationsOut = iteration;
		  result_code = 0;
		  break;
		}

		//test proportional x: beta*beta <= gamma*gamma*phi*phiTilde
		const T beta_norm = _beta.norm();
		assert_eq(beta_norm, beta_norm);
		assert_eq(_g,_g);
		if(beta_norm*beta_norm <= _Gamma*_Gamma*_projector.PHITPHI(result,_alphaBar,_phi)){

		  //prepare conjugate gradient
		  _A.multiply(_p,AP);
		  const T pd = _p.dot(AP);
		  assert_eq(pd,pd);
		  assert_gt(pd,0); // pd = p^t*A*p > 0
		  alphaCG = (_z.dot(_g)) / pd;
		  assert_eq(alphaCG, alphaCG);
		  assert_ge(alphaCG,0.0f);
		  alphaF = _projector.stepLimit(result,_p,alphaCG);
		  assert_eq(alphaF, alphaF);
		  assert_ge(alphaF,0.0f);

		  DEBUG_LOG("alphaCG: "<<alphaCG);
		  DEBUG_LOG("alphaF: "<<alphaF);
		  DEBUG_LOG("alphaF-alphaCG: "<<alphaF-alphaCG);

		  if(alphaCG <= alphaF){
			//conjugate gradient step
			DEBUG_LOG("cg step");
			num_cg ++;
			assert_ge(alphaCG,0.0f);
			result = result-alphaCG*_p;
			_g -= alphaCG*AP;

			assert_eq(_g,_g);
			_projector.PHI(_g, _phi);
			assert_ge(_g.dot(_phi),-ScalarUtil<T>::scalar_eps);
			_precond.solve(_phi,_z);
			assert_ge(_g.dot(_z),-ScalarUtil<T>::scalar_eps);
			beta = (_z.dot(AP)) / (_p.dot(AP));
			_p = _z-beta*_p;
			assert_eq(_p,_p);
{'fraction_non_alphanumeric': 0.1693548387096774, 'fraction_numerical': 0.012096774193548387, 'mean_word_length': 4.140636565507032, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5663533', 'n_tokens_mistral': 3205, 'n_tokens_neox': 2833, 'n_words': 461}