diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/MILU0Preconditioner.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/MILU0Preconditioner.cs
new file mode 100644
index 00000000..bc6970f9
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Solvers/MILU0Preconditioner.cs
@@ -0,0 +1,267 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+using System;
+using MathNet.Numerics.LinearAlgebra.Solvers;
+using MathNet.Numerics.LinearAlgebra.Storage;
+using MathNet.Numerics.Properties;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
+{
+
+#if NOSYSNUMERICS
+ using Complex = Numerics.Complex;
+#else
+ using Complex = System.Numerics.Complex;
+#endif
+
+ ///
+ /// A simple milu(0) preconditioner.
+ ///
+ ///
+ /// Original Fortran code by Youcef Saad (07 January 2004)
+ ///
+ public sealed class MILU0Preconditioner : IPreconditioner
+ {
+ // Matrix stored in Modified Sparse Row (MSR) format containing the L and U
+ // factors together.
+
+ // The diagonal (stored in alu(0:n-1) ) is inverted. Each i-th row of the matrix
+ // contains the i-th row of L (excluding the diagonal entry = 1) followed by
+ // the i-th row of U.
+ private Complex[] _alu;
+
+ // The row pointers (stored in jlu(0:n) ) and column indices to off-diagonal elements.
+ private int[] _jlu;
+
+ // Pointer to the diagonal elements in MSR storage (for faster LU solving).
+ private int[] _diag;
+
+ /// Use modified or standard ILU(0)
+ public MILU0Preconditioner(bool modified = true)
+ {
+ UseModified = modified;
+ }
+
+ ///
+ /// Gets or sets a value indicating whether to use modified or standard ILU(0).
+ ///
+ public bool UseModified { get; set; }
+
+ ///
+ /// Gets a value indicating whether the preconditioner is initialized.
+ ///
+ public bool IsInitialized { get; private set; }
+
+ ///
+ /// Initializes the preconditioner and loads the internal data structures.
+ ///
+ /// The matrix upon which the preconditioner is based.
+ /// If is .
+ /// If is not a square or is not an
+ /// instance of SparseCompressedRowMatrixStorage.
+ public void Initialize(Matrix matrix)
+ {
+ var csr = matrix.Storage as SparseCompressedRowMatrixStorage;
+ if (csr == null)
+ {
+ throw new ArgumentException("Matrix must be in sparse storage format", "matrix");
+ }
+
+ // Dimension of matrix
+ int n = csr.RowCount;
+ if (n != csr.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix");
+ }
+
+ // Original matrix compressed sparse row storage.
+ Complex[] a = csr.Values;
+ int[] ja = csr.ColumnIndices;
+ int[] ia = csr.RowPointers;
+
+ _alu = new Complex[ia[n] + 1];
+ _jlu = new int[ia[n] + 1];
+ _diag = new int[n];
+
+ int code = Compute(n, a, ja, ia, _alu, _jlu, _diag, UseModified);
+ if (code > -1)
+ {
+ throw new Exception("Zero pivot encountered on row " + code + " during ILU process");
+ }
+
+ IsInitialized = true;
+ }
+
+ ///
+ /// Approximates the solution to the matrix equation Ax = b.
+ ///
+ /// The right hand side vector b.
+ /// The left hand side vector x.
+ public void Approximate(Vector input, Vector result)
+ {
+ if (_alu == null)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDoesNotExist);
+ }
+
+ if ((result.Count != input.Count) || (result.Count != _diag.Length))
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ int n = _diag.Length;
+
+ // Forward solve.
+ for (int i = 0; i < n; i++)
+ {
+ result[i] = input[i];
+ for (int k = _jlu[i]; k < _diag[i]; k++)
+ {
+ result[i] = result[i] - _alu[k] * result[_jlu[k]];
+ }
+ }
+
+ // Backward solve.
+ for (int i = n - 1; i >= 0; i--)
+ {
+ for (int k = _diag[i]; k < _jlu[i + 1]; k++)
+ {
+ result[i] = result[i] - _alu[k] * result[_jlu[k]];
+ }
+ result[i] = _alu[i] * result[i];
+ }
+ }
+
+ ///
+ /// MILU0 is a simple milu(0) preconditioner.
+ ///
+ /// Order of the matrix.
+ /// Matrix values in CSR format (input).
+ /// Column indices (input).
+ /// Row pointers (input).
+ /// Matrix values in MSR format (output).
+ /// Row pointers and column indices (output).
+ /// Pointer to diagonal elements (output).
+ /// Returns 0 on success or k > 0 if a zero pivot was encountered at step k.
+ private int Compute(int n, Complex[] a, int[] ja, int[] ia, Complex[] alu, int[] jlu, int[] ju, bool modified)
+ {
+ var iw = new int[n];
+ int i;
+
+ // Set initial pointer value.
+ int p = n + 1;
+ jlu[0] = p;
+
+ // Initialize work vector.
+ for (i = 0; i < n; i++)
+ {
+ iw[i] = -1;
+ }
+
+ // The main loop.
+ for (i = 0; i < n; i++)
+ {
+ int pold = p;
+
+ // Generating row i of L and U.
+ int j;
+ for (j = ia[i]; j < ia[i + 1]; j++)
+ {
+ // Copy row i of A, JA, IA into row i of ALU, JLU (LU matrix).
+ int jcol = ja[j];
+
+ if (jcol == i)
+ {
+ alu[i] = a[j];
+ iw[jcol] = i;
+ ju[i] = p;
+ }
+ else
+ {
+ alu[p] = a[j];
+ jlu[p] = ja[j];
+ iw[jcol] = p;
+ p = p + 1;
+ }
+ }
+
+ jlu[i + 1] = p;
+
+ Complex s = Complex.Zero;
+
+ int k;
+ for (j = pold; j < ju[i]; j++)
+ {
+ int jrow = jlu[j];
+ Complex tl = alu[j] * alu[jrow];
+ alu[j] = tl;
+
+ // Perform linear combination.
+ for (k = ju[jrow]; k < jlu[jrow + 1]; k++)
+ {
+ int jw = iw[jlu[k]];
+ if (jw != -1)
+ {
+ alu[jw] = alu[jw] - tl * alu[k];
+ }
+ else
+ {
+ // Accumulate fill-in values.
+ s = s + tl * alu[k];
+ }
+ }
+ }
+
+ if (modified)
+ {
+ alu[i] = alu[i] - s;
+ }
+
+ if (alu[i] == Complex.Zero)
+ {
+ return i;
+ }
+
+ // Invert and store diagonal element.
+ alu[i] = 1.0d / alu[i];
+
+ // Reset pointers in work array.
+ iw[i] = -1;
+ for (k = pold; k < p; k++)
+ {
+ iw[jlu[k]] = -1;
+ }
+ }
+
+ return -1;
+ }
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/MILU0Preconditioner.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/MILU0Preconditioner.cs
new file mode 100644
index 00000000..bfba1c97
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/MILU0Preconditioner.cs
@@ -0,0 +1,262 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+using System;
+using MathNet.Numerics.LinearAlgebra.Solvers;
+using MathNet.Numerics.LinearAlgebra.Storage;
+using MathNet.Numerics.Properties;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
+{
+ using Numerics;
+
+ ///
+ /// A simple milu(0) preconditioner.
+ ///
+ ///
+ /// Original Fortran code by Youcef Saad (07 January 2004)
+ ///
+ public sealed class MILU0Preconditioner : IPreconditioner
+ {
+ // Matrix stored in Modified Sparse Row (MSR) format containing the L and U
+ // factors together.
+
+ // The diagonal (stored in alu(0:n-1) ) is inverted. Each i-th row of the matrix
+ // contains the i-th row of L (excluding the diagonal entry = 1) followed by
+ // the i-th row of U.
+ private Complex32[] _alu;
+
+ // The row pointers (stored in jlu(0:n) ) and column indices to off-diagonal elements.
+ private int[] _jlu;
+
+ // Pointer to the diagonal elements in MSR storage (for faster LU solving).
+ private int[] _diag;
+
+ /// Use modified or standard ILU(0)
+ public MILU0Preconditioner(bool modified = true)
+ {
+ UseModified = modified;
+ }
+
+ ///
+ /// Gets or sets a value indicating whether to use modified or standard ILU(0).
+ ///
+ public bool UseModified { get; set; }
+
+ ///
+ /// Gets a value indicating whether the preconditioner is initialized.
+ ///
+ public bool IsInitialized { get; private set; }
+
+ ///
+ /// Initializes the preconditioner and loads the internal data structures.
+ ///
+ /// The matrix upon which the preconditioner is based.
+ /// If is .
+ /// If is not a square or is not an
+ /// instance of SparseCompressedRowMatrixStorage.
+ public void Initialize(Matrix matrix)
+ {
+ var csr = matrix.Storage as SparseCompressedRowMatrixStorage;
+ if (csr == null)
+ {
+ throw new ArgumentException("Matrix must be in sparse storage format", "matrix");
+ }
+
+ // Dimension of matrix
+ int n = csr.RowCount;
+ if (n != csr.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix");
+ }
+
+ // Original matrix compressed sparse row storage.
+ Complex32[] a = csr.Values;
+ int[] ja = csr.ColumnIndices;
+ int[] ia = csr.RowPointers;
+
+ _alu = new Complex32[ia[n] + 1];
+ _jlu = new int[ia[n] + 1];
+ _diag = new int[n];
+
+ int code = Compute(n, a, ja, ia, _alu, _jlu, _diag, UseModified);
+ if (code > -1)
+ {
+ throw new Exception("Zero pivot encountered on row " + code + " during ILU process");
+ }
+
+ IsInitialized = true;
+ }
+
+ ///
+ /// Approximates the solution to the matrix equation Ax = b.
+ ///
+ /// The right hand side vector b.
+ /// The left hand side vector x.
+ public void Approximate(Vector input, Vector result)
+ {
+ if (_alu == null)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDoesNotExist);
+ }
+
+ if ((result.Count != input.Count) || (result.Count != _diag.Length))
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ int n = _diag.Length;
+
+ // Forward solve.
+ for (int i = 0; i < n; i++)
+ {
+ result[i] = input[i];
+ for (int k = _jlu[i]; k < _diag[i]; k++)
+ {
+ result[i] = result[i] - _alu[k] * result[_jlu[k]];
+ }
+ }
+
+ // Backward solve.
+ for (int i = n - 1; i >= 0; i--)
+ {
+ for (int k = _diag[i]; k < _jlu[i + 1]; k++)
+ {
+ result[i] = result[i] - _alu[k] * result[_jlu[k]];
+ }
+ result[i] = _alu[i] * result[i];
+ }
+ }
+
+ ///
+ /// MILU0 is a simple milu(0) preconditioner.
+ ///
+ /// Order of the matrix.
+ /// Matrix values in CSR format (input).
+ /// Column indices (input).
+ /// Row pointers (input).
+ /// Matrix values in MSR format (output).
+ /// Row pointers and column indices (output).
+ /// Pointer to diagonal elements (output).
+ /// Returns 0 on success or k > 0 if a zero pivot was encountered at step k.
+ private int Compute(int n, Complex32[] a, int[] ja, int[] ia, Complex32[] alu, int[] jlu, int[] ju, bool modified)
+ {
+ var iw = new int[n];
+ int i;
+
+ // Set initial pointer value.
+ int p = n + 1;
+ jlu[0] = p;
+
+ // Initialize work vector.
+ for (i = 0; i < n; i++)
+ {
+ iw[i] = -1;
+ }
+
+ // The main loop.
+ for (i = 0; i < n; i++)
+ {
+ int pold = p;
+
+ // Generating row i of L and U.
+ int j;
+ for (j = ia[i]; j < ia[i + 1]; j++)
+ {
+ // Copy row i of A, JA, IA into row i of ALU, JLU (LU matrix).
+ int jcol = ja[j];
+
+ if (jcol == i)
+ {
+ alu[i] = a[j];
+ iw[jcol] = i;
+ ju[i] = p;
+ }
+ else
+ {
+ alu[p] = a[j];
+ jlu[p] = ja[j];
+ iw[jcol] = p;
+ p = p + 1;
+ }
+ }
+
+ jlu[i + 1] = p;
+
+ Complex32 s = Complex32.Zero;
+
+ int k;
+ for (j = pold; j < ju[i]; j++)
+ {
+ int jrow = jlu[j];
+ Complex32 tl = alu[j] * alu[jrow];
+ alu[j] = tl;
+
+ // Perform linear combination.
+ for (k = ju[jrow]; k < jlu[jrow + 1]; k++)
+ {
+ int jw = iw[jlu[k]];
+ if (jw != -1)
+ {
+ alu[jw] = alu[jw] - tl * alu[k];
+ }
+ else
+ {
+ // Accumulate fill-in values.
+ s = s + tl * alu[k];
+ }
+ }
+ }
+
+ if (modified)
+ {
+ alu[i] = alu[i] - s;
+ }
+
+ if (alu[i] == Complex32.Zero)
+ {
+ return i;
+ }
+
+ // Invert and store diagonal element.
+ alu[i] = 1.0f / alu[i];
+
+ // Reset pointers in work array.
+ iw[i] = -1;
+ for (k = pold; k < p; k++)
+ {
+ iw[jlu[k]] = -1;
+ }
+ }
+
+ return -1;
+ }
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/MILU0Preconditioner.cs b/src/Numerics/LinearAlgebra/Single/Solvers/MILU0Preconditioner.cs
new file mode 100644
index 00000000..c81d1bfd
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Single/Solvers/MILU0Preconditioner.cs
@@ -0,0 +1,260 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+using System;
+using MathNet.Numerics.LinearAlgebra.Solvers;
+using MathNet.Numerics.LinearAlgebra.Storage;
+using MathNet.Numerics.Properties;
+
+namespace MathNet.Numerics.LinearAlgebra.Single.Solvers
+{
+ ///
+ /// A simple milu(0) preconditioner.
+ ///
+ ///
+ /// Original Fortran code by Youcef Saad (07 January 2004)
+ ///
+ public sealed class MILU0Preconditioner : IPreconditioner
+ {
+ // Matrix stored in Modified Sparse Row (MSR) format containing the L and U
+ // factors together.
+
+ // The diagonal (stored in alu(0:n-1) ) is inverted. Each i-th row of the matrix
+ // contains the i-th row of L (excluding the diagonal entry = 1) followed by
+ // the i-th row of U.
+ private float[] _alu;
+
+ // The row pointers (stored in jlu(0:n) ) and column indices to off-diagonal elements.
+ private int[] _jlu;
+
+ // Pointer to the diagonal elements in MSR storage (for faster LU solving).
+ private int[] _diag;
+
+ /// Use modified or standard ILU(0)
+ public MILU0Preconditioner(bool modified = true)
+ {
+ UseModified = modified;
+ }
+
+ ///
+ /// Gets or sets a value indicating whether to use modified or standard ILU(0).
+ ///
+ public bool UseModified { get; set; }
+
+ ///
+ /// Gets a value indicating whether the preconditioner is initialized.
+ ///
+ public bool IsInitialized { get; private set; }
+
+ ///
+ /// Initializes the preconditioner and loads the internal data structures.
+ ///
+ /// The matrix upon which the preconditioner is based.
+ /// If is .
+ /// If is not a square or is not an
+ /// instance of SparseCompressedRowMatrixStorage.
+ public void Initialize(Matrix matrix)
+ {
+ var csr = matrix.Storage as SparseCompressedRowMatrixStorage;
+ if (csr == null)
+ {
+ throw new ArgumentException("Matrix must be in sparse storage format", "matrix");
+ }
+
+ // Dimension of matrix
+ int n = csr.RowCount;
+ if (n != csr.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix");
+ }
+
+ // Original matrix compressed sparse row storage.
+ float[] a = csr.Values;
+ int[] ja = csr.ColumnIndices;
+ int[] ia = csr.RowPointers;
+
+ _alu = new float[ia[n] + 1];
+ _jlu = new int[ia[n] + 1];
+ _diag = new int[n];
+
+ int code = Compute(n, a, ja, ia, _alu, _jlu, _diag, UseModified);
+ if (code > -1)
+ {
+ throw new Exception("Zero pivot encountered on row " + code + " during ILU process");
+ }
+
+ IsInitialized = true;
+ }
+
+ ///
+ /// Approximates the solution to the matrix equation Ax = b.
+ ///
+ /// The right hand side vector b.
+ /// The left hand side vector x.
+ public void Approximate(Vector input, Vector result)
+ {
+ if (_alu == null)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDoesNotExist);
+ }
+
+ if ((result.Count != input.Count) || (result.Count != _diag.Length))
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ int n = _diag.Length;
+
+ // Forward solve.
+ for (int i = 0; i < n; i++)
+ {
+ result[i] = input[i];
+ for (int k = _jlu[i]; k < _diag[i]; k++)
+ {
+ result[i] = result[i] - _alu[k] * result[_jlu[k]];
+ }
+ }
+
+ // Backward solve.
+ for (int i = n - 1; i >= 0; i--)
+ {
+ for (int k = _diag[i]; k < _jlu[i + 1]; k++)
+ {
+ result[i] = result[i] - _alu[k] * result[_jlu[k]];
+ }
+ result[i] = _alu[i] * result[i];
+ }
+ }
+
+ ///
+ /// MILU0 is a simple milu(0) preconditioner.
+ ///
+ /// Order of the matrix.
+ /// Matrix values in CSR format (input).
+ /// Column indices (input).
+ /// Row pointers (input).
+ /// Matrix values in MSR format (output).
+ /// Row pointers and column indices (output).
+ /// Pointer to diagonal elements (output).
+ /// Returns 0 on success or k > 0 if a zero pivot was encountered at step k.
+ private int Compute(int n, float[] a, int[] ja, int[] ia, float[] alu, int[] jlu, int[] ju, bool modified)
+ {
+ var iw = new int[n];
+ int i;
+
+ // Set initial pointer value.
+ int p = n + 1;
+ jlu[0] = p;
+
+ // Initialize work vector.
+ for (i = 0; i < n; i++)
+ {
+ iw[i] = -1;
+ }
+
+ // The main loop.
+ for (i = 0; i < n; i++)
+ {
+ int pold = p;
+
+ // Generating row i of L and U.
+ int j;
+ for (j = ia[i]; j < ia[i + 1]; j++)
+ {
+ // Copy row i of A, JA, IA into row i of ALU, JLU (LU matrix).
+ int jcol = ja[j];
+
+ if (jcol == i)
+ {
+ alu[i] = a[j];
+ iw[jcol] = i;
+ ju[i] = p;
+ }
+ else
+ {
+ alu[p] = a[j];
+ jlu[p] = ja[j];
+ iw[jcol] = p;
+ p = p + 1;
+ }
+ }
+
+ jlu[i + 1] = p;
+
+ float s = 0.0f;
+
+ int k;
+ for (j = pold; j < ju[i]; j++)
+ {
+ int jrow = jlu[j];
+ float tl = alu[j] * alu[jrow];
+ alu[j] = tl;
+
+ // Perform linear combination.
+ for (k = ju[jrow]; k < jlu[jrow + 1]; k++)
+ {
+ int jw = iw[jlu[k]];
+ if (jw != -1)
+ {
+ alu[jw] = alu[jw] - tl * alu[k];
+ }
+ else
+ {
+ // Accumulate fill-in values.
+ s = s + tl * alu[k];
+ }
+ }
+ }
+
+ if (modified)
+ {
+ alu[i] = alu[i] - s;
+ }
+
+ if (alu[i] == 0.0f)
+ {
+ return i;
+ }
+
+ // Invert and store diagonal element.
+ alu[i] = 1.0f / alu[i];
+
+ // Reset pointers in work array.
+ iw[i] = -1;
+ for (k = pold; k < p; k++)
+ {
+ iw[jlu[k]] = -1;
+ }
+ }
+
+ return -1;
+ }
+ }
+}
diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 870394e0..4a99efb3 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -122,8 +122,11 @@
+
+
+