function x = m2fun(b, B, C, G, H, m, n)
% Solve the system L^T x = b.

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

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