kill(all)$

maxima:myLU :decomposition A=L.U ( = function) , version 1.0,part 3
peter.vlasschaert@gmail.com, 28/08/2020

see website:maxima:LU :decomposition A=L.U , version 1.0,part 1,2
rem : very good example programming in  Maxima.( square matrix)

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

load(lrats)$ ratprint:false$

A,A1,..,Ai = matrix representation
==========================

v:4;

m[i,j]:= if i=j then f[i,j] else 1;
Ai : genmatrix(m,v,v); 
p1a:list_matrix_entries(%); /*math*/
p1b:sort(p1a);              
p1c:listofvars(p1b);        /*find variables*/

How to make function : myLU[A]
=========================

ref :maxima:LU :decomposition A=L.U , version 1.0,part 1

hl[i,j]:=if i > j then l[i,j] elseif i=j then 1 else 0;      ( L :entries)
hu[i,j]:=if i <= j then u[i,j] else 0 ;                           ( U: entries)

ref :maxima:LU :decomposition A=L.U , version 1.0,part 2

A:matrix(
		[1,	1,	2],
		[2,	-1,	2],
		[4,	1,	4]
	);

first : try ( theoretical : LU function)

myLUtheo(A):= block([v,L,U],local(hl,hu),
v:length(A),
hl[i,j]:=if i > j then l[i,j] elseif i=j then 1 else 0,
L: genmatrix(hl,v,v),   
hu[i,j]:=if i <= j then u[i,j] else 0,
U: genmatrix(hu,v,v),
    [L,U]    
)$

optimized the code
---------------------------

grind(myLUtheo);

p2a:myLUtheo(A);

LL:p2a[1];
UU:p2a[2];

How get :matrix entrie ' first matrix'

p2b:p2a[1][2,1];

p2b1:list_matrix_entries(LL.UU-A);

p2b2:solve(%);

p2b3:flatten(p2b2);

LL  ( = lower triangle of LU decomposition)
-----------------------------------------------------------

p2b4a:lratsubst(p2b3,LL);

UU  ( = upper triangle of LU decomposition)
-----------------------------------------------------------

p2b4b:lratsubst(p2b3,UU);

check : L.U-A  = 0 ?
---------------------------

p2b4c:p2b4a.p2b4b-A;

info :myLUtheo(A)
-------------------------

p2c1:[L,U];  /* only known inside the function */

p2c2:hl[2,1]; /*only known inside the function */

second :try ( practical: LU function)
------------------------------------------------
use : principal 'p2b1,p2b2,p2b4a,p2b4b'

myLU(A):= block([v,L,U],local(hl,hu),
/* theoretical (first part) */    
v:length(A),
hl[i,j]:=if i > j then l[i,j] elseif i=j then 1 else 0,
L: genmatrix(hl,v,v),   
hu[i,j]:=if i <= j then u[i,j] else 0,
U: genmatrix(hu,v,v),
/******************************************************/    
/* practical   (second part) */
q1:list_matrix_entries(L.U-A),/*equation :rhs(not aviable)*/
q2:flatten(solve(q1)),        /* input : use lratsubst    */
q3:lratsubst(q2,[L,U])
/******************************************************/      
)$

optimized code.

grind(myLU);

result :myLU(A)
============

myLU(A);

 L = "Lower triangle ( LU decomposition)"
================================

L:myLU(A)[1];

U="Upper triangle ( LU decomposition)"
===============================

U:myLU(A)[2];

info :
===
ref :maxima:LU :decomposition A=L.U , version 1.0,part 2

How to do this with : Maxima 'commands'
---------------------------------------------------------

lu_factor(A,generalring)$
[DT,LT,UT]:get_lu_factors(%);

1e)  A = LT.D.UT  
=============

LT.DT.UT;

Matrix decomposition :wikipedia
https://en.wikipedia.org/wiki/Matrix_decomposition
A=LDU , A=LUP  ( P= permutation matrix)