function [x, r, iter] = pregmres(A,b,x0,m,tol,maxiter,reorthog)
%pregmres Solve Ax = b using preconditioned GMRES
%
%   [x r iter] = pregmres(A,b,x0,m,tol,maxiter).
%   x0 is the initial approximation to the solution of the sparse
%   system Ax = b, d integer m < n is the dimension of the Krylov subspace
%   used. tol is the error tolerance, and maxiter is the maximum number of
%   iterations to perform. r is the final residual, and iter is the
%   number of iterations required. iter = -1 if the error tolerance
%   was not met. If this method fails, try the normal equations
%   with either cg or precg.

if nargin == 6
   reorthog = 1;
end

[L, U] = ilub(A);
b_bar = L\b;
B = L\(A/U);
[x_bar , ~, iter] = gmresb(B,b_bar,x0,m,tol,maxiter,reorthog);
x = backsolve(U,x_bar);
r = norm(A*x-b);
