How to stop MPSolve iterations prematurely if any of converged roots have large magnitude

[advanpix]

This is not a bug or issue report, but cry for help.

I use MPSolve to find roots for a huge number of different polynomials - to select those polynomials with roots magnitude less than some threshold (|z| < T for all roots z).

To speed-up the search, I would like MPSolve to stop internal Ehrich-Aberth iterations as soon as at least one root |z| > T. So that such polynomials is quickly skipped.

I really need your help on how to implement this properly with MPSolve.

MPSolve stops iterations based on conditions are checked in mps_check_stop - it checks various combinations of root statuses and inclusion properties of the root.

Could you please advise what combinations of root status and inclusion indicate that root value is already reliable (won't change much) so that we can check if |z| < T?

I included following code to mps_check_stop :

for (i = 0; i < s->n; i++)
{
    if(MPS_ROOT_STATUS_IS_COMPUTED(s->root[i]->status))
    {
          if(s->lastphase == float_phase)
          {
                double r = hypot(cplx_Re(s->root[i]->fvalue),cplx_Im(s->root[i]->fvalue));
                if(r > threshold)
                {
                    return true; /* return 'true' to stop further iterations */
                }
          }
          else if(s->lastphase == dpe_phase)
          {
                double r = hypot(rdpe_get_d(cdpe_Re(s->root[i]->dvalue)),rdpe_get_d(cdpe_Im(s->root[i]->dvalue)));
                if(r > threshold)
                {
                    return true; /* return 'true' to stop further iterations */
                }
          }
          else if(s->lastphase == mp_phase)
          {
                double r = hypot(mpf_get_d(mpc_Re(s->root[i]->mvalue)),mpf_get_d(mpc_Im(s->root[i]->mvalue)));
                if(r > threshold)
                {
                    return true; /* return 'true' to stop further iterations */
                }
          }
    }
}

I use MPS_OUTPUT_GOAL_APPROXIMATE in mps_context_set_output_goal. Should I include any other conditions, e.g. s->root[i]->inclusion == MPS_ROOT_INCLUSION_OUT?

[robol]

I'll try to offer some advice, let me know if you manage to get this done. For what you want, you may safely ignore the s->root[i]->status variable. That only concern the meeting of the target to output an approximation, namely "Newton-isolation", which can be checked with MPS_ROOT_STATUS_IS_ISOLATED, and approsimation, which is (almost) never true until a very last stage, when the approximation is refined with Newton steps.

At any moment, you can check the roots and radii, and use those to exclude the roots. In your case you can check if abs(s->root[i]->fvalue) - s->root[i]->frad > T, and exit, without further checks. This needs to be done in the three versions you have above, because MPSolve might be working with the three different arithmetics, but I see you already got that part. If I remember correctly, the radii are set to infinite at the beginning, and then applied later on, but the inclusions should hold throughout the algorithm.

s->root[i]->inclusion is related to the case where the users only wants to compute roots inside some domain, and is used to track if those root are or not inside the domain of interest.

[advanpix]

Thank you for your quick help! I appreciate it very much.

As I see there are radii for double precision and for dpe, but no for mp:

struct mps_approximation {
  cplx_t fvalue;
  cdpe_t dvalue;
  mpc_t mvalue;

  double frad;
  rdpe_t drad;

  mps_boolean approximated;
  mps_boolean again;
  long int wp;
  ...
}

What radius we should use for the multiprecision phase?

[robol]

The radius for the multiprecision phase is also DPE; there is no need for the increased accuracy on the radius itself in higher precision, the only reason for switching to DPE is the larger exponent.

[advanpix]

Great, I will use DPE radius for multiprecision phase. Thank you very much for your quick help!