Schedule two computations together

[cordovan66]

Sorry for the long post but I'm a beginner and have been struggling for a few days.

I'm trying to schedule an algorithm so I can compute two quantities in a single pass through the input data. I have it working in two passes but want to make it faster. I pared my program down to only keep salient structures that may be important for the schedule. I'm computing G2 and P2 which are defined as the gaussian filtered versions of G * G and P * P respectively where G and P are input color images.
Gb is the repeat_edge version of G.
G2t is Gb*Gb
G2 is Gaussian filtered G2t
Similarly for P, Pb, P2t, P2
Gaussian filtering is done in a separable way.
I define Q=G2+8 and M=P2+3 and output = Q+M. These are simple placeholders for the actual computations in my program which is why I have not simply combined them. For example output=Q+M is a placeholder for a complicated operation that needs Q and M defined on identical domains.
The working schedules for Q and M parallelize, vectorize and unroll as shown. I've been trying for a few days to schedule P and Q in the same loop structure so it wouldn't have to go through the input data twice but nothing has worked. I tried
M.compute_at(Q,x) and M.compute_with(Q,x) and many versions of these but get "computed at invalid location" errors.
Here's the Generator code

class Cord : public Halide::Generator<Cord> {
public:
   Input<Buffer<int32_t, 3>> P{"P"};
   Input<Buffer<int32_t, 3>> G{ "G" };
   Output<Buffer<float, 3>> output{"output"};

   Cord():Q("Q"), M("M") {}

   void generate() {

       // Algorithm
   	// Gaussian 1D kernel
       float _ker[11] = { 0.0088,    0.0271,    0.0651,    0.1216,    0.1771,    0.2006,    0.1771,    0.1216,    0.0651,    0.0271,    0.0088 };
       Buffer<float, 1> ker(11);
       for (int i = 0; i < 11; i++) {
           ker(i) = _ker[i];
       }

       Func Gb("Gb"), G2("G2"), G2t("G2t"), Pb("Pb"), P2("P2"), P2t("P2t");

       Gb = BoundaryConditions::repeat_edge(G);
       G2t(x, y, c) = Gb(x, y, c) * Gb(x, y, c);
       G2 = filterGaussianSep(G2t, ker1, Q);
       
       Pb = BoundaryConditions::repeat_edge(P);
       P2t(x, y, c) = Pb(x, y, c) * Pb(x, y, c);
       P2 = filterGaussianSep(P2t, ker1, M);

       Q(x, y, c) = G2(x, y, c) + 8;
       M(x, y, c) = P2(x, y, c) + 3;

       output(x, y, c) = Q(x, y, c) + M(x, y, c);

       // Schedule
       if (!using_autoscheduler()) {
           int strip_size = 480;
           int vec = get_target().natural_vector_size(UInt(16));

           Q.bound(c, 0, 3);
           M.bound(c, 0, 3);

           G2t.compute_at(Q, yi).store_at(Q, yo).fold_storage(y, 11)
               .vectorize(x, vec);
           P2t.compute_at(M, yi).store_at(M, yo).fold_storage(y, 11)
               .vectorize(x, vec);

           Q
               .compute_root()
               .reorder(c, x, y)
               .unroll(c)
               .split(y, yo, yi, strip_size / 2)
               .vectorize(x, 2 * vec, TailStrategy::RoundUp)
               .parallel(yo);

           M
               //.compute_at(Q, x);
               //.compute_with(Q, x);
               .compute_root()
               .reorder(c, x, y)
               .unroll(c)
               .split(y, yo, yi, strip_size / 2)
               .vectorize(x, 2 * vec, TailStrategy::RoundUp)
               .parallel(yo);
       }
   }

private:
   Func filterGaussianSep(Func input, Buffer<float, 1> ker, Func parent);
   Func Q, M;
};

Func Cord::filterGaussianSep(Func input, Buffer<float, 1> ker, Func parent) {
   // Separable Gaussian filter
   Func filteredHor("filteredHor"), filtered("filtered");
   filteredHor(x, y, c) = cast<float>(0);
   RDom dxGauss(-5, 11, "dxGauss");
   filteredHor(x, y, c) += input(x + dxGauss, y, c) * ker(5 + dxGauss);

   filtered(x, y, c) = cast<float>(0);
   RDom dyGauss(-5, 11, "dyGauss");
   filtered(x, y, c) += filteredHor(x, y + dyGauss, c) * ker(5 + dyGauss);

   int vec = get_target().natural_vector_size(UInt(16));

   filteredHor
       .store_at(parent, yo)
       .fold_storage(y, 11)
       .compute_at(parent, yi)
       .vectorize(x, vec);

   return filtered;
}