[Technical] tweak internal plan structure to simplify code

[mreineck]

Currently, we have a lot of the following data members in FINUFFT_PLAN_T:

  BIGINT ms;               // number of modes in x (1) dir (historical CMCL name) = N1
  BIGINT mt;               // number of modes in y (2) direction = N2
  BIGINT mu;               // number of modes in z (3) direction = N3

 [...]

  BIGINT nf1 = 1;          // size of internal fine grid in x (1) direction
  BIGINT nf2 = 1;          // " y (2)
  BIGINT nf3 = 1;          // " z (3)

 [...]

  std::vector<TF> phiHat1; // FT of kernel in t1,2, on x-axis mode grid
  std::vector<TF> phiHat2; // " y-axis.
  std::vector<TF> phiHat3; // " z-axis.

 [...]

  std::vector<TF> Xp, Yp, Zp; // internal primed NU points (x'_j, etc)
  std::vector<TF> Sp, Tp, Up; // internal primed targs (s'_k, etc)

This gives explicit names to variables referring dimensions 1, 2, 3. I suggest to combine the dimensions in the following way:

  std::array<BIGINT,3> m;    // number of modes in x, y, and z dirs (historical CMCL names N1, N2, N3)

 [...]

  std::array<BIGINT,3> nf = {1,1,1};          // size of internal fine grid in x, y, and z direction

 [...]

  std::array<std::vector<TF>,3> phiHat; // FT of kernel in t1,2, on xyz-axis mode grid

 [...]

  std::array<std::vector<TF>,3> XYZp; // internal primed NU points (x'_j etc)
  std::array<std::vector<TF>,3> STUp; // internal primed targs (s'_k, etc)

The reason is not (only) that this is more compact, but also because it enables many simplifications in the program code.

For example, this current snippet of setpts

    // pick x, s intervals & shifts & # fine grid pts (nf) in each dim...
    TF S1 = 0, S2 = 0, S3 = 0; // get half-width X, center C, which contains {x_j}...
    arraywidcen(nj, xj, &(t3P.X1), &(t3P.C1));
    arraywidcen(nk, s, &S1, &(t3P.D1)); // same D, S, but for {s_k}
    set_nhg_type3(S1, t3P.X1, opts, spopts, &(nf1), &(t3P.h1),
                  &(t3P.gam1));         // applies twist i)
    t3P.C2 = 0.0;                       // their defaults if dim 2 unused, etc
    t3P.D2 = 0.0;
    if (d > 1) {
      arraywidcen(nj, yj, &(t3P.X2), &(t3P.C2)); // {y_j}
      arraywidcen(nk, t, &S2, &(t3P.D2));        // {t_k}
      set_nhg_type3(S2, t3P.X2, opts, spopts, &(nf2), &(t3P.h2), &(t3P.gam2));
    }
    t3P.C3 = 0.0;
    t3P.D3 = 0.0;
    if (d > 2) {
      arraywidcen(nj, zj, &(t3P.X3), &(t3P.C3)); // {z_j}
      arraywidcen(nk, u, &S3, &(t3P.D3));        // {u_k}
      set_nhg_type3(S3, t3P.X3, opts, spopts, &(nf3), &(t3P.h3), &(t3P.gam3));
    }

could become something like

    // pick x, s intervals & shifts & # fine grid pts (nf) in each dim...
    std::array<TF,.3> S = {0,0,0}; // get half-width X, center C, which contains {x_j}...
    for (size_t idim=0; idim<dim; ++idim) {
      arraywidcen(nj, xyzj[idim], &(t3P.X[idim]), &(t3P.C[idim]));
      arraywidcen(nk, stuk[idim], &S[idim], &(t3P.D[idim])); // same D, S, but for {s_k}
      set_nhg_type3(S[idim], t3P.X[idim], opts, spopts, &(nf[idim]), &(t3P.h[idim]),
                  &(t3P.gam[idim]));         // applies twist i)
    }

There are quite a few locations where this kind of shortening is possible.

[ahbarnett]

I am all in favor of this. It's more abstract but shorter and more maintainable.

[mreineck]

See #592 for a possible implementation.

[DiamonDinoia]

Hi Martin,

How about using a multidimensional span together with this? This way it is possible to allocate the data in one go and reinterpret it as a multidimensions data structure at runtime, data can be 1,2 or 3 dimensional.

Thanks,
Marco

[mreineck]

Spans should work nicely here! But do we already have them available in our current set of dependencies? I didn't think it's worth pulling in another external library for this purpose.
(ducc's mav structures could be useful here; but that would make ducc a mandatory dependency, which again may not be desirable.)

[DiamonDinoia]

I have nothing against meaningful dependencies. How about others?

[mreineck]

Hmm, if spans are introduced, I'd prefer if it is not me doing it. I wouldn't know which external package to choose ...

[ahbarnett]

I feel that this is a minor clean-up, so we shouldn't introduce a new dependency. Given the pain of dependencies, we should only do it if there is a large performance gain. This PR and discussion is more just about tidying some mem mgmt. Anyway, best wishes to you both, Alex