% analyze_SVD.m

%

% function [x,V_KA,iflag] = analyze_SVD(A,b,iprint);

%

% This MATLAB function uses SVD to analyze the existence of

% solutions to a linear problem with a (possibly) singular

% matrix.

%

% K. Beers

% MIT ChE

% 7/10/2002

function [x,V_KA,iflag] = analyze_SVD(A,b,iprint);

iflag = 0;

cutoff = 10*eps;

N = length(b);

if((size(A,1) ~= N) | (size(A,2) ~= N))

    iflag = -1;

    error('analyze_SVD: A is incorrectly dimensioned');

end

% First, check to see if the determinant of A is

% greater in magnitude than the cutoff value.  If so,

% simply return the unique solution using elimination

% methods.

det_A = det(A);

if(abs(det_A) >= cutoff)

    x = A\b;

    V_KA = 0;

    iflag = 1;

    return;

end

% Otherwise, use SVD analysis to "solve" the linear

% problem.

% First, perform a singular value decomposition.

[W,S,V] = svd(A);

if(iprint)

    W,

    S,

    V,

end

% Count the number of near-zero singular values

index_KA = find(diag(S) < cutoff);

dim_KA = length(index_KA);

if(iprint)

    index_KA,

    dim_KA,

end

% Check to see if b is in the range of A

b_in_RA = 1;

for k=1:dim_KA

    j = index_KA(k);

    var1 = dot(b,W(:,j));

    if(var1 > cutoff)

        b_in_RA = 0;

        break;

    end

end

if(iprint)

    b_in_RA,

end

% If b is not in the range of A, there is no solution

if(~b_in_RA)

    iflag = -1;

    x = 0;

    V_KA = 0;

    if(iprint)

        disp('analyze_SVD: b is not in the range of A');

    end

    return;

end

% If b is in the range of A, we form the particular

% solution of Ax = b that is of minimum length.

S_inv_corr = S;

for k=1:N

    if(S(k,k) >= cutoff)

        S_inv_corr(k,k) = 1/S(k,k);

    else

        S_inv_corr(k,k) = 0;

    end

end

x = V*S_inv_corr*W'*b;

if(iprint)

    S_inv_corr,

    x,

end

% We now return a column matrix whos columns are the basis

% vectors for the kernel of A.

V_KA = zeros(N,dim_KA);

for k=1:dim_KA

    j = index_KA(k);

    V_KA(:,k) = V(:,j);

end

if(iprint)

    V_KA,

end

iflag = 1;

return;