adding residual_improvement_factor_; remove use of auto

pghysels · Jul 17, 2024 · cafb7e8 · cafb7e8
1 parent 530ae0a
commit cafb7e8
Showing 1 changed file with 18 additions and 11 deletions.
diff --git a/src/Algorithm/IpPDFullSpaceSolver.cpp b/src/Algorithm/IpPDFullSpaceSolver.cpp
@@ -1152,31 +1152,31 @@ 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<Vector> 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();
- auto tmp_slack_x_U = slack_x_U.MakeNewCopy();
+ SmartPtr<Vector> 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();
- auto tmp_slack_s_L = slack_s_L.MakeNewCopy();
+ SmartPtr<Vector> 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();
- auto tmp_slack_s_U = slack_s_U.MakeNewCopy();
+ SmartPtr<Vector> 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<Number>::max();
+ Number NRBE_inf = std::numeric_limits<Number>::max();
  while ( !done )
  {
  std::vector<Number> givens_c(restart), givens_s(restart),
  b_(restart+1), H(restart*(restart+1));
- std::vector<SmartPtr<IteratesVector>> V, Z;
+ std::vector<SmartPtr<IteratesVector> > 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<Number>::epsilon() )
+ rho < std::numeric_limits<Number>::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;