function x = lsolve(L,y,W,NAA)
%LSOLVE Utility routine for "preconditioned" iterative methods.
%
% x = lsolve(L,y,W,NAA)
%
% Computes the vector
%    x = L_p*y
% where L_p is the A-weighted generalized inverse of L.
%
% Typically, L is a p-by-n band matrix with bandwidth n-p+1, W holds
% a basis for the null space of L, and NAA is a utility matrix which
% should be computed by routine pinit.
%
% Alternatively, L is square and dense, and W and NAA are not needed.
%
% Notice that x and y may be matrices, in which case
%    x(:,i) = L_p*y(:,i) .

% Reference: M. Hanke, "Regularization with differential operators.
% An iterative approach", J. Numer. Funct. Anal. Optim. 13 (1992),
% 523-540.

% Per Christian Hansen, UNI-C, and Martin Hanke, Institut fuer
% Praktische Mathematik, Universitaet Karlsruhe, 05/26/93.

% Initialization.
[p,n] = size(L); nu = n-p; [py,ly] = size(y);

% Special treatment of square L.
if (nu==0), x = L\y; return; end
 
% Compute a particular solution
x = [[eye(nu),zeros(nu,p)];L]\[zeros(nu,ly);y];

% Perform the necessary projection.
x = x - W*(NAA*x);
