track()

int PathTracker::track

(

const Matrix *

start_sols

)

Definition at line 1640 of file NAG.cpp.

{
  double the_smallest_number = 1e-13;
  double epsilon2 = mpfr_get_d(epsilon, MPFR_RNDN);
  epsilon2 *= epsilon2;                           // epsilon^2
  double t_step = mpfr_get_d(init_dt, MPFR_RNDN);  // initial step
  double dt_min_dbl = mpfr_get_d(min_dt, MPFR_RNDN);
  double dt_increase_factor_dbl = mpfr_get_d(dt_increase_factor, MPFR_RNDN);
  double dt_decrease_factor_dbl = mpfr_get_d(dt_decrease_factor, MPFR_RNDN);
  double infinity_threshold2 = mpfr_get_d(infinity_threshold, MPFR_RNDN);
  infinity_threshold2 *= infinity_threshold2;
  double end_zone_factor_dbl = mpfr_get_d(end_zone_factor, MPFR_RNDN);
 
  if (C == nullptr)
    C = cast_to_CCC(
        start_sols->get_ring());  // fixes the problem for PrecookedSLPs
 
  int n = n_coords = start_sols->n_cols();
  n_sols = start_sols->n_rows();
 
  if (M2_numericalAlgebraicGeometryTrace > 1)
    printf(
        "epsilon2 = %e, t_step = %lf, dt_min_dbl = %lf, dt_increase_factor_dbl "
        "= %lf, dt_decrease_factor_dbl = %lf\n",
        epsilon2,
        t_step,
        dt_min_dbl,
        dt_increase_factor_dbl,
        dt_decrease_factor_dbl);
 
  // memory distribution for arrays
  complex* s_sols = newarray_atomic(complex, n * n_sols);
  raw_solutions = newarray(Solution, n_sols);
  complex* x0t0 = newarray_atomic(complex, n + 1);
  complex* x0 = x0t0;
  complex* t0 = x0t0 + n;
  complex* x1t1 = newarray_atomic(complex, n + 1);
  //  complex* x1 =  x1t1;
  //  complex* t1 = x1t1+n;
  complex* dxdt = newarray_atomic(complex, n + 1);
  complex* dx = dxdt;
  complex* dt = dxdt + n;
  complex* Hxt = newarray_atomic(complex, (n + 1) * n);
  complex* HxtH = newarray_atomic(complex, (n + 2) * n);
  complex* HxH = newarray_atomic(complex, (n + 1) * n);
  complex *LHS, *RHS;
  complex one_half(0.5, 0);
  complex* xt = newarray_atomic(complex, n + 1);
  complex* dx1 = newarray_atomic(complex, n);
  complex* dx2 = newarray_atomic(complex, n);
  complex* dx3 = newarray_atomic(complex, n);
  complex* dx4 = newarray_atomic(complex, n);
 
  // read solutions: rows are solutions
  int i, j;
  complex* c = s_sols;
  for (i = 0; i < n_sols; i++)
    for (j = 0; j < n; j++, c++)
      *c = complex(toBigComplex(C, start_sols->elem(i, j)));
 
  Solution* t_s = raw_solutions;  // current target solution
  complex* s_s = s_sols;          // current start solution
 
  for (int sol_n = 0; sol_n < n_sols; sol_n++, s_s += n, t_s++)
    {
      t_s->make(n, s_s);  // cook a Solution
      t_s->status = PROCESSING;
      bool end_zone = false;
      double tol2 =
          epsilon2;  // current tolerance squared, will change in end zone
      copy_complex_array<ComplexField>(n, s_s, x0);
      *t0 = complex(0, 0);
 
      *dt = complex(t_step);
      int predictor_successes = 0;
      int count = 0;  // number of steps
      while (t_s->status == PROCESSING &&
             1 - t0->getreal() > the_smallest_number)
        {
          if (1 - t0->getreal() <= end_zone_factor_dbl + the_smallest_number &&
              !end_zone)
            {
              end_zone = true;
              // to do: see if this path coincides with any other path on entry
              // to the end zone
            }
          if (end_zone)
            {
              if (dt->getreal() > 1 - t0->getreal())
                *dt = complex(1 - t0->getreal());
            }
          else
            {
              if (dt->getreal() > 1 - end_zone_factor_dbl - t0->getreal())
                *dt = complex(1 - end_zone_factor_dbl - t0->getreal());
            }
 
          bool LAPACK_success = false;
 
          if (is_projective)  // normalize
            multiply_complex_array_scalar<ComplexField>(
                n, x0, 1 / sqrt(norm2_complex_array<ComplexField>(n, x0)));
 
          // PREDICTOR in: x0t0,dt,pred_type
          //           out: dx
          switch (pred_type)
            {
              case TANGENT:
                {
                  evaluate_slpHxt(n, x0t0, Hxt);
                  LHS = Hxt;
                  RHS = Hxt + n * n;
                  multiply_complex_array_scalar<ComplexField>(n, RHS, -*dt);
                  LAPACK_success = solve_via_lapack_without_transposition(
                      n, LHS, 1, RHS, dx);
                }
                break;
              case EULER:
                {
                  evaluate_slpHxtH(n, x0t0, HxtH);  // for Euler "H" is attached
                  LHS = HxtH;
                  RHS = HxtH + n * (n + 1);  // H
                  complex* Ht = RHS - n;
                  multiply_complex_array_scalar<ComplexField>(n, Ht, *dt);
                  add_to_complex_array<ComplexField>(n, RHS, Ht);
                  negate_complex_array<ComplexField>(n, RHS);
                  LAPACK_success = solve_via_lapack_without_transposition(
                      n, LHS, 1, RHS, dx);
                }
                break;
              case RUNGE_KUTTA:
                {
                  copy_complex_array<ComplexField>(n + 1, x0t0, xt);
 
                  // dx1
                  evaluate_slpHxt(n, xt, Hxt);
                  LHS = Hxt;
                  RHS = Hxt + n * n;
                  //
                  negate_complex_array<ComplexField>(n, RHS);
                  solve_via_lapack_without_transposition(n, LHS, 1, RHS, dx1);
 
                  // dx2
                  multiply_complex_array_scalar<ComplexField>(
                      n, dx1, one_half * (*dt));
                  add_to_complex_array<ComplexField>(
                      n, xt, dx1);            // x0+.5dx1*dt
                  xt[n] += one_half * (*dt);  // t0+.5dt
                  //
                  evaluate_slpHxt(n, xt, Hxt);
                  LHS = Hxt;
                  RHS = Hxt + n * n;
                  //
                  negate_complex_array<ComplexField>(n, RHS);
                  LAPACK_success = solve_via_lapack_without_transposition(
                      n, LHS, 1, RHS, dx2);
 
                  // dx3
                  multiply_complex_array_scalar<ComplexField>(
                      n, dx2, one_half * (*dt));
                  copy_complex_array<ComplexField>(n, x0t0, xt);  // spare t
                  add_to_complex_array<ComplexField>(
                      n, xt, dx2);  // x0+.5dx2*dt
                  // xt[n] += one_half*(*dt); // t0+.5dt (SAME)
                  //
                  evaluate_slpHxt(n, xt, Hxt);
                  LHS = Hxt;
                  RHS = Hxt + n * n;
                  //
                  negate_complex_array<ComplexField>(n, RHS);
                  LAPACK_success =
                      LAPACK_success && solve_via_lapack_without_transposition(
                                            n, LHS, 1, RHS, dx3);
 
                  // dx4
                  multiply_complex_array_scalar<ComplexField>(n, dx3, *dt);
                  copy_complex_array<ComplexField>(n + 1, x0t0, xt);
                  add_to_complex_array<ComplexField>(n, xt, dx3);  // x0+dx3*dt
                  xt[n] += *dt;                                    // t0+dt
                  //
                  evaluate_slpHxt(n, xt, Hxt);
                  LHS = Hxt;
                  RHS = Hxt + n * n;
                  //
                  negate_complex_array<ComplexField>(n, RHS);
                  LAPACK_success =
                      LAPACK_success && solve_via_lapack_without_transposition(
                                            n, LHS, 1, RHS, dx4);
 
                  // "dx1" = .5*dx1*dt, "dx2" = .5*dx2*dt, "dx3" = dx3*dt
                  multiply_complex_array_scalar<ComplexField>(n, dx4, *dt);
                  multiply_complex_array_scalar<ComplexField>(n, dx1, 2);
                  multiply_complex_array_scalar<ComplexField>(n, dx2, 4);
                  multiply_complex_array_scalar<ComplexField>(n, dx3, 2);
                  add_to_complex_array<ComplexField>(n, dx4, dx1);
                  add_to_complex_array<ComplexField>(n, dx4, dx2);
                  add_to_complex_array<ComplexField>(n, dx4, dx3);
                  multiply_complex_array_scalar<ComplexField>(n, dx4, 1.0 / 6);
                  copy_complex_array<ComplexField>(n, dx4, dx);
                }
                break;
              case PROJECTIVE_NEWTON:
                {
                  evaluate_slpHxt(n, x0t0, Hxt);
                  LHS = Hxt;
                  RHS = Hxt + n * n;
                  LAPACK_success = solve_via_lapack_without_transposition(
                      n, LHS, 1, RHS, dx);  // dx used as temp
                  double chi2 = sqrt(norm2_complex_array<ComplexField>(n, RHS) +
                                     norm2_complex_array<ComplexField>(n, dx));
                  double chi1;
 
                  // evaluate again: LHS destroyed by solve_... !!!
                  evaluate_slpHxt(n, x0t0, Hxt);
                  LHS = Hxt;
 
                  for (j = 0; j < n; j++)
                    for (i = 0; i < n; i++)
                      *(LHS + n * i + j) =
                          *(LHS + n * i + j) *
                          DMforPN[j];  // multiply j-th column by sqrt(degree)
                  // print_complex_matrix(n,(double*)LHS); //!!!
                  norm_of_inverse_via_svd(n, LHS, chi1);
                  // printf("chi1=%lf\n", chi1);
                  *dt = 0.04804448 / (maxDegreeTo3halves * chi1 * chi2 * bigT);
                  if (dt->getreal() < dt_min_dbl) t_s->status = MIN_STEP_FAILED;
                  if (dt->getreal() > 1 - t0->getreal())
                    *dt = complex(1 - t0->getreal());
                  zero_complex_array<ComplexField>(n, dx);
                }
                break;
              default:
                ERROR("unknown predictor");
            };
 
          // make prediction
          copy_complex_array<ComplexField>(n + 1, x0t0, x1t1);
          add_to_complex_array<ComplexField>(n + 1, x1t1, dxdt);
 
          // CORRECTOR
          int n_corr_steps = 0;
          bool is_successful;
          do
            {
              n_corr_steps++;
              //
              evaluate_slpHxH(n, x1t1, HxH);
              LHS = HxH;
              RHS = HxH + n * n;  // i.e., H
              //
              negate_complex_array<ComplexField>(n, RHS);
              LAPACK_success =
                  LAPACK_success &&
                  solve_via_lapack_without_transposition(n, LHS, 1, RHS, dx);
              add_to_complex_array<ComplexField>(n, x1t1, dx);
              is_successful =
                  (pred_type == PROJECTIVE_NEWTON) ||
                  (norm2_complex_array<ComplexField>(n, dx) <
                   tol2 * norm2_complex_array<ComplexField>(n, x1t1));
              // printf("c: |dx|^2 = %lf\n",
              // norm2_complex_array<ComplexField>(n,dx));
            }
          while (!is_successful and n_corr_steps < max_corr_steps);
 
          if (!is_successful)
            {
              // predictor failure
              predictor_successes = 0;
              *dt = complex(dt_decrease_factor_dbl) * (*dt);
              if (dt->getreal() < dt_min_dbl) t_s->status = MIN_STEP_FAILED;
            }
          else
            {
              // predictor success
              predictor_successes = predictor_successes + 1;
              copy_complex_array<ComplexField>(n + 1, x1t1, x0t0);
              count++;
              if (is_successful &&
                  predictor_successes >= num_successes_before_increase)
                {
                  predictor_successes = 0;
                  *dt = complex(dt_increase_factor_dbl) * (*dt);
                }
            }
          if (norm2_complex_array<ComplexField>(n, x0) > infinity_threshold2)
            t_s->status = INFINITY_FAILED;
          if (!LAPACK_success) t_s->status = SINGULAR;
        }
      // record the solution
      copy_complex_array<ComplexField>(n, x0, t_s->x);
      t_s->t = t0->getreal();
      if (t_s->status == PROCESSING) t_s->status = REGULAR;
      evaluate_slpHxH(n, x0t0, HxH);
      cond_number_via_svd(n, HxH /*Hx*/, t_s->cond);
      t_s->num_steps = count;
      if (M2_numericalAlgebraicGeometryTrace > 0)
        {
          if (sol_n % 50 == 0) printf("\n");
          switch (t_s->status)
            {
              case REGULAR:
                printf(".");
                break;
              case INFINITY_FAILED:
                printf("I");
                break;
              case MIN_STEP_FAILED:
                printf("M");
                break;
              case SINGULAR:
                printf("S");
                break;
              default:
                printf("-");
            }
          fflush(stdout);
        }
    }
  if (M2_numericalAlgebraicGeometryTrace > 0) printf("\n");
 
  // clear arrays
  // freemem(t_sols); // do not delete (same as raw_solutions)
  freemem(s_sols);
  freemem(x0t0);
  freemem(x1t1);
  freemem(dxdt);
  freemem(xt);
  freemem(dx1);
  freemem(dx2);
  freemem(dx3);
  freemem(dx4);
  freemem(Hxt);
  freemem(HxtH);
  freemem(HxH);
 
  return n_sols;
}

References add_to_complex_array(), bigT, C, cast_to_CCC(), Solution::cond, cond_number_via_svd(), copy_complex_array(), DMforPN, dt_decrease_factor, dt_increase_factor, Matrix::elem(), end_zone_factor, epsilon, ERROR, EULER, evaluate_slpHxH(), evaluate_slpHxt(), evaluate_slpHxtH(), freemem(), Matrix::get_ring(), complex::getreal(), INFINITY_FAILED, infinity_threshold, init_dt, is_projective, M2_numericalAlgebraicGeometryTrace, Solution::make(), Matrix, max_corr_steps, maxDegreeTo3halves, min_dt, MIN_STEP_FAILED, multiply_complex_array_scalar(), Matrix::n_cols(), n_coords, Matrix::n_rows(), n_sols, negate_complex_array(), newarray, newarray_atomic, norm2_complex_array(), norm_of_inverse_via_svd(), Solution::num_steps, num_successes_before_increase, pred_type, PROCESSING, PROJECTIVE_NEWTON, raw_solutions, REGULAR, RUNGE_KUTTA, SINGULAR, solve_via_lapack_without_transposition(), Solution::status, Solution::t, TANGENT, toBigComplex(), Solution::x, and zero_complex_array().

Referenced by rawLaunchPT().