(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



LU[A_,0]:=Module[{LU=A, n=Length[A],d,e,r,i,j},
Do[
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]


LU[A_,1]:=Module[{LU=A, n=Length[A],d,e,r,i,j},
Do[
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
LU[[i,r+1;;n]]-= e* LU[[r,r+1;;n]],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]


LU[A_,2]:=Module[{LU=A, n=Length[A],d,e,r,i,j},
Do[
d=1/LU[[r,r]];
e = d LU[[r+1;;n,r]];
LU[[r+1;;n,r]]=e;
LU[[r+1;;n,r+1;;n]]-=Outer[Times,e,LU[[r,r+1;;n]]],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]


SparsePattern[AIn_SparseArray]:= Module[
{RulesForA},
RulesForA=ArrayRules[AIn];
SparseArray[
Append[Table[
RulesForA[[i,1]]-> 1,
{i,1, Length[RulesForA]-1}],
RulesForA[[-1]]
]
]
]


Clear[ILU,A]
ILU[A_,2]:=Module[
{LU=A, Patt,n=Length[A],d,e,r, Skeleton},
Patt=SparsePattern[A];
Do[
d=1/LU[[r,r]];
e = d LU[[r+1;;n,r]];
LU[[r+1;;n,r]]=e;
LU[[r+1;;n,r+1;;n]]-=Patt[[r+1;;n,r+1;;n]]*Outer[Times,e,LU[[r,r+1;;n]]];,
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]



