(* FATTORIZZAZIONE LU DI UNA MATRICE *)

Clear[a];n=4;
a={{1,-3,4,5},{2,3,4,-1},{5,-5,2,3},{-1,-2,7,6}};
MatrixForm[a]

a1=a;
Do[m[j][1]=a1[[j]][[1]]/a1[[1]][[1]],{j,2,n}];
a2={a1[[1]],a1[[2]]-m[2][1]*a1[[1]],a1[[3]]-m[3][1]*a1[[1]],a1[[4]]-m[4][1]*a1[[1]]};
Do[m[j][2]=a2[[j]][[2]]/a2[[2]][[2]],{j,3,n}];a3={a2[[1]],a2[[2]],a2[[3]]-m[3][2]*a2[[2]],a2[[4]]-m[4][2]*a2[[2]]};
Do[m[j][3]=a3[[j]][[3]]/a3[[3]][[3]],{j,4,n}]
a4={a3[[1]],a3[[2]],a3[[3]],a3[[4]]-m[4][3]*a3[[3]]};
MatrixForm[a2]
MatrixForm[a3]
MatrixForm[a4]
U=a4;

Do[ll[i][j]=0,{i,1,n},{j,1,n}];
Do[ll[i][i]=1,{i,1,n}];
Do[ll[i][j]=m[i][j],{i,2,n},{j,1,i-1}];

L=Table[ll[i][j],{i,1,n},{j,1,n}];
MatrixForm[L]

LU=Dot[L,U];
MatrixForm[LU]
MatrixForm[a]

(* FATTORIZZAZIONE CON SCAMBIO *)

Clear[a,a1,a2,a3,a4]
n=4;
a={{1,-3,4,5},{2,-6,4,-1},{5,-5,2,3},{-1,-2,7,6}};
MatrixForm[a]

a1=a;
Do[m[j][1]=a1[[j]][[1]]/a1[[1]][[1]],{j,2,n}];
a2={a1[[1]],a1[[2]]-m[2][1]*a1[[1]],a1[[3]]-m[3][1]*a1[[1]],a1[[4]]-m[4][1]*a1[[1]]};
MatrixForm[a2]

P={{1,0,0,0},{0,0,1,0},{0,1,0,0},{0,0,0,1}};
a1=Dot[P,a];
MatrixForm[a1]

Do[m[j][1]=a1[[j]][[1]]/a1[[1]][[1]],{j,2,n}]

a2={a1[[1]],a1[[2]]-m[2][1]*a1[[1]],a1[[3]]-m[3][1]*a1[[1]],a1[[4]]-m[4][1]*a1[[1]]};
Do[m[j][2]=a2[[j]][[2]]/a2[[2]][[2]],{j,3,n}];a3={a2[[1]],a2[[2]],a2[[3]]-m[3][2]*a2[[2]],a2[[4]]-m[4][2]*a2[[2]]};
Do[m[j][3]=a3[[j]][[3]]/a3[[3]][[3]],{j,4,n}]
a4={a3[[1]],a3[[2]],a3[[3]],a3[[4]]-m[4][3]*a3[[3]]};
MatrixForm[a2]
MatrixForm[a3]
MatrixForm[a4]
U=a4;

Do[ll[i][j]=0,{i,1,n},{j,1,n}];
Do[ll[i][i]=1,{i,1,n}];
Do[ll[i][j]=m[i][j],{i,2,n},{j,1,i-1}];

L=Table[ll[i][j],{i,1,n},{j,1,n}];
MatrixForm[L]

LU=Dot[L,U];
MatrixForm[LU]

MatrixForm[Dot[P,a]]

(* ESERCIZIO 6.5.2 *)

Clear[a,a1,a2,a3,a4]
n=3;
a={{1,2,-1},{2,4,0},{0,1,-1}};
a=a;
MatrixForm[a]

P={{1,0,0},{0,0,1},{0,1,0}};
MatrixForm[P]

a1=Dot[P,a];
MatrixForm[a1]

m[3][1]=a1[[3]][[1]]/a1[[1]][[1]];

a2={a1[[1]],a1[[2]],a1[[3]]-m[3][1]*a1[[1]]};
MatrixForm[a2]
U=a2;

Do[ll[i][j]=0,{i,1,n},{j,1,n}];
Do[ll[i][i]=1,{i,1,n}];
Do[ll[i][j]=m[i][j],{i,3,n},{j,1,i-1}];

L=Table[ll[i][j],{i,1,n},{j,1,n}];
MatrixForm[L]

LU=Dot[L,U];
MatrixForm[LU]
MatrixForm[a1]

(* ESERCIZIO 6.5.4 *)

Clear[a,a1,a2,a3,a4];
n=3;
a={{2,-1,1},{3,3,9},{3,3,5}};
MatrixForm[a]

a1=a;
Do[m[j][1]=a1[[j]][[1]]/a1[[1]][[1]],{j,2,n}];
a2={a1[[1]],a1[[2]]-m[2][1]*a1[[1]],a1[[3]]-m[3][1]*a1[[1]]};
Do[m[j][2]=a2[[j]][[2]]/a2[[2]][[2]],{j,3,n}];a3={a2[[1]],a2[[2]],a2[[3]]-m[3][2]*a2[[2]]};
MatrixForm[a2]
MatrixForm[a3]
U=a3;

Do[ll[i][j]=0,{i,1,n},{j,1,n}];
Do[ll[i][i]=1,{i,1,n}];
Do[ll[i][j]=m[i][j],{i,2,n},{j,1,i-1}];

L=Table[ll[i][j],{i,1,n},{j,1,n}];
MatrixForm[L]

LU=Dot[L,U];
MatrixForm[LU]
MatrixForm[a]

b={-1,0,4};
x={x1,x2,x3};
y={y1,y2,y3};
sys=Dot[L,y]==b;
Thread[sys]

y1sol=b[[1]]/L[[1]][[1]]
y2sol=(b[[2]]-L[[2]][[1]]*y1sol)/L[[2]][[2]]
y3sol=(b[[3]]-L[[3]][[1]]*y1sol-L[[3]][[2]]*y2sol)/L[[3]][[3]]

sol={y1->y1sol,y2->y2sol,y3->y3sol}

sys1=Dot[U,x]==y/.sol

x3sol=y3sol/U[[n]][[n]]
x2sol=(y2sol-U[[2]][[3]]*x3sol)/U[[n-1]][[n-1]]
x1sol=(y1sol-U[[1]][[3]]*x3sol-U[[1]][[2]]*x2sol)/U[[n-2]][[n-2]]

Solve[Dot[a,x]==b,x]