kill(all)$

maxima:LU :decomposition A=L.U , version 1.0,part 1
peter.vlasschaert@gmail.com, 27/08/2020

load packages,option
=================

load(lrats)$ ratprint:false$

parameter
---------------

n:3; /* change n the example can't be used */

rem : n must be integer

LU decomposition
===============
used : part 2 'solve linear equations '    A.x = b , L.U.x = L.(U.x) = L.y = b

A:genmatrix(h,n,n);

p1:list_matrix_entries(A);

A:Matrix
L :Lower Triangle , with ones on diagonal.
U :Upper Triangle 

1e)  How construct  L,structure
========================

hl[i,j]:=if i > j then l[i,j] elseif i=j then 1 else 0;

for i:1 thru n do(
  for j:1 thru n do(
h[i,j]:hl[i,j]    
    )
)$

syntax : if condition
==============

if .cond1 then value1 elseif cond2 then value2 ..........elseif condi then valuei else valuef
rem : 1e) valuef = 0 in hl[i,j]

L:ev(A,nouns);

p2:list_matrix_entries(L);

p2a:listofvars(p2);

2e)  How construct  U,structure
========================

hu[i,j]:=if i <= j then u[i,j] else 0 ;

for i:1 thru n do(
  for j:1 thru n do(
h[i,j]:hu[i,j]    
    )
)$

U:ev(A,nouns);

p3:list_matrix_entries(U);

p3a:sort(p3);

a)  remove doubles,triples ( list  convert  set)
-------------------------------------------------------------

p3b:setify(p3a);

b) how to remove zeo 'p3b' (set  convert list)
-------------------------------------------------------------

p3c:listify(p3b);

c) start list b) from  2
-------------------------------------------------------------

p3d:makelist(p3c[i],i,2,length(p3c));

3e) How construct A
================
A = L.U

B:L.U;

4e) combination matrix(A)= matrix(B)
=============================

p4:A=B;

5e)  find  l[i,j] & u[i,j] as function f(h[i,j])
==============================

a) find all variables
---------------------------

p5a:append(p2a,p3d);

b)  solve f(h[i,j])     A  matrix replaced  C matrix => h by H
-----------------------------------------------------------------------------
reason :   h : definition 

C:genmatrix(H,n,n);

p5:C = B;

p5b:list_matrix_entries(C);

p5c:list_matrix_entries(B);

p5d:makelist(p5b[i]=p5c[i],i,1,n*n);

p5e:solve(p5d,p5a)[1];

final : L=> L1,U => U1  as function H[i,j]
==============================
use : 'p5e'

L :matrix
------------

L1:lratsubst(p5e,L);

U:matrix
------------

U1:lratsubst(p5e,U);

example :
-------------

D:matrix([1,1,-1],[1,2,-2],[-2,1,1]);

p6:[H[1,1]=1,H[1,2]=1,H[1,3]=-1,H[2,1]=1,H[2,2]=2,H[2,3]=-2,H[3,1]=-2,
    H[3,2]=1,H[3,3]=1];

L1 → L2, U1  U1→ U2
--------------------------------

p6a:L2:lratsubst(p6,L1);
p6b:U2:lratsubst(p6,U1);

p6c:L2.U2;

check : D against p6c
-------------------------------

is(D=p6c);

We can use this doc. generate for general  LU decomposition ,different 'integer =n'
( because also verified by the:example).