function Y = mldivide(A,B)
% Embedded MATLAB Library function.
%
% Limitations:
% 1) No NaN & Inf inputs.
% 2) Can only handle square matrices A.

% $Revision: 1.1.6.1 $  $Date: 2004/12/16 21:55:22 $
% Copyright 2002-2004 The MathWorks, Inc.

eml_assert(nargin >= 2, 'error', 'Not enough input arguments.');
eml_assert(isfloat(A), 'error', ['Operation ''mldivide'' is not defined for values of class ''' class(A) '''.']);
eml_assert(size(A,1)==size(A,2),'Non-square matrices are not currently supported by mldivide.');
eml_assert(isscalar(A) || size(B,1)==size(A,2),'Matrix dimensions must agree.');

if isscalar(A)
    Y = B ./ A;
    return
end

if (isempty(B))
    Y = zeros(size(A,2),size(B,2));
    return
end

[X,pivot] = eml_lu(A);

if (~isreal(A) && isreal(B))
    Y = complex(B(pivot,:));
else
    Y = B(pivot,:);
end

n = size(A,2);
hasNotWarned = true;

for k = 1:size(B,2)
    % Solve L*Y = B
    for j = 1:n
        for i = j+1:n
            Y(i,k) = Y(i,k) - Y(j,k)*X(i,j);
        end
    end
    % Solve U*X = Y;
    for j = n:-1:1
        if hasNotWarned && X(j,j) == 0 
            warning('Matrix is singular to working precision.');
            hasNotWarned = false;
        end
        Y(j,k) = Y(j,k) / X(j,j);
        for i = 1:j-1
            Y(i,k) = Y(i,k) - Y(j,k)*X(i,j);
        end
    end
end

function [X,pivot] = eml_lu(A)
% Embedded MATLAB Library function.
%
% Limitations:
% 1) No NaN & Inf inputs.

% $INCLUDE(DOC) toolbox/eml/lib/matlab/lu.m $
% Copyright 2002-2004 The MathWorks, Inc.
% $Revision: 1.1.6.1 $  $Date: 2004/12/16 21:55:22 $

eml_assert(nargin > 0, 'error', 'Not enough input arguments.');
eml_assert(isfloat(A), 'error', ['Function ''lu'' is not defined for values of class ''' class(A) '''.']);

m = size(A,1);
n = size(A,2);

if isempty(A) || m == 1
    X = A;
    pivot = 1:m;
    return
end

% Init
X = A;

% Initialize row-pivot indices
pivot = 1:m;

for k = 1:min(m-1,n)
    p = k;    
    
    maxval = abs(X(k,k));
    % Scan the lower triangular part of this column only
    % MATLAB idiom:
    % [maxval,p] = max(abs(X(k:m,k)));
    % p = p+k-1;
    for i = k+1:m
        if(abs(X(i,k)) > maxval)
            maxval = abs(X(i,k));
            p = i;
        end
    end
    
    % Swap rows if max value if not in row k
    if(p ~= k)
        % MATLAB idiom:
        % X([p,k],:) = X([k,p],:); 
        for a = 1:n
            temp = X(p,a);
            X(p,a) = X(k,a);
            X(k,a) = temp;
        end
        % Swap pivot row indices
        % MATLAB idiom:
        % pivot([p,k]) = pivot([k,p]); 
        ptemp = pivot(k);
        pivot(k) = pivot(p);
        pivot(p) = ptemp;
    end

    % Divide lower triangular part of column by max
    Adiag = X(k,k);
    if Adiag ~= 0
        % MATLAB idiom:
        % X(k+1:m,k) = X(k+1:m,k) / X(k,k);
        for i = k+1:m
            X(i,k) = X(i,k) / Adiag;
        end
        % Subtract multiple of column from remaining columns
        % MATLAB idiom:
        % X(k+1:m,k+1:n) = X(k+1:m,k+1:n) - X(k+1:m,k)*X(k,k+1:n);
        for j = k+1:n
            for i = k+1:m
                X(i,j) = X(i,j) - X(i,k)*X(k,j);
            end
        end
    end
end