function x = m1fun(b, B, C, G, H, m, n)
% Solve the system Ltil x = b.

n2 = n*n;
x = zeros(m*n2, 1);

for i = 1 : m
   ii = (i-1) * n + (1:n);
   ii2 = (i-1) * n2 + (1:n2);
   rhs = b(ii2);
   for j = 1 : i-1
      jj = (j-1) * n + (1:n);
      jj2 = (j-1) * n2 + (1:n2);
      y = kronUTSolveV(G(jj,:)', H(jj,:)', x(jj2));
      z = KPmult(B(ii,jj), C(ii,jj), y);
      rhs = rhs - z;
   end
   x(ii2) = kronLTSolve(G(ii,:), H(ii,:), rhs);
end
