function [T,x_0] = pinit(W,A,b)
%PINIT Utility init.-procedure for "preconditioned" iterative methods.
% T = pinit(W,A)
% [T,x_0] = pinit(W,A,b)
% Initialization for `preconditioning' of general-form problems.
% Here, W holds a basis for the null space of L.
% Determines the matrix T needed in the iterative routines for
% treating regularization problems in general form.
% If b is also specified then x_0, the component of the solution in
% the null space of L, is also computed.
% Reference: P. C. Hansen, "Rank-Deficient and Discrete Ill-Posed Problems.
% Numerical Aspects of Linear Inversion", SIAM, Philadelphia, 1997.
% Per Christian Hansen, IMM, 07/29/97.
[n,nu] = size(W);
% Special treatment of square L.
if (nu==0), T = ; x_0 = zeros(n,1); return, end
% Compute T.
S = pinv(A*W);
T = S*A;
% If required, also compute x_0.
if (nargin==3), x_0 = W*(S*b); end