Factorización LU

format long
fmt = format

warning: using the gnuplot graphics toolkit is discouraged

The gnuplot graphics toolkit is not actively maintained and has a number
of limitations that are unlikely to be fixed.  Communication with gnuplot
uses a one-directional pipe and limited information is passed back to the
Octave interpreter so most changes made interactively in the plot window
will not be reflected in the graphics properties managed by Octave.  For
example, if the plot window is closed with a mouse click, Octave will not
be notified and will not update its internal list of open figure windows.
The qt toolkit is recommended instead.

fmt = long

printf('Octave version: %s.', version)

Octave version: 9.1.0.

LUFactSimple.m

function [L, U, err] = LUFactSimple(A)
    % LU factorization of a square matrix.
    %
    % Input
    %   A(1:n,1:n) : the matrix to be factorized
    %
    % Output
    %   L(1:n,1:n), U(1:n,1:n) : the factors in A = LU, if Gaussian
    %                            elimination proceeds to completion
    %   err : = 0 if Gaussian elimination proceeds to completion
    %         = 1 if a zero pivot is encountered
    n = size(A, 1);
    U = A;
    L = eye(n);
    err = 0;

    for k = 1:n - 1

        for j = k + 1:n

            if abs(U(k, k)) > eps(U(k, k))
                L(j, k) = U(j, k) / U(k, k);
            else
                err = 1;
                return;
            end

            for t = k:n
                U(j, t) = U(j, t) - L(j, k) * U(k, t);
            end

        end

    end

end

Prueba de LUFactSimple.m

clearvars

n = 3;
A = rand(n, n)

tic
[L, U, err] = LUFactSimple(A)
toc

L * U

A =

   0.180410419500211   0.224341716470787   0.101449379297359
   0.681066884770393   0.271073279923097   0.747127494468375
   0.947437255029046   0.645668259768920   0.897827760160164

L =

   1.000000000000000                   0                   0
   3.775097284608863   1.000000000000000                   0
   5.251566165932764   0.924698652781662   1.000000000000000

U =

   0.180410419500211   0.224341716470787   0.101449379297359
                   0  -0.575838524750263   0.364146218157659
                   0                   0   0.028334114941347

err = 0

Elapsed time is 1.01985 seconds.

ans =

   0.180410419500211   0.224341716470787   0.101449379297359
   0.681066884770393   0.271073279923097   0.747127494468375
   0.947437255029046   0.645668259768920   0.897827760160164

TriDiagonal.m

function [x, err] = TriDiagonal(a, b, c, f)
    % Solution of a linear system of equations with a tridiagonal
    % matrix.
    %
    %   a(k) x(k-1) + b(k) x(k) + c(k) x(k+1) = f(k)
    %
    % with a(1) = c(n) = 0.
    %
    % Input
    %   a(1:n), b(1:n), c(1:n) : the coefficients of the system
    %                            matrix
    %   f(1:n) : the right hand side vector
    %
    % Output
    %   x(1:n) : the solution to the linear system of equations, if
    %            no division by zero occurs
    %   err : = 0, if no division by zero occurs
    %         = 1, if division by zero is encountered
    n = length(f);
    alpha = zeros(n, 1);
    beta = zeros(n, 1);
    err = 0;

    if abs(b(1)) > eps(b(1))
        alpha(1) = -c(1) / b(1);
        beta(1) = f(1) / b(1);
    else
        err = 1;
        return;
    end

    for k = 2:n
        denominator = a(k) * alpha(k - 1) + b(k);

        if abs(denominator) > eps(denominator)
            alpha(k) =- c(k) / denominator;
            beta(k) = (f(k) - a(k) * beta(k - 1)) / denominator;
        else
            err = 1;
            return;
        end

    end

    if abs(a(n) * alpha(n - 1) + b(n)) > eps(b(n))
        x(n) = (f(n) - a(n) * beta(n - 1)) / (a(n) * alpha(n - 1) + b(n));
    else
        err = 1;
        return;
    end

    for k = n - 1:-1:1
        x(k) = alpha(k) * x(k + 1) + beta(k);
    end

end

Prueba de TriDiagonal.m

clearvars

n = 5;
a = rand(n, 1)
b = rand(n, 1)
c = rand(n, 1)
f = rand(n, 1)

tic
[x, err] = TriDiagonal(a, b, c, f)
toc

lhs = @(k) a(k) * x(k - 1) + b(k) * x(k) + c(k) * x(k + 1);
rhs = @(k) f(k);

printf('k \t lhs(k) \t rhs(k)\n')
for k=2:n-1
    printf('k = %d\t%f\t%f\n', k, lhs(k), rhs(k))
end

a =

   5.190900154323299e-01
   7.015485663290333e-02
   7.680906530740392e-01
   2.296632816662104e-01
   8.611027785182878e-01

b =

   0.135650021055739
   0.385678086499345
   0.930896574152246
   0.429490435313398
   0.417994154478562

c =

   0.481433364216400
   0.369132687352825
   0.391021583338371
   0.254886212183741
   0.838640087272484

f =

   0.633109112276637
   0.449677053795114
   0.690511626181794
   0.418507575180630
   0.205226954598441

x =

 Columns 1 through 3:

  -5.507820233026685e+01   1.683403980794051e+01  -5.902592465957056e+00

 Columns 4 and 5:

  -1.724931355710016e+01   3.602600329407923e+01

err = 0

Elapsed time is 1.01978 seconds.

k 	 lhs(k) 	 rhs(k)

k = 2	0.449677	0.449677
k = 3	0.690512	0.690512
k = 4	0.418508	0.418508

CyclicallyTriDiagonal.m

function [x, err] = CyclicallyTriDiagonal(a, b, c, f)
    % Solution of a cyclically tridiagonal linear system of equations:
    %
    %   b(1) x(  1) + c(1) x(  2) + a(1) x(  n) = f(1),
    %   a(k) x(k-1) + b(k) x(  k) + c(k) x(k+1) = f(k),
    %   c(n) x(  1) + a(n) x(n-1) + b(n) x(  n) = f(n).
    %
    % Input
    %   a(1:n), b(1:n), c(1:n) : the coefficients of the system
    %                            matrix
    %   f(1:n) : the right hand side vector
    %
    % Output
    %   x(1:n) : the solution to the linear system of equations, if
    %            no division by zero is encountered
    %   err : = 0, if division by zero does not occur
    %         = 1, if division by zero occurs
    n = length(f);
    err = 0;
    newa = a(2:n);
    newa(1) = 0;
    newb = b(2:n);
    newc = c(2:n);
    newc(n - 1) = 0;
    newf = f(2:n);
    [u, err] = TriDiagonal(newa, newb, newc, newf);

    if err == 1
        return;
    end

    newf = zeros(n - 1, 1);
    newf(1) =- a(2);
    newf(n - 1) = -c(n);
    [v, err] = TriDiagonal(newa, newb, newc, newf);

    if err == 1
        return;
    end

    x = zeros(n, 1);
    denominator = b(1) + c(1) * v(1) + a(1) * v(n - 1);

    if abs(denominator) > eps(denominator)
        x(1) = (f(1) - c(1) * u(1) - a(1) * u(n - 1)) / denominator;
    else
        err = 1;
        return;
    end

    for k = 2:n
        x(k) = u(k - 1) + x(1) * v(k - 1);
    end

end

Prueba de CyclicallyTriDiagonal.m

clearvars

n = 5;
a = rand(n, 1)
b = rand(n, 1)
c = rand(n, 1)
f = rand(n, 1)

tic
[x, err] = CyclicallyTriDiagonal(a, b, c, f)
toc

lhs = @(k) a(k) * x(k - 1) + b(k) * x(k) + c(k) * x(k + 1);
rhs = @(k) f(k);

printf('k \t lhs(k) \t rhs(k)\n')
for k=2:n-1
    printf('k = %d\t%f\t%f\n', k, lhs(k), rhs(k))
end

a =

   0.994922875798857
   0.546348650244543
   0.177979837790046
   0.503483593329326
   0.816705345802433

b =

   7.866691507468093e-02
   4.012236119204076e-01
   2.661330928485295e-02
   1.775326498797364e-01
   3.655123349684792e-01

c =

   1.913542320563613e-01
   3.694490042524899e-01
   5.012223555645414e-01
   1.455424413559246e-01
   7.844040061393354e-02

f =

   0.481998773759005
   0.361029494105122
   0.646435105394834
   0.937804084591927
   0.889149543980786

x =

  -0.777796629719652
   0.671938132965005
   1.397703559295639
   0.976904048545877
   0.416723185119980

err = 0

Elapsed time is 1.02028 seconds.

k 	 lhs(k) 	 rhs(k)

k = 2	0.361029	0.361029
k = 3	0.646435	0.646435
k = 4	0.937804	0.937804

LUFactEfficient.m

function [Afact, swaps, err] = LUFactEfficient(A)
    % Efficient implementation of LU factorization with partial
    % pivoting.
    %
    % Input
    %   A(1:n,1:n) : the matrix to be factorized into A = P'LU
    %
    % Output
    %   Afact(1:n,1:n) : a matrix containing the matrices L and U
    %                    below and above the diagonal, respectively
    %   swaps : the indices that indicate the permutations
    %           (row swaps)
    %   err : = 0 if no the factorization finished successfully
    %         = 1 if a division by zero was encountered
    n = size(A, 1);
    swaps = 1:n;
    err = 0;

    for k = 1:n - 1
        i = k;

        for t = k:n

            if abs(A(t, k)) > abs(A(i, k))
                i = t;
            end

        end

        tt = swaps(k);
        swaps(k) = swaps(i);
        swaps(i) = tt;

        if abs(A(swaps(k), k)) <= eps(A(swaps(k), k))
            err = 1;
            return;
        end

        for i = k + 1:n
            xmult = A(swaps(i), k) / A(swaps(k), k);
            A(swaps(i), k) = xmult;

            for j = k + 1:n
                A(swaps(i), j) = A(swaps(i), j) - xmult * A(swaps(k), j);
            end

        end

    end

    Afact = A;
end

Prueba de LUFactEfficient.m

clearvars

n = 3;
A = rand(n, n)

tic
[Afact, swaps, err] = LUFactEfficient(A)
toc

A =

   5.354339928600693e-01   1.181666427160966e-02   6.556479993425937e-01
   2.134985312932789e-01   2.332058611025711e-02   1.599254811493648e-01
   4.649105358541434e-01   8.878330017306825e-01   6.909137354845939e-01

Afact =

   5.354339928600693e-01   1.181666427160966e-02   6.556479993425937e-01
   3.987392174203530e-01   2.120487311289858e-02  -1.040860872703211e-01
   8.682872997487171e-01   8.775727422182493e-01   1.216229045497647e-01

swaps =

   1   3   2

err = 0

Elapsed time is 1.02017 seconds.

CholeskyDecomposition.m

function [R, err] = CholeskyDecomposition(A)
    % Choleksy factorization of the HPD matrix A:
    %   A = R'*R
    %
    % Input:
    %   A(1:n,1:n) : an HPD matrix
    %
    % Output:
    %   R(1:n,1:n) : The upper triangular matrix satisfying A = R'*R
    %   err : = 0 if the decomposition finished successfully
    %         = 1 if an error occurred
    err = 0;
    n = size(A, 1);
    R = zeros(n, n);
    R(1, 1) = sqrt(A(1, 1));

    for i = 2:n

        for j = 1:i - 1
            sum = A(i, j);

            for k = 1:j - 1
                sum = sum - R(k, i) * conj(R(k, j));
            end

            if abs(R(j, j)) > eps(R(j, j))
                R(j, i) = sum / R(j, j);
            else
                err = 1;
                return;
            end

        end

        sum = A(i, i);

        for k = 1:i - 1
            sum = sum - R(k, i) * conj(R(k, i));
        end

        R(i, i) = sqrt(sum);
    end

end

Prueba de CholeskyDecomposition.m

clearvars

A = [1 1/2; 1/2 1]

tic
[R, err] = CholeskyDecomposition(A)
toc

R' * R

A =

   1.000000000000000   0.500000000000000
   0.500000000000000   1.000000000000000

R =

   1.000000000000000   0.500000000000000
                   0   0.866025403784439

err = 0

Elapsed time is 1.01975 seconds.

ans =

   1.000000000000000   0.500000000000000
   0.500000000000000   1.000000000000000

Solve.m

function [x, err] = Solve(A, f)
    % Solves the system Ax = f.
    %
    % Input
    %   A(1:n,1:n) : the coefficient matrix
    %   f(1:n) : the RHS vector
    %
    % Output:
    %   x(1:n) : the solution vector
    %   err : = 0 if the solution was found successfully
    %         = 1 if an error occurred during the process
    n = length(f);
    err = 0;
    [AA, swaps, err] = LUFactEfficient(A);

    if err == 1
        return
    end

    for k = 1:n - 1

        for i = k + 1:n
            f(swaps(i)) = f(swaps(i)) - AA(swaps(i), k) ...
                *f(swaps(k));
        end

    end

    if abs(AA(swaps(n), n)) > eps(AA(swaps(n), n))
        x(n) = f(swaps(n)) / AA(swaps(n), n);
    else
        err = 1;
        return;
    end

    for i = n - 1:-1:1
        xsum = f(swaps(i));

        for j = i + 1:n
            xsum = xsum - AA(swaps(i), j) * x(j);
        end

        if abs(AA(swaps(i), i)) > eps(AA(swaps(i), i))
            x(i) = xsum / AA(swaps(i), i);
        else
            err = 1;
            return;
        end

    end

end

Prueba de Solve.m

n = 10
A = rand(n, n);
f = rand(n, 1)

tic
[x, err] = Solve(A, f)
toc

A * x'

n = 10

f =

   3.917672960586684e-01
   9.370306322946349e-01
   1.262360299890020e-01
   8.792291497767564e-01
   7.539948459868598e-01
   1.572338968638454e-02
   9.100074826828775e-02
   6.306062946539864e-01
   9.160734040589031e-01
   9.394105707293133e-01

x =

 Columns 1 through 3:

   0.345304046709965   0.443695225859445  -0.507124860287545

 Columns 4 through 6:

   0.319498371625244  -0.461587550737034   0.452225343631443

 Columns 7 through 9:

   1.050964176714675  -0.208399561579305  -0.227501064451041

 Column 10:

  -0.358390153448041

err = 0

Elapsed time is 1.0227 seconds.

ans =

   3.917672960586683e-01
   9.370306322946350e-01
   1.262360299890022e-01
   8.792291497767565e-01
   7.539948459868596e-01
   1.572338968638456e-02
   9.100074826828769e-02
   6.306062946539864e-01
   9.160734040589033e-01
   9.394105707293132e-01