function [A,CN]=svdinvert(B)
%
% Usage: [A] = svdinvert(B)
% This function preforms an SVD decomposition on the matrix B.  The
% inverse is returned in A.
%
% [A,CN] = svdinvert(B)
% The condition number is also calculated and returned in CN.
%

% Written by Jeffrey A Bloom


disp('   svdinvert')

MaxCondNum = 10e10;	% largest allowed condition number

[U,S,V]=svd(B);		% calculate the SVD

% B = U*S*V'  Where S is a diagonal matrix with diagonal elements s.
% A = B^(-1) = V*T*U'  Where T is a diagonal matrix with diagonal elements 1/s.

maxs = max(diag(S));    % largest eigenvalue
mins = min(diag(S));	% smallest eigenvalue
minw = maxs/MaxCondNum; % smallest allowed eigenvalue
smallw = find(diag(S)<minw); % find any zero eigenvalues

if(~isempty(smallw))
  disp('SVDINVERT WARNING: This matrix is too poorly conditioned for accurate results!')
  newvec = diag(S);	% get diag
  newvec(smallw) = inf*ones(size(smallw));  % replace zeros with infinity
  S = diag(newvec);	% replace S with new diag matrix
end

T = diag(1./diag(S));
A = V*T*U';

if nargout==2
  CN = maxs/mins;
end

return

% end function