% Computes pseudoinverse with Tikhonov-Phillips regulariztion for matrices whose
% conditioning can be deduced from their dimensions, i.e. m<n => purely underdetermined
% m>n => purely overdetermined,  m=n => nonsingular
% Inversion is applied on Gram-matrices (numerically not optimal, but faster (?),
% and, especially if only few rows of the pseudoinverse are demanded, memory saving
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      

% A: Matrix to be inverted (constraints see above)
% lambda (optional): Tikhonov-Phillips regularization parameter (Gram(A)+lambda*I is inverted)
% rows() (optional): Vector with rows of pseudoinverse G which shall be output

% OH 30.04.98




function  G = pinv_tikh(A, lambda, rows);

if nargin==0,
    disp(' G = pinv_tikh(A, lambda, rows); ');
    return;
end;


[m n] = size(A);

p = min(m,n);

if (nargin==1) | isempty(lambda), 
    lambda = 0.0; 
end;

if (m<n),
    [U,S,V] = svd(A', 0);		% SVD "economy size"
    MaxDiagVal=max(diag(S))
    S = S+eye(p)*(lambda.*MaxDiagVal);
    S = inv(S);
    G = V*S*U';
    if (nargin==1) | (nargin==2) | isempty(rows),		% if no special rows selected
        G = U*S*V';
	else
        G = U(rows,:)*S*V';
	end;
end;


if (m>n),
    [U,S,V] = svd(A, 0);		% SVD "economy size"
    MaxDiagVal=max(diag(S))
    S = S+eye(p)*(lambda.*MaxDiagVal);
    S = inv(S);
    G = V*S*U';
    if (nargin==1) | (nargin==2) | isempty(rows),		% if no special rows selected
        G = V*S*U';
	else
        G = V(rows,:)*S*U';
	end;
end;


if (m==n),
    [U,S,V] = svd(A, 0);		% SVD "economy size"
    MaxDiagVal=max(diag(S))    
    S = S+eye(p)*(lambda.*MaxDiagVal);
    S = inv(S);
    G = V*S*U';
    if (nargin==1) | (nargin==2) | isempty(rows),		% if no special rows selected
        G = V*S*U';
	else
        G = V(rows,:)*S*U';
	end;
end;
