NOTAS:

LU

A = Table[Random[], {i, 7}, {j, 7}];

n = Length[A];
L = IdentityMatrix[n];   (*Doolitle*)
U = 
        L;
For[i = 1,i ≤ n,i++,
      U[[i, i]] = (A[[i, i]] - 
                Sum[L[[i, k]]*U[[k, i]], {k, 1, i - 1}]);  
      If[U[[i, i]] == 0, 
            Print["No admite descomp. LU"];  Break[]];
       For[j = i + 1, 
            j ≤ 
              n, j++, 
                  L[[j, i]] = (A[[j, i]] - 
                    Sum[L[[j, k]]*U[[k, i]], {k, 1, i - 1}])/
                U[[i, 
                    i]]; 
                  U[[i, j]] = (A[[i, j]] - 
                  Sum[L[[i, k]]*U[[k, j]], {k, 1, 
                      i - 1}]);
                 ]
      ]
L // 
      MatrixForm
U // 
      MatrixForm
A - 
        L . U//MatrixForm

(*ejemplo*)
b = 
          Table[Random[], {i, n}];
y = 
        LinearSolve[L, b];
x = 
        LinearSolve[U, y];
A . x - b

Sistema Triangular Superior

U = 
        Table[If[i ≤ j, Random[], 0], {i, 7}, {j, 
            7}];
b = Table[Random[], {i, 7}];

n = Length[U];
x = Table[0, {i, n}];
Do[x[[i]] = (b[[i]] - 
              Sum[U[[i, j]]*x[[j]], {j, i + 1, n}])/
          U[[i, i]], {i, n, 1, -1}]   
U . x - b // MatrixForm

Gauss Seidel

n = 7;
A = Table[Random[], {i, 7}, {j, 7}] + n*IdentityMatrix[n];   
b = 
        Table[Random[], {i, 7}];

(*primera forma,  
          pasos fijos sin criterio de parada*)
n = 
          Length[A];
x = 
        Table[0, {i, n}];
For[i = 1,i ≤ 10,i++,    
        Do[
            x[[k]] = (b[[k]] - 
                    Sum[A[[k, j]]*x[[j]], {j, 1, k - 1}] - 
                    Sum[A[[k, j]]*x[[j]], {j, k + 1, n}])/
                A[[k, k]], {k, 1, n}];
     ]
A . x - b

(*segunda forma,  con criterio de parada*)
n = 
          Length[A];
x = 
        Table[0, {i, n}];
error = 
        1;
contador = 
        0;
\(While[
      error > 10^\((\(-5\))\), 
\ \(contador++\); 
	y = x; 
	
      Do[x[\([k]\)] = \((b[\([k]\)] - 
                Sum[A[\([k, j]\)]*x[\([j]\)], {j, 1, k - 1}] - 
                Sum[A[\([k, j]\)]*y[\([j]\)], {j, k + 1, n}])\)/
            A[\([k, k]\)], {k, 1, n}]; 
	
      error = Max[
          Abs[y - x]];
]\)
Print[
      error,  ", no. pasos = ", 
      contador]
A . x - b