function [x_lambda,rho,eta] = tikhcstr(A,b,G,d,L,lambda,x_0)%TIKHCSTR Tikhonov regularization with linear inequality constraints.% % [x_lambda,rho,eta] = tikhcstr(A,b,G,d,L,lambda,x_0)% % Computes the constrained Tikhonov regularized solution x_lambda
% that solves the problem%    min { || A x - b ||^2 + lambda^2 || L (x - x_0) ||^2 }% subject to the linear inequality constraints
%    G x >= d . %% If x_0 is not specified, then x_0 = 0 is used.% If L is specified as the empty matrix [], then L = I is used.%% If lambda is a vector, then x_lambda is a matrix such that%   x_lambda = [ x_lambda(1), x_lambda(2), ... ] .%% The solution and residual norms are returned in eta and rho.
%
% If G and d define an empty set of constraints, then NaNs are returned.

% Ann-Charlotte Berglund, IMM & UNI-C, June 28, 1999.
% Initialization.[m,n] = size(A);if nargin < 6  error('Too few input arguments'),endif nargin < 7  x_0 = zeros(n,1);endif min(lambda) < 0  error('Illegal regularization parameter lambda'),endif isempty(L)  L = speye(n);end
% If no constraints then compute the ordinary Tikhonov solution.if isempty(G) | isempty(d)
  if norm(L-speye(n),'fro')==0
    [U,s,V] = csvd(A);    [x_lambda,rho,eta] = tikhonov(U,s,V,b,lambda,x_0);
  else
    [U,sm,X] = cgsvd(A,L);
    [x_lambda,rho,eta] = tikhonov(U,sm,X,b,lambda,x_0);
  end
  return;end
% Otherwise prepare for computing the constrainted solution.b = b(:);d = d(:);x_0 = x_0(:);[mG,nG] = size(G);% First try to find a feasible point.[p,status] = fpoint([],G,d,x_0);if ~strcmp(status,'feasible')  % There exists no feasible point for the given constraints.  x_lambda = NaN*ones(n,length(lambda));
  rho = NaN*ones(1,length(lambda));  eta = NaN*ones(1,length(lambda));  returnend% Constrained Tikhonov regularization.for i = 1:length(lambda)
    x_lambda(:,i) = tcsub([A;lambda(i)*L],[b;lambda(i)*L*x_0],...
                        G,d,x_0,0,'tikhcstr',mG,nG);  % Residual and solution norms.
  if nargout > 1, eta(i,1) = norm(L*x_lambda(:,i)); end
  if nargout == 3, rho(i,1) = norm(A*x_lambda(:,i)-b); end  
end % Loop over lambda-values.