%=======================================================================
% ch08_matrix_inverse.m
%
% Compute the inverse of a small matrix by Doolittle LU factorization.
%
% Origin:
%   MATLAB translation of ch08_matrix_inverse.f90 by Alexander Godunov,
%   based on Inverse.f90. The no-pivot LU factorization and repeated
%   triangular solves are retained.
%
% Example:
%   A = [3 2 4; 2 -3 1; 1 1 2].
%   Solve A*x=e_k for each column e_k of the identity, then place
%   those solutions into the columns of A^{-1}.
%
% Explicit inversion is seldom needed to solve A*x=b. If only x is
% needed, solve the system directly. The example shows how one LU
% factorization can be reused for several right-hand sides.
%=======================================================================

clear; clc;

A = [3.0,  2.0, 4.0;
     2.0, -3.0, 1.0;
     1.0,  1.0, 2.0];
n = size(A,1);

fprintf('Matrix inverse from Doolittle LU factorization\n');
fprintf('\nMatrix A:\n');
for i = 1:n
    fprintf('%12.6f%12.6f%12.6f\n', A(i,:));
end

% A_work is changed during elimination; keep A for verification.
[A_inverse, A_work, info] = inverse_lu(A);
if info ~= 0
    if info < 0
        error('Invalid matrix size or shape.');
    end
    error('Zero pivot at step %d. This no-pivot LU method cannot invert A.', info);
end

fprintf('\nInverse matrix A^{-1}:\n');
for i = 1:n
    fprintf('%12.6f%12.6f%12.6f\n', A_inverse(i,:));
end

% Multiply A and its computed inverse explicitly. Compare each entry
% with the corresponding entry of the identity matrix.
max_error = 0.0;
for i = 1:n
    for j = 1:n
        element = 0.0;
        for k = 1:n
            element = element + A(i,k)*A_inverse(k,j);
        end
        if i == j
            element = element - 1.0;
        end
        max_error = max(max_error, abs(element));
    end
end
fprintf('\nMaximum |A*A^{-1}-I| = %.4e\n', max_error);


function [C, A, info] = inverse_lu(A)
%=======================================================================
% Form C = A^{-1} by Doolittle LU and repeated triangular solves.
%
% Input:
%   A - square coefficient matrix, overwritten during elimination
%
% Output:
%   C    - inverse matrix if info = 0
%   A    - modified matrix after factorization
%   info - 0 on success; -1 for invalid size or shape;
%          positive factorization step of an exact zero pivot
%
% Factor A = L*U with unit diagonal in L. For each unit vector e_k,
% solve L*d=e_k and U*x=d, and store x in column k of C. The same
% LU factors are reused for each column.
%
% No row pivoting is performed, as in Inverse.f90. A nonsingular
% matrix may still have a zero pivot without row swaps. A small
% identity-check error alone does not establish an accurate inverse
% for an ill-conditioned matrix.
%=======================================================================

n = size(A,1);
C = zeros(n,n);
info = 0;
if n < 1 || size(A,2) ~= n
    info = -1;
    return
end
lower = zeros(n,n);
upper = zeros(n,n);

% Record the multipliers in L and update the trailing entries of A.
for k = 1:n-1
    if A(k,k) == 0.0
        info = k;
        return
    end
    for i = k+1:n
        coeff = A(i,k)/A(k,k);
        lower(i,k) = coeff;
        for j = k+1:n
            A(i,j) = A(i,j) - coeff*A(k,j);
        end
    end
end
if A(n,n) == 0.0
    info = n;
    return
end

% L has unit diagonal; U is the upper triangle after elimination.
for i = 1:n
    lower(i,i) = 1.0;
end
for j = 1:n
    for i = 1:j
        upper(i,j) = A(i,j);
    end
end

% Solve one system for each column of the identity matrix.
for k = 1:n
    rhs = zeros(n,1);
    rhs(k) = 1.0;
    d = zeros(n,1);
    x = zeros(n,1);

    % Forward substitution: L*d=e_k.
    for i = 1:n
        d(i) = rhs(i);
        for j = 1:i-1
            d(i) = d(i) - lower(i,j)*d(j);
        end
    end

    % Back substitution: U*x=d.
    for i = n:-1:1
        x(i) = d(i);
        for j = i+1:n
            x(i) = x(i) - upper(i,j)*x(j);
        end
        x(i) = x(i)/upper(i,i);
    end

    % The solution of A*x=e_k becomes column k of the inverse.
    for i = 1:n
        C(i,k) = x(i);
    end
end
end