The Dooliitle Decomposition algoirithm as below: $$A=LU$$ where

$A= \begin{bmatrix} a_{11} & a_{12} & \cdots a_{1n}\\ a_{21} & a_{22} & \cdots a_{2n}\\ \vdots & \vdots & \vdots \\ a_{n1} & a_{n2} & \cdots a_{nn}\\ \end{bmatrix} L=\begin{bmatrix} 1 & 0 & \cdots 0\\ l_{21} & 1 & \cdots 0\\ \vdots & \vdots & \vdots \\ l_{n1} & l_{n2} & \cdots 1\\ \end{bmatrix} $ $ U=\begin{bmatrix} u_{11} & u_{12} & \cdots u_{1n}\\ 0 & a_{12} & \cdots u_{12}\\ \vdots & \vdots & \vdots \\ 0 & 0 & \cdots u_{nn}\\ \end{bmatrix} $

And

\begin{equation} \begin{cases} u_{kj}=a_{kj}-\sum_{m=1}^{k-1}l_{km}u_{mj} && (j=k,k+1,\cdots,n)\\ l_{ik}=1/u_{kk} (a_{ik}-\sum_{m=1}^{k-1}l_{im}u_{mk}) && (i=k+1,k+2,\cdots,n) \end{cases} \end{equation}

Here, $k=1,2,\cdots,n$

###Implementation

oneStep[mat_?MatrixQ, {intermediate_, index_}] := Module[{temp = intermediate, k = index, row}, row = Length@mat; Table[ temp[[k, j]] = mat[[k, j]] - Sum[temp[[k, m]] temp[[m, j]], {m, 1, k - 1}], {j, k, row}]; Table[ temp[[i, k]] = 1/temp[[k, k]] (mat[[i, k]] - Sum[temp[[i, m]] temp[[m, k]], {m, 1, k - 1}]), {i, k + 1, row}]; {temp, k+1} ] (*==============================*) doolittleDecomposite[mat_?MatrixQ] := Module[{intermediate, row}, intermediate = ConstantArray[0, Dimensions@mat]; row = Length@mat; First@ Nest[oneStep[mat, #] &, {intermediate, 1}, row] ] 

###Test

test={{1, 6, 1, 9}, {4, 4, 9, 7}, {10, 4, 10, 2}, {10, 3, 10, 2}}; doolittleDecomposite[test] 

$\left( \begin{array}{cccc} 1 & 6 & 1 & 9 \\ 4 & -20 & 5 & -29 \\ 10 & \frac{14}{5} & -14 & -\frac{34}{5} \\ 10 & \frac{57}{20} & \frac{57}{56} & \frac{11}{7} \\ \end{array} \right)$

Namely,

$L=\left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 4 & 1 & 0 & 0 \\ 10 & \frac{14}{5} & 1 & 0 \\ 10 & \frac{57}{20} & \frac{57}{56} & 1 \\ \end{array} \right) U=\left( \begin{array}{cccc} 1 & 6 & 1 & 9 \\ 0 & -20 & 5 & -29 \\ 0 & 0 & -14 & -\frac{34}{5} \\ 0 & 0 & 0 & \frac{11}{7} \\ \end{array} \right)$

###Performance testing

SeedRandom[1234]; testReal = RandomReal[{1., 100.}, {200, 200}]; doolittleDecomposite[testReal]; // AbsoluteTiming LUDecomposition[testReal]; // AbsoluteTiming 

{6.102349, Null} {0.034002, Null} (*big performance difference*) 

SeedRandom[12] testInteger = RandomInteger[{1, 100}, {100, 100}]; doolittleDecomposite[testInteger]; // AbsoluteTiming LUDecomposition[testInteger]; // AbsoluteTiming 

{6.173353, Null} {4.481256, Null} (*a bit difference*) 

###Question

• Is it possible to speed up my implementation code?
• Could someone give me some hints(like how to organize the data) to implement this algorithm elegantly in Mathematica?Namely, I think my method is not good:)