diff --git a/.idea/.idea.TrainAR/.idea/workspace.xml b/.idea/.idea.TrainAR/.idea/workspace.xml
index 3944e2e9..23b98578 100644
--- a/.idea/.idea.TrainAR/.idea/workspace.xml
+++ b/.idea/.idea.TrainAR/.idea/workspace.xml
@@ -10,6 +10,9 @@
+
+
+
@@ -101,6 +104,7 @@
+
diff --git a/Assembly-CSharp-Editor.csproj b/Assembly-CSharp-Editor.csproj
index 8e01de22..2133fe27 100644
--- a/Assembly-CSharp-Editor.csproj
+++ b/Assembly-CSharp-Editor.csproj
@@ -865,6 +865,10 @@
{8e530897-0605-6cd4-3bfe-114dbb14bed6}Whinarn.UnityMeshSimplifier.Runtime
+
+ {ff750720-481d-f431-2cbf-6b672d06727f}
+ UnityMeshDecimation
+
restart with the original vector
+ // - calculation diverged --> restart with the original vector
+ // - Some unknown status occurred --> To be safe restart.
+ input.CopyTo(internalInput);
+ }
+ }
+ }
+ }
+}
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs.meta
new file mode 100644
index 00000000..1bef5c6d
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: 1fd7d5b601567f34a9ceab68e52b34f1
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/DiagonalPreconditioner.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/DiagonalPreconditioner.cs
new file mode 100644
index 00000000..035880d3
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/DiagonalPreconditioner.cs
@@ -0,0 +1,107 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2010 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;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// A diagonal preconditioner. The preconditioner uses the inverse
+ /// of the matrix diagonal as preconditioning values.
+ ///
+ public sealed class DiagonalPreconditioner : IPreconditioner
+ {
+ ///
+ /// The inverse of the matrix diagonal.
+ ///
+ Complex[] _inverseDiagonals;
+
+ ///
+ /// Returns the decomposed matrix diagonal.
+ ///
+ /// The matrix diagonal.
+ internal DiagonalMatrix DiagonalEntries()
+ {
+ var result = new DiagonalMatrix(_inverseDiagonals.Length);
+ for (var i = 0; i < _inverseDiagonals.Length; i++)
+ {
+ result.At(i, i, 1/_inverseDiagonals[i]);
+ }
+
+ return result;
+ }
+
+ ///
+ /// Initializes the preconditioner and loads the internal data structures.
+ ///
+ ///
+ /// The upon which this preconditioner is based.
+ /// If is .
+ /// If is not a square matrix.
+ public void Initialize(Matrix matrix)
+ {
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException("Matrix must be square.", nameof(matrix));
+ }
+
+ _inverseDiagonals = new Complex[matrix.RowCount];
+ for (var i = 0; i < matrix.RowCount; i++)
+ {
+ _inverseDiagonals[i] = 1/matrix.At(i, i);
+ }
+ }
+
+ ///
+ /// Approximates the solution to the matrix equation Ax = b.
+ ///
+ /// The right hand side vector.
+ /// The left hand side vector. Also known as the result vector.
+ public void Approximate(Vector rhs, Vector lhs)
+ {
+ if (_inverseDiagonals == null)
+ {
+ throw new ArgumentException("The requested matrix does not exist.");
+ }
+
+ if ((lhs.Count != rhs.Count) || (lhs.Count != _inverseDiagonals.Length))
+ {
+ throw new ArgumentException("All vectors must have the same dimensionality.", nameof(rhs));
+ }
+
+ for (var i = 0; i < _inverseDiagonals.Length; i++)
+ {
+ lhs[i] = rhs[i]*_inverseDiagonals[i];
+ }
+ }
+ }
+}
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/DiagonalPreconditioner.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/DiagonalPreconditioner.cs.meta
new file mode 100644
index 00000000..b0ccde44
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/DiagonalPreconditioner.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: cbd21724cfc5ab84f90b4e7af3fd4cf9
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs
new file mode 100644
index 00000000..eb870acb
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs
@@ -0,0 +1,373 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// 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;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// A Generalized Product Bi-Conjugate Gradient iterative matrix solver.
+ ///
+ ///
+ ///
+ /// The Generalized Product Bi-Conjugate Gradient (GPBiCG) solver is an
+ /// alternative version of the Bi-Conjugate Gradient stabilized (CG) solver.
+ /// Unlike the CG solver the GPBiCG solver can be used on
+ /// non-symmetric matrices.
+ /// Note that much of the success of the solver depends on the selection of the
+ /// proper preconditioner.
+ ///
+ ///
+ /// The GPBiCG algorithm was taken from:
+ /// GPBiCG(m,l): A hybrid of BiCGSTAB and GPBiCG methods with
+ /// efficiency and robustness
+ ///
+ /// S. Fujino
+ ///
+ /// Applied Numerical Mathematics, Volume 41, 2002, pp 107 - 117
+ ///
+ ///
+ ///
+ /// The example code below provides an indication of the possible use of the
+ /// solver.
+ ///
+ ///
+ public sealed class GpBiCg : IIterativeSolver
+ {
+ ///
+ /// Indicates the number of BiCGStab steps should be taken
+ /// before switching.
+ ///
+ int _numberOfBiCgStabSteps = 1;
+
+ ///
+ /// Indicates the number of GPBiCG steps should be taken
+ /// before switching.
+ ///
+ int _numberOfGpbiCgSteps = 4;
+
+ ///
+ /// Gets or sets the number of steps taken with the BiCgStab algorithm
+ /// before switching over to the GPBiCG algorithm.
+ ///
+ public int NumberOfBiCgStabSteps
+ {
+ get => _numberOfBiCgStabSteps;
+
+ set
+ {
+ if (value < 0)
+ {
+ throw new ArgumentOutOfRangeException(nameof(value));
+ }
+
+ _numberOfBiCgStabSteps = value;
+ }
+ }
+
+ ///
+ /// Gets or sets the number of steps taken with the GPBiCG algorithm
+ /// before switching over to the BiCgStab algorithm.
+ ///
+ public int NumberOfGpBiCgSteps
+ {
+ get => _numberOfGpbiCgSteps;
+
+ set
+ {
+ if (value < 0)
+ {
+ throw new ArgumentOutOfRangeException(nameof(value));
+ }
+
+ _numberOfGpbiCgSteps = value;
+ }
+ }
+
+ ///
+ /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax
+ ///
+ /// Instance of the A.
+ /// Residual values in .
+ /// Instance of the x.
+ /// Instance of the b.
+ static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+ {
+ // -Ax = residual
+ matrix.Multiply(x, residual);
+ residual.Multiply(-1, residual);
+
+ // residual + b
+ residual.Add(b, residual);
+ }
+
+ ///
+ /// Decide if to do steps with BiCgStab
+ ///
+ /// Number of iteration
+ /// true if yes, otherwise false
+ bool ShouldRunBiCgStabSteps(int iterationNumber)
+ {
+ // Run the first steps as BiCGStab
+ // The number of steps past a whole iteration set
+ var difference = iterationNumber % (_numberOfBiCgStabSteps + _numberOfGpbiCgSteps);
+
+ // Do steps with BiCGStab if:
+ // - The difference is zero or more (i.e. we have done zero or more complete cycles)
+ // - The difference is less than the number of BiCGStab steps that should be taken
+ return (difference >= 0) && (difference < _numberOfBiCgStabSteps);
+ }
+
+ ///
+ /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
+ /// solution vector and x is the unknown vector.
+ ///
+ /// The coefficient matrix, A.
+ /// The solution vector, b
+ /// The result vector, x
+ /// The iterator to use to control when to stop iterating.
+ /// The preconditioner to use for approximations.
+ public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator, IPreconditioner preconditioner)
+ {
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException("Matrix must be square.", nameof(matrix));
+ }
+
+ if (input.Count != matrix.RowCount || result.Count != input.Count)
+ {
+ throw Matrix.DimensionsDontMatch(matrix, input, result);
+ }
+
+ if (iterator == null)
+ {
+ iterator = new Iterator();
+ }
+
+ if (preconditioner == null)
+ {
+ preconditioner = new UnitPreconditioner();
+ }
+
+ preconditioner.Initialize(matrix);
+
+ // x_0 is initial guess
+ // Take x_0 = 0
+ var xtemp = new DenseVector(input.Count);
+
+ // r_0 = b - Ax_0
+ // This is basically a SAXPY so it could be made a lot faster
+ var residuals = new DenseVector(matrix.RowCount);
+ CalculateTrueResidual(matrix, residuals, xtemp, input);
+
+ // Define the temporary scalars
+ Complex beta = 0;
+
+ // Define the temporary vectors
+ // rDash_0 = r_0
+ var rdash = DenseVector.OfVector(residuals);
+
+ // t_-1 = 0
+ var t = new DenseVector(residuals.Count);
+ var t0 = new DenseVector(residuals.Count);
+
+ // w_-1 = 0
+ var w = new DenseVector(residuals.Count);
+
+ // Define the remaining temporary vectors
+ var c = new DenseVector(residuals.Count);
+ var p = new DenseVector(residuals.Count);
+ var s = new DenseVector(residuals.Count);
+ var u = new DenseVector(residuals.Count);
+ var y = new DenseVector(residuals.Count);
+ var z = new DenseVector(residuals.Count);
+
+ var temp = new DenseVector(residuals.Count);
+ var temp2 = new DenseVector(residuals.Count);
+ var temp3 = new DenseVector(residuals.Count);
+
+ // for (k = 0, 1, .... )
+ var iterationNumber = 0;
+ while (iterator.DetermineStatus(iterationNumber, xtemp, input, residuals) == IterationStatus.Continue)
+ {
+ // p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1))
+ p.Subtract(u, temp);
+
+ temp.Multiply(beta, temp2);
+ residuals.Add(temp2, p);
+
+ // Solve M b_k = p_k
+ preconditioner.Approximate(p, temp);
+
+ // s_k = A b_k
+ matrix.Multiply(temp, s);
+
+ // alpha_k = (r*_0 * r_k) / (r*_0 * s_k)
+ var alpha = rdash.ConjugateDotProduct(residuals)/rdash.ConjugateDotProduct(s);
+
+ // y_k = t_(k-1) - r_k - alpha_k * w_(k-1) + alpha_k s_k
+ s.Subtract(w, temp);
+ t.Subtract(residuals, y);
+
+ temp.Multiply(alpha, temp2);
+ y.Add(temp2, temp3);
+ temp3.CopyTo(y);
+
+ // Store the old value of t in t0
+ t.CopyTo(t0);
+
+ // t_k = r_k - alpha_k s_k
+ s.Multiply(-alpha, temp2);
+ residuals.Add(temp2, t);
+
+ // Solve M d_k = t_k
+ preconditioner.Approximate(t, temp);
+
+ // c_k = A d_k
+ matrix.Multiply(temp, c);
+ var cdot = c.ConjugateDotProduct(c);
+
+ // cDot can only be zero if c is a zero vector
+ // We'll set cDot to 1 if it is zero to prevent NaN's
+ // Note that the calculation should continue fine because
+ // c.DotProduct(t) will be zero and so will c.DotProduct(y)
+ if (cdot.Real.AlmostEqualNumbersBetween(0, 1) && cdot.Imaginary.AlmostEqualNumbersBetween(0, 1))
+ {
+ cdot = 1.0;
+ }
+
+ // Even if we don't want to do any BiCGStab steps we'll still have
+ // to do at least one at the start to initialize the
+ // system, but we'll only have to take special measures
+ // if we don't do any so ...
+ var ctdot = c.ConjugateDotProduct(t);
+ Complex eta;
+ Complex sigma;
+ if (((_numberOfBiCgStabSteps == 0) && (iterationNumber == 0)) || ShouldRunBiCgStabSteps(iterationNumber))
+ {
+ // sigma_k = (c_k * t_k) / (c_k * c_k)
+ sigma = ctdot/cdot;
+
+ // eta_k = 0
+ eta = 0;
+ }
+ else
+ {
+ var ydot = y.ConjugateDotProduct(y);
+
+ // yDot can only be zero if y is a zero vector
+ // We'll set yDot to 1 if it is zero to prevent NaN's
+ // Note that the calculation should continue fine because
+ // y.DotProduct(t) will be zero and so will c.DotProduct(y)
+ if (ydot.Real.AlmostEqualNumbersBetween(0, 1) && ydot.Imaginary.AlmostEqualNumbersBetween(0, 1))
+ {
+ ydot = 1.0;
+ }
+
+ var ytdot = y.ConjugateDotProduct(t);
+ var cydot = c.ConjugateDotProduct(y);
+
+ var denom = (cdot*ydot) - (cydot*cydot);
+
+ // sigma_k = ((y_k * y_k)(c_k * t_k) - (y_k * t_k)(c_k * y_k)) / ((c_k * c_k)(y_k * y_k) - (y_k * c_k)(c_k * y_k))
+ sigma = ((ydot*ctdot) - (ytdot*cydot))/denom;
+
+ // eta_k = ((c_k * c_k)(y_k * t_k) - (y_k * c_k)(c_k * t_k)) / ((c_k * c_k)(y_k * y_k) - (y_k * c_k)(c_k * y_k))
+ eta = ((cdot*ytdot) - (cydot*ctdot))/denom;
+ }
+
+ // u_k = sigma_k s_k + eta_k (t_(k-1) - r_k + beta_(k-1) u_(k-1))
+ u.Multiply(beta, temp2);
+ t0.Add(temp2, temp);
+
+ temp.Subtract(residuals, temp3);
+ temp3.CopyTo(temp);
+ temp.Multiply(eta, temp);
+
+ s.Multiply(sigma, temp2);
+ temp.Add(temp2, u);
+
+ // z_k = sigma_k r_k +_ eta_k z_(k-1) - alpha_k u_k
+ z.Multiply(eta, z);
+ u.Multiply(-alpha, temp2);
+ z.Add(temp2, temp3);
+ temp3.CopyTo(z);
+
+ residuals.Multiply(sigma, temp2);
+ z.Add(temp2, temp3);
+ temp3.CopyTo(z);
+
+ // x_(k+1) = x_k + alpha_k p_k + z_k
+ p.Multiply(alpha, temp2);
+ xtemp.Add(temp2, temp3);
+ temp3.CopyTo(xtemp);
+
+ xtemp.Add(z, temp3);
+ temp3.CopyTo(xtemp);
+
+ // r_(k+1) = t_k - eta_k y_k - sigma_k c_k
+ // Copy the old residuals to a temp vector because we'll
+ // need those in the next step
+ residuals.CopyTo(t0);
+
+ y.Multiply(-eta, temp2);
+ t.Add(temp2, residuals);
+
+ c.Multiply(-sigma, temp2);
+ residuals.Add(temp2, temp3);
+ temp3.CopyTo(residuals);
+
+ // beta_k = alpha_k / sigma_k * (r*_0 * r_(k+1)) / (r*_0 * r_k)
+ // But first we check if there is a possible NaN. If so just reset beta to zero.
+ beta = (!sigma.Real.AlmostEqualNumbersBetween(0, 1) || !sigma.Imaginary.AlmostEqualNumbersBetween(0, 1)) ? alpha/sigma*rdash.ConjugateDotProduct(residuals)/rdash.ConjugateDotProduct(t0) : 0;
+
+ // w_k = c_k + beta_k s_k
+ s.Multiply(beta, temp2);
+ c.Add(temp2, w);
+
+ // Get the real value
+ preconditioner.Approximate(xtemp, result);
+
+ // Now check for convergence
+ if (iterator.DetermineStatus(iterationNumber, result, input, residuals) != IterationStatus.Continue)
+ {
+ // Recalculate the residuals and go round again. This is done to ensure that
+ // we have the proper residuals.
+ CalculateTrueResidual(matrix, residuals, result, input);
+ }
+
+ // Next iteration.
+ iterationNumber++;
+ }
+ }
+ }
+}
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs.meta
new file mode 100644
index 00000000..50ecf0b3
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: c0594509afa168f4fb460a5109e21204
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILU0Preconditioner.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILU0Preconditioner.cs
new file mode 100644
index 00000000..77fc9f90
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILU0Preconditioner.cs
@@ -0,0 +1,222 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2010 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;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// An incomplete, level 0, LU factorization preconditioner.
+ ///
+ ///
+ /// The ILU(0) algorithm was taken from:
+ /// Iterative methods for sparse linear systems
+ /// Yousef Saad
+ /// Algorithm is described in Chapter 10, section 10.3.2, page 275
+ ///
+ public sealed class ILU0Preconditioner : IPreconditioner
+ {
+ ///
+ /// The matrix holding the lower (L) and upper (U) matrices. The
+ /// decomposition matrices are combined to reduce storage.
+ ///
+ SparseMatrix _decompositionLU;
+
+ ///
+ /// Returns the upper triagonal matrix that was created during the LU decomposition.
+ ///
+ /// A new matrix containing the upper triagonal elements.
+ internal Matrix UpperTriangle()
+ {
+ var result = new SparseMatrix(_decompositionLU.RowCount);
+ for (var i = 0; i < _decompositionLU.RowCount; i++)
+ {
+ for (var j = i; j < _decompositionLU.ColumnCount; j++)
+ {
+ result[i, j] = _decompositionLU[i, j];
+ }
+ }
+
+ return result;
+ }
+
+ ///
+ /// Returns the lower triagonal matrix that was created during the LU decomposition.
+ ///
+ /// A new matrix containing the lower triagonal elements.
+ internal Matrix LowerTriangle()
+ {
+ var result = new SparseMatrix(_decompositionLU.RowCount);
+ for (var i = 0; i < _decompositionLU.RowCount; i++)
+ {
+ for (var j = 0; j <= i; j++)
+ {
+ if (i == j)
+ {
+ result[i, j] = 1.0;
+ }
+ else
+ {
+ result[i, j] = _decompositionLU[i, j];
+ }
+ }
+ }
+
+ return result;
+ }
+
+ ///
+ /// Initializes the preconditioner and loads the internal data structures.
+ ///
+ /// The matrix upon which the preconditioner is based.
+ /// If is .
+ /// If is not a square matrix.
+ public void Initialize(Matrix matrix)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException(nameof(matrix));
+ }
+
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException("Matrix must be square.", nameof(matrix));
+ }
+
+ _decompositionLU = SparseMatrix.OfMatrix(matrix);
+
+ // M == A
+ // for i = 2, ... , n do
+ // for k = 1, .... , i - 1 do
+ // if (i,k) == NZ(Z) then
+ // compute z(i,k) = z(i,k) / z(k,k);
+ // for j = k + 1, ...., n do
+ // if (i,j) == NZ(Z) then
+ // compute z(i,j) = z(i,j) - z(i,k) * z(k,j)
+ // end
+ // end
+ // end
+ // end
+ // end
+ for (var i = 0; i < _decompositionLU.RowCount; i++)
+ {
+ for (var k = 0; k < i; k++)
+ {
+ if (_decompositionLU[i, k] != 0.0)
+ {
+ var t = _decompositionLU[i, k]/_decompositionLU[k, k];
+ _decompositionLU[i, k] = t;
+ if (_decompositionLU[k, i] != 0.0)
+ {
+ _decompositionLU[i, i] = _decompositionLU[i, i] - (t*_decompositionLU[k, i]);
+ }
+
+ for (var j = k + 1; j < _decompositionLU.RowCount; j++)
+ {
+ if (j == i)
+ {
+ continue;
+ }
+
+ if (_decompositionLU[i, j] != 0.0)
+ {
+ _decompositionLU[i, j] = _decompositionLU[i, j] - (t*_decompositionLU[k, j]);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ ///
+ /// Approximates the solution to the matrix equation Ax = b.
+ ///
+ /// The right hand side vector.
+ /// The left hand side vector. Also known as the result vector.
+ public void Approximate(Vector rhs, Vector lhs)
+ {
+ if (_decompositionLU == null)
+ {
+ throw new ArgumentException("The requested matrix does not exist.");
+ }
+
+ if ((lhs.Count != rhs.Count) || (lhs.Count != _decompositionLU.RowCount))
+ {
+ throw new ArgumentException("All vectors must have the same dimensionality.");
+ }
+
+ // Solve:
+ // Lz = y
+ // Which gives
+ // for (int i = 1; i < matrix.RowLength; i++)
+ // {
+ // z_i = l_ii^-1 * (y_i - SUM_(j -1; i--)
+ // {
+ // x_i = u_ii^-1 * (z_i - SUM_(j > i) u_ij * x_j)
+ // }
+ for (var i = _decompositionLU.RowCount - 1; i > -1; i--)
+ {
+ _decompositionLU.Row(i, rowValues);
+
+ var sum = Complex.Zero;
+ for (var j = _decompositionLU.RowCount - 1; j > i; j--)
+ {
+ sum += rowValues[j]*lhs[j];
+ }
+
+ lhs[i] = 1/rowValues[i]*(lhs[i] - sum);
+ }
+ }
+ }
+}
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILU0Preconditioner.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILU0Preconditioner.cs.meta
new file mode 100644
index 00000000..9c45f31c
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILU0Preconditioner.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: ead7dde852fa3d54eafdd56cab4dc2fb
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILUTPPreconditioner.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILUTPPreconditioner.cs
new file mode 100644
index 00000000..ebf1ba99
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILUTPPreconditioner.cs
@@ -0,0 +1,864 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2010 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 System.Collections.Generic;
+using MathNet.Numerics.LinearAlgebra.Solvers;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// This class performs an Incomplete LU factorization with drop tolerance
+ /// and partial pivoting. The drop tolerance indicates which additional entries
+ /// will be dropped from the factorized LU matrices.
+ ///
+ ///
+ /// The ILUTP-Mem algorithm was taken from:
+ /// ILUTP_Mem: a Space-Efficient Incomplete LU Preconditioner
+ ///
+ /// Tzu-Yi Chen, Department of Mathematics and Computer Science,
+ /// Pomona College, Claremont CA 91711, USA
+ /// Published in:
+ /// Lecture Notes in Computer Science
+ /// Volume 3046 / 2004
+ /// pp. 20 - 28
+ /// Algorithm is described in Section 2, page 22
+ ///
+ public sealed class ILUTPPreconditioner : IPreconditioner
+ {
+ ///
+ /// The default fill level.
+ ///
+ public const double DefaultFillLevel = 200.0;
+
+ ///
+ /// The default drop tolerance.
+ ///
+ public const double DefaultDropTolerance = 0.0001;
+
+ ///
+ /// The decomposed upper triangular matrix.
+ ///
+ SparseMatrix _upper;
+
+ ///
+ /// The decomposed lower triangular matrix.
+ ///
+ SparseMatrix _lower;
+
+ ///
+ /// The array containing the pivot values.
+ ///
+ int[] _pivots;
+
+ ///
+ /// The fill level.
+ ///
+ double _fillLevel = DefaultFillLevel;
+
+ ///
+ /// The drop tolerance.
+ ///
+ double _dropTolerance = DefaultDropTolerance;
+
+ ///
+ /// The pivot tolerance.
+ ///
+ double _pivotTolerance;
+
+ ///
+ /// Initializes a new instance of the class with the default settings.
+ ///
+ public ILUTPPreconditioner()
+ {
+ }
+
+ ///
+ /// Initializes a new instance of the class with the specified settings.
+ ///
+ ///
+ /// The amount of fill that is allowed in the matrix. The value is a fraction of
+ /// the number of non-zero entries in the original matrix. Values should be positive.
+ ///
+ ///
+ /// The absolute drop tolerance which indicates below what absolute value an entry
+ /// will be dropped from the matrix. A drop tolerance of 0.0 means that no values
+ /// will be dropped. Values should always be positive.
+ ///
+ ///
+ /// The pivot tolerance which indicates at what level pivoting will take place. A
+ /// value of 0.0 means that no pivoting will take place.
+ ///
+ public ILUTPPreconditioner(double fillLevel, double dropTolerance, double pivotTolerance)
+ {
+ if (fillLevel < 0)
+ {
+ throw new ArgumentOutOfRangeException(nameof(fillLevel));
+ }
+
+ if (dropTolerance < 0)
+ {
+ throw new ArgumentOutOfRangeException(nameof(dropTolerance));
+ }
+
+ if (pivotTolerance < 0)
+ {
+ throw new ArgumentOutOfRangeException(nameof(pivotTolerance));
+ }
+
+ _fillLevel = fillLevel;
+ _dropTolerance = dropTolerance;
+ _pivotTolerance = pivotTolerance;
+ }
+
+ ///
+ /// Gets or sets the amount of fill that is allowed in the matrix. The
+ /// value is a fraction of the number of non-zero entries in the original
+ /// matrix. The standard value is 200.
+ ///
+ ///
+ ///
+ /// Values should always be positive and can be higher than 1.0. A value lower
+ /// than 1.0 means that the eventual preconditioner matrix will have fewer
+ /// non-zero entries as the original matrix. A value higher than 1.0 means that
+ /// the eventual preconditioner can have more non-zero values than the original
+ /// matrix.
+ ///
+ ///
+ /// Note that any changes to the FillLevel after creating the preconditioner
+ /// will invalidate the created preconditioner and will require a re-initialization of
+ /// the preconditioner.
+ ///
+ ///
+ /// Thrown if a negative value is provided.
+ public double FillLevel
+ {
+ get => _fillLevel;
+ set
+ {
+ if (value < 0)
+ {
+ throw new ArgumentOutOfRangeException(nameof(value));
+ }
+
+ _fillLevel = value;
+ }
+ }
+
+ ///
+ /// Gets or sets the absolute drop tolerance which indicates below what absolute value
+ /// an entry will be dropped from the matrix. The standard value is 0.0001.
+ ///
+ ///
+ ///
+ /// The values should always be positive and can be larger than 1.0. A low value will
+ /// keep more small numbers in the preconditioner matrix. A high value will remove
+ /// more small numbers from the preconditioner matrix.
+ ///
+ ///
+ /// Note that any changes to the DropTolerance after creating the preconditioner
+ /// will invalidate the created preconditioner and will require a re-initialization of
+ /// the preconditioner.
+ ///
+ ///
+ /// Thrown if a negative value is provided.
+ public double DropTolerance
+ {
+ get => _dropTolerance;
+ set
+ {
+ if (value < 0)
+ {
+ throw new ArgumentOutOfRangeException(nameof(value));
+ }
+
+ _dropTolerance = value;
+ }
+ }
+
+ ///
+ /// Gets or sets the pivot tolerance which indicates at what level pivoting will
+ /// take place. The standard value is 0.0 which means pivoting will never take place.
+ ///
+ ///
+ ///
+ /// The pivot tolerance is used to calculate if pivoting is necessary. Pivoting
+ /// will take place if any of the values in a row is bigger than the
+ /// diagonal value of that row divided by the pivot tolerance, i.e. pivoting
+ /// will take place if row(i,j) > row(i,i) / PivotTolerance for
+ /// any j that is not equal to i.
+ ///
+ ///
+ /// Note that any changes to the PivotTolerance after creating the preconditioner
+ /// will invalidate the created preconditioner and will require a re-initialization of
+ /// the preconditioner.
+ ///
+ ///
+ /// Thrown if a negative value is provided.
+ public double PivotTolerance
+ {
+ get => _pivotTolerance;
+ set
+ {
+ if (value < 0)
+ {
+ throw new ArgumentOutOfRangeException(nameof(value));
+ }
+
+ _pivotTolerance = value;
+ }
+ }
+
+ ///
+ /// Returns the upper triagonal matrix that was created during the LU decomposition.
+ ///
+ ///
+ /// This method is used for debugging purposes only and should normally not be used.
+ ///
+ /// A new matrix containing the upper triagonal elements.
+ internal Matrix UpperTriangle()
+ {
+ return _upper.Clone();
+ }
+
+ ///
+ /// Returns the lower triagonal matrix that was created during the LU decomposition.
+ ///
+ ///
+ /// This method is used for debugging purposes only and should normally not be used.
+ ///
+ /// A new matrix containing the lower triagonal elements.
+ internal Matrix LowerTriangle()
+ {
+ return _lower.Clone();
+ }
+
+ ///
+ /// Returns the pivot array. This array is not needed for normal use because
+ /// the preconditioner will return the solution vector values in the proper order.
+ ///
+ ///
+ /// This method is used for debugging purposes only and should normally not be used.
+ ///
+ /// The pivot array.
+ internal int[] Pivots()
+ {
+ var result = new int[_pivots.Length];
+ for (var i = 0; i < _pivots.Length; i++)
+ {
+ result[i] = _pivots[i];
+ }
+
+ return result;
+ }
+
+ ///
+ /// Initializes the preconditioner and loads the internal data structures.
+ ///
+ ///
+ /// The upon which this preconditioner is based. Note that the
+ /// method takes a general matrix type. However internally the data is stored
+ /// as a sparse matrix. Therefore it is not recommended to pass a dense matrix.
+ ///
+ /// If is .
+ /// If is not a square matrix.
+ public void Initialize(Matrix matrix)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException(nameof(matrix));
+ }
+
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException("Matrix must be square.", nameof(matrix));
+ }
+
+ var sparseMatrix = matrix as SparseMatrix ?? SparseMatrix.OfMatrix(matrix);
+
+ // The creation of the preconditioner follows the following algorithm.
+ // spaceLeft = lfilNnz * nnz(A)
+ // for i = 1, .. , n
+ // {
+ // w = a(i,*)
+ // for j = 1, .. , i - 1
+ // {
+ // if (w(j) != 0)
+ // {
+ // w(j) = w(j) / a(j,j)
+ // if (w(j) < dropTol)
+ // {
+ // w(j) = 0;
+ // }
+ // if (w(j) != 0)
+ // {
+ // w = w - w(j) * U(j,*)
+ // }
+ // }
+ // }
+ //
+ // for j = i, .. ,n
+ // {
+ // if w(j) <= dropTol * ||A(i,*)||
+ // {
+ // w(j) = 0
+ // }
+ // }
+ //
+ // spaceRow = spaceLeft / (n - i + 1) // Determine the space for this row
+ // lfil = spaceRow / 2 // space for this row of L
+ // l(i,j) = w(j) for j = 1, .. , i -1 // only the largest lfil elements
+ //
+ // lfil = spaceRow - nnz(L(i,:)) // space for this row of U
+ // u(i,j) = w(j) for j = i, .. , n // only the largest lfil - 1 elements
+ // w = 0
+ //
+ // if max(U(i,i + 1: n)) > U(i,i) / pivTol then // pivot if necessary
+ // {
+ // pivot by swapping the max and the diagonal entries
+ // Update L, U
+ // Update P
+ // }
+ // spaceLeft = spaceLeft - nnz(L(i,:)) - nnz(U(i,:))
+ // }
+ // Create the lower triangular matrix
+ _lower = new SparseMatrix(sparseMatrix.RowCount);
+
+ // Create the upper triangular matrix and copy the values
+ _upper = new SparseMatrix(sparseMatrix.RowCount);
+
+ // Create the pivot array
+ _pivots = new int[sparseMatrix.RowCount];
+ for (var i = 0; i < _pivots.Length; i++)
+ {
+ _pivots[i] = i;
+ }
+
+ var workVector = new DenseVector(sparseMatrix.RowCount);
+ var rowVector = new DenseVector(sparseMatrix.ColumnCount);
+ var indexSorting = new int[sparseMatrix.RowCount];
+
+ // spaceLeft = lfilNnz * nnz(A)
+ var spaceLeft = (int) _fillLevel*sparseMatrix.NonZerosCount;
+
+ // for i = 1, .. , n
+ for (var i = 0; i < sparseMatrix.RowCount; i++)
+ {
+ // w = a(i,*)
+ sparseMatrix.Row(i, workVector);
+
+ // pivot the row
+ PivotRow(workVector);
+ var vectorNorm = workVector.InfinityNorm();
+
+ // for j = 1, .. , i - 1)
+ for (var j = 0; j < i; j++)
+ {
+ // if (w(j) != 0)
+ // {
+ // w(j) = w(j) / a(j,j)
+ // if (w(j) < dropTol)
+ // {
+ // w(j) = 0;
+ // }
+ // if (w(j) != 0)
+ // {
+ // w = w - w(j) * U(j,*)
+ // }
+ if (workVector[j] != 0.0)
+ {
+ // Calculate the multiplication factors that go into the L matrix
+ workVector[j] = workVector[j]/_upper[j, j];
+ if (workVector[j].Magnitude < _dropTolerance)
+ {
+ workVector[j] = 0.0;
+ }
+
+ // Calculate the addition factor
+ if (workVector[j] != 0.0)
+ {
+ // vector update all in one go
+ _upper.Row(j, rowVector);
+
+ // zero out columnVector[k] because we don't need that
+ // one anymore for k = 0 to k = j
+ for (var k = 0; k <= j; k++)
+ {
+ rowVector[k] = 0.0;
+ }
+
+ rowVector.Multiply(workVector[j], rowVector);
+ workVector.Subtract(rowVector, workVector);
+ }
+ }
+ }
+
+ // for j = i, .. ,n
+ for (var j = i; j < sparseMatrix.RowCount; j++)
+ {
+ // if w(j) <= dropTol * ||A(i,*)||
+ // {
+ // w(j) = 0
+ // }
+ if (workVector[j].Magnitude <= _dropTolerance*vectorNorm)
+ {
+ workVector[j] = 0.0;
+ }
+ }
+
+ // spaceRow = spaceLeft / (n - i + 1) // Determine the space for this row
+ var spaceRow = spaceLeft/(sparseMatrix.RowCount - i + 1);
+
+ // lfil = spaceRow / 2 // space for this row of L
+ var fillLevel = spaceRow/2;
+ FindLargestItems(0, i - 1, indexSorting, workVector);
+
+ // l(i,j) = w(j) for j = 1, .. , i -1 // only the largest lfil elements
+ var lowerNonZeroCount = 0;
+ var count = 0;
+ for (var j = 0; j < i; j++)
+ {
+ if ((count > fillLevel) || (indexSorting[j] == -1))
+ {
+ break;
+ }
+
+ _lower[i, indexSorting[j]] = workVector[indexSorting[j]];
+ count += 1;
+ lowerNonZeroCount += 1;
+ }
+
+ FindLargestItems(i + 1, sparseMatrix.RowCount - 1, indexSorting, workVector);
+
+ // lfil = spaceRow - nnz(L(i,:)) // space for this row of U
+ fillLevel = spaceRow - lowerNonZeroCount;
+
+ // u(i,j) = w(j) for j = i + 1, .. , n // only the largest lfil - 1 elements
+ var upperNonZeroCount = 0;
+ count = 0;
+ for (var j = 0; j < sparseMatrix.RowCount - i; j++)
+ {
+ if ((count > fillLevel - 1) || (indexSorting[j] == -1))
+ {
+ break;
+ }
+
+ _upper[i, indexSorting[j]] = workVector[indexSorting[j]];
+ count += 1;
+ upperNonZeroCount += 1;
+ }
+
+ // Simply copy the diagonal element. Next step is to see if we pivot
+ _upper[i, i] = workVector[i];
+
+ // if max(U(i,i + 1: n)) > U(i,i) / pivTol then // pivot if necessary
+ // {
+ // pivot by swapping the max and the diagonal entries
+ // Update L, U
+ // Update P
+ // }
+
+ // Check if we really need to pivot. If (i+1) >=(mCoefficientMatrix.Rows -1) then
+ // we are working on the last row. That means that there is only one number
+ // And pivoting is useless. Also the indexSorting array will only contain
+ // -1 values.
+ if ((i + 1) < (sparseMatrix.RowCount - 1))
+ {
+ if (workVector[i].Magnitude < _pivotTolerance*workVector[indexSorting[0]].Magnitude)
+ {
+ // swap columns of u (which holds the values of A in the
+ // sections that haven't been partitioned yet.
+ SwapColumns(_upper, i, indexSorting[0]);
+
+ // Update P
+ var temp = _pivots[i];
+ _pivots[i] = _pivots[indexSorting[0]];
+ _pivots[indexSorting[0]] = temp;
+ }
+ }
+
+ // spaceLeft = spaceLeft - nnz(L(i,:)) - nnz(U(i,:))
+ spaceLeft -= lowerNonZeroCount + upperNonZeroCount;
+ }
+
+ for (var i = 0; i < _lower.RowCount; i++)
+ {
+ _lower[i, i] = 1.0;
+ }
+ }
+
+ ///
+ /// Pivot elements in the according to internal pivot array
+ ///
+ /// Row to pivot in
+ void PivotRow(Vector row)
+ {
+ var knownPivots = new Dictionary();
+
+ // pivot the row
+ for (var i = 0; i < row.Count; i++)
+ {
+ if ((_pivots[i] != i) && (!PivotMapFound(knownPivots, i)))
+ {
+ // store the pivots in the hashtable
+ knownPivots.Add(_pivots[i], i);
+
+ var t = row[i];
+ row[i] = row[_pivots[i]];
+ row[_pivots[i]] = t;
+ }
+ }
+ }
+
+ ///
+ /// Was pivoting already performed
+ ///
+ /// Pivots already done
+ /// Current item to pivot
+ /// true if performed, otherwise false
+ bool PivotMapFound(Dictionary knownPivots, int currentItem)
+ {
+ if (knownPivots.TryGetValue(_pivots[currentItem], out var knownPivot) && knownPivot.Equals(currentItem))
+ {
+ return true;
+ }
+
+ if (knownPivots.TryGetValue(currentItem, out var pivot) && pivot.Equals(_pivots[currentItem]))
+ {
+ return true;
+ }
+
+ return false;
+ }
+
+ ///
+ /// Swap columns in the
+ ///
+ /// Source .
+ /// First column index to swap
+ /// Second column index to swap
+ static void SwapColumns(Matrix matrix, int firstColumn, int secondColumn)
+ {
+ for (var i = 0; i < matrix.RowCount; i++)
+ {
+ var temp = matrix[i, firstColumn];
+ matrix[i, firstColumn] = matrix[i, secondColumn];
+ matrix[i, secondColumn] = temp;
+ }
+ }
+
+ ///
+ /// Sort vector descending, not changing vector but placing sorted indices to
+ ///
+ /// Start sort form
+ /// Sort till upper bound
+ /// Array with sorted vector indices
+ /// Source
+ static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+ {
+ // Copy the indices for the values into the array
+ for (var i = 0; i < upperBound + 1 - lowerBound; i++)
+ {
+ sortedIndices[i] = lowerBound + i;
+ }
+
+ for (var i = upperBound + 1 - lowerBound; i < sortedIndices.Length; i++)
+ {
+ sortedIndices[i] = -1;
+ }
+
+ // Sort the first set of items.
+ // Sorting starts at index 0 because the index array
+ // starts at zero
+ // and ends at index upperBound - lowerBound
+ ILUTPElementSorter.SortDoubleIndicesDecreasing(0, upperBound - lowerBound, sortedIndices, values);
+ }
+
+ ///
+ /// Approximates the solution to the matrix equation Ax = b.
+ ///
+ /// The right hand side vector.
+ /// The left hand side vector. Also known as the result vector.
+ public void Approximate(Vector rhs, Vector lhs)
+ {
+ if (_upper == null)
+ {
+ throw new ArgumentException("The requested matrix does not exist.");
+ }
+
+ if ((lhs.Count != rhs.Count) || (lhs.Count != _upper.RowCount))
+ {
+ throw new ArgumentException("All vectors must have the same dimensionality.", nameof(rhs));
+ }
+
+ // Solve equation here
+ // Pivot(vector, result);
+ // Solve L*Y = B(piv,:)
+ var rowValues = new DenseVector(_lower.RowCount);
+ for (var i = 0; i < _lower.RowCount; i++)
+ {
+ _lower.Row(i, rowValues);
+
+ var sum = Complex.Zero;
+ for (var j = 0; j < i; j++)
+ {
+ sum += rowValues[j]*lhs[j];
+ }
+
+ lhs[i] = rhs[i] - sum;
+ }
+
+ // Solve U*X = Y;
+ for (var i = _upper.RowCount - 1; i > -1; i--)
+ {
+ _upper.Row(i, rowValues);
+
+ var sum = Complex.Zero;
+ for (var j = _upper.RowCount - 1; j > i; j--)
+ {
+ sum += rowValues[j]*lhs[j];
+ }
+
+ lhs[i] = 1/rowValues[i]*(lhs[i] - sum);
+ }
+
+ // We have a column pivot so we only need to pivot the
+ // end result not the incoming right hand side vector
+ var temp = lhs.Clone();
+
+ Pivot(temp, lhs);
+ }
+
+ ///
+ /// Pivot elements in according to internal pivot array
+ ///
+ /// Source .
+ /// Result after pivoting.
+ void Pivot(Vector vector, Vector result)
+ {
+ for (var i = 0; i < _pivots.Length; i++)
+ {
+ result[i] = vector[_pivots[i]];
+ }
+ }
+ }
+
+ ///
+ /// An element sort algorithm for the class.
+ ///
+ ///
+ /// This sort algorithm is used to sort the columns in a sparse matrix based on
+ /// the value of the element on the diagonal of the matrix.
+ ///
+ internal static class ILUTPElementSorter
+ {
+ ///
+ /// Sorts the elements of the vector in decreasing
+ /// fashion. The vector itself is not affected.
+ ///
+ /// The starting index.
+ /// The stopping index.
+ /// An array that will contain the sorted indices once the algorithm finishes.
+ /// The that contains the values that need to be sorted.
+ public static void SortDoubleIndicesDecreasing(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+ {
+ // Move all the indices that we're interested in to the beginning of the
+ // array. Ignore the rest of the indices.
+ if (lowerBound > 0)
+ {
+ for (var i = 0; i < (upperBound - lowerBound + 1); i++)
+ {
+ Exchange(sortedIndices, i, i + lowerBound);
+ }
+
+ upperBound -= lowerBound;
+ lowerBound = 0;
+ }
+
+ HeapSortDoublesIndices(lowerBound, upperBound, sortedIndices, values);
+ }
+
+ ///
+ /// Sorts the elements of the vector in decreasing
+ /// fashion using heap sort algorithm. The vector itself is not affected.
+ ///
+ /// The starting index.
+ /// The stopping index.
+ /// An array that will contain the sorted indices once the algorithm finishes.
+ /// The that contains the values that need to be sorted.
+ private static void HeapSortDoublesIndices(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+ {
+ var start = ((upperBound - lowerBound + 1) / 2) - 1 + lowerBound;
+ var end = (upperBound - lowerBound + 1) - 1 + lowerBound;
+
+ BuildDoubleIndexHeap(start, upperBound - lowerBound + 1, sortedIndices, values);
+
+ while (end >= lowerBound)
+ {
+ Exchange(sortedIndices, end, lowerBound);
+ SiftDoubleIndices(sortedIndices, values, lowerBound, end);
+ end -= 1;
+ }
+ }
+
+ ///
+ /// Build heap for double indices
+ ///
+ /// Root position
+ /// Length of
+ /// Indices of
+ /// Target
+ private static void BuildDoubleIndexHeap(int start, int count, int[] sortedIndices, Vector values)
+ {
+ while (start >= 0)
+ {
+ SiftDoubleIndices(sortedIndices, values, start, count);
+ start -= 1;
+ }
+ }
+
+ ///
+ /// Sift double indices
+ ///
+ /// Indices of
+ /// Target
+ /// Root position
+ /// Length of
+ private static void SiftDoubleIndices(int[] sortedIndices, Vector values, int begin, int count)
+ {
+ var root = begin;
+
+ while (root * 2 < count)
+ {
+ var child = root * 2;
+ if ((child < count - 1) && (values[sortedIndices[child]].Magnitude > values[sortedIndices[child + 1]].Magnitude))
+ {
+ child += 1;
+ }
+
+ if (values[sortedIndices[root]].Magnitude <= values[sortedIndices[child]].Magnitude)
+ {
+ return;
+ }
+
+ Exchange(sortedIndices, root, child);
+ root = child;
+ }
+ }
+
+ ///
+ /// Sorts the given integers in a decreasing fashion.
+ ///
+ /// The values.
+ public static void SortIntegersDecreasing(int[] values)
+ {
+ HeapSortIntegers(values, values.Length);
+ }
+
+ ///
+ /// Sort the given integers in a decreasing fashion using heapsort algorithm
+ ///
+ /// Array of values to sort
+ /// Length of
+ private static void HeapSortIntegers(int[] values, int count)
+ {
+ var start = (count / 2) - 1;
+ var end = count - 1;
+
+ BuildHeap(values, start, count);
+
+ while (end >= 0)
+ {
+ Exchange(values, end, 0);
+ Sift(values, 0, end);
+ end -= 1;
+ }
+ }
+
+ ///
+ /// Build heap
+ ///
+ /// Target values array
+ /// Root position
+ /// Length of
+ private static void BuildHeap(int[] values, int start, int count)
+ {
+ while (start >= 0)
+ {
+ Sift(values, start, count);
+ start -= 1;
+ }
+ }
+
+ ///
+ /// Sift values
+ ///
+ /// Target value array
+ /// Root position
+ /// Length of
+ private static void Sift(int[] values, int start, int count)
+ {
+ var root = start;
+
+ while (root * 2 < count)
+ {
+ var child = root * 2;
+ if ((child < count - 1) && (values[child] > values[child + 1]))
+ {
+ child += 1;
+ }
+
+ if (values[root] > values[child])
+ {
+ Exchange(values, root, child);
+ root = child;
+ }
+ else
+ {
+ return;
+ }
+ }
+ }
+
+ ///
+ /// Exchange values in array
+ ///
+ /// Target values array
+ /// First value to exchange
+ /// Second value to exchange
+ private static void Exchange(int[] values, int first, int second)
+ {
+ var t = values[first];
+ values[first] = values[second];
+ values[second] = t;
+ }
+ }
+}
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILUTPPreconditioner.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILUTPPreconditioner.cs.meta
new file mode 100644
index 00000000..cd1227d1
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/ILUTPPreconditioner.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: 79bdca8be0ff5294eae74e938739a419
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MILU0Preconditioner.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MILU0Preconditioner.cs
new file mode 100644
index 00000000..aab88c5b
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MILU0Preconditioner.cs
@@ -0,0 +1,261 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// 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;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// A simple milu(0) preconditioner.
+ ///
+ ///
+ /// Original Fortran code by Yousef 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", nameof(matrix));
+ }
+
+ // Dimension of matrix
+ int n = csr.RowCount;
+ if (n != csr.ColumnCount)
+ {
+ throw new ArgumentException("Matrix must be square.", nameof(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 NumericalBreakdownException("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("The requested matrix does not exist.");
+ }
+
+ if ((result.Count != input.Count) || (result.Count != _diag.Length))
+ {
+ throw new ArgumentException("All vectors must have the same dimensionality.");
+ }
+
+ 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).
+ /// True if the modified/MILU algorithm should be used (recommended)
+ /// 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/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MILU0Preconditioner.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MILU0Preconditioner.cs.meta
new file mode 100644
index 00000000..4ec22164
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MILU0Preconditioner.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: 7a51132033f393f498c4fc6c6d141f2f
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs
new file mode 100644
index 00000000..04232d17
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs
@@ -0,0 +1,538 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// 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 System.Collections.Generic;
+using System.Diagnostics;
+using System.Linq;
+using MathNet.Numerics.Distributions;
+using MathNet.Numerics.LinearAlgebra.Solvers;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// A Multiple-Lanczos Bi-Conjugate Gradient stabilized iterative matrix solver.
+ ///
+ ///
+ ///
+ /// The Multiple-Lanczos Bi-Conjugate Gradient stabilized (ML(k)-BiCGStab) solver is an 'improvement'
+ /// of the standard BiCgStab solver.
+ ///
+ ///
+ /// The algorithm was taken from:
+ /// ML(k)BiCGSTAB: A BiCGSTAB variant based on multiple Lanczos starting vectors
+ ///
+ /// Man-Chung Yeung and Tony F. Chan
+ ///
+ /// SIAM Journal of Scientific Computing
+ ///
+ /// Volume 21, Number 4, pp. 1263 - 1290
+ ///
+ ///
+ /// The example code below provides an indication of the possible use of the
+ /// solver.
+ ///
+ ///
+ public sealed class MlkBiCgStab : IIterativeSolver
+ {
+ ///
+ /// The default number of starting vectors.
+ ///
+ const int DefaultNumberOfStartingVectors = 50;
+
+ ///
+ /// The collection of starting vectors which are used as the basis for the Krylov sub-space.
+ ///
+ IList> _startingVectors;
+
+ ///
+ /// The number of starting vectors used by the algorithm
+ ///
+ int _numberOfStartingVectors = DefaultNumberOfStartingVectors;
+
+ ///
+ /// Gets or sets the number of starting vectors.
+ ///
+ ///
+ /// Must be larger than 1 and smaller than the number of variables in the matrix that
+ /// for which this solver will be used.
+ ///
+ public int NumberOfStartingVectors
+ {
+ [DebuggerStepThrough]
+ get => _numberOfStartingVectors;
+
+ [DebuggerStepThrough]
+ set
+ {
+ if (value <= 1)
+ {
+ throw new ArgumentOutOfRangeException(nameof(value));
+ }
+
+ _numberOfStartingVectors = value;
+ }
+ }
+
+ ///
+ /// Resets the number of starting vectors to the default value.
+ ///
+ public void ResetNumberOfStartingVectors()
+ {
+ _numberOfStartingVectors = DefaultNumberOfStartingVectors;
+ }
+
+ ///
+ /// Gets or sets a series of orthonormal vectors which will be used as basis for the
+ /// Krylov sub-space.
+ ///
+ public IList> StartingVectors
+ {
+ [DebuggerStepThrough]
+ get => _startingVectors;
+
+ [DebuggerStepThrough]
+ set
+ {
+ if ((value == null) || (value.Count == 0))
+ {
+ _startingVectors = null;
+ }
+ else
+ {
+ _startingVectors = value;
+ }
+ }
+ }
+
+ ///
+ /// Gets the number of starting vectors to create
+ ///
+ /// Maximum number
+ /// Number of variables
+ /// Number of starting vectors to create
+ static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables)
+ {
+ // Create no more starting vectors than the size of the problem - 1
+ return Math.Min(maximumNumberOfStartingVectors, (numberOfVariables - 1));
+ }
+
+ ///
+ /// Returns an array of starting vectors.
+ ///
+ /// The maximum number of starting vectors that should be created.
+ /// The number of variables.
+ ///
+ /// An array with starting vectors. The array will never be larger than the
+ /// but it may be smaller if
+ /// the is smaller than
+ /// the .
+ ///
+ static IList> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
+ {
+ // Create no more starting vectors than the size of the problem - 1
+ // Get random values and then orthogonalize them with
+ // modified Gramm - Schmidt
+ var count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables);
+
+ // Get a random set of samples based on the standard normal distribution with
+ // mean = 0 and sd = 1
+ var distribution = new Normal();
+
+ var matrix = new DenseMatrix(numberOfVariables, count);
+ for (var i = 0; i < matrix.ColumnCount; i++)
+ {
+ var samples = new Complex[matrix.RowCount];
+ var samplesRe = distribution.Samples().Take(matrix.RowCount).ToArray();
+ var samplesIm = distribution.Samples().Take(matrix.RowCount).ToArray();
+ for (int j = 0; j < matrix.RowCount; j++)
+ {
+ samples[j] = new Complex(samplesRe[j], samplesIm[j]);
+ }
+
+ // Set the column
+ matrix.SetColumn(i, samples);
+ }
+
+ // Compute the orthogonalization.
+ var gs = matrix.GramSchmidt();
+ var orthogonalMatrix = gs.Q;
+
+ // Now transfer this to vectors
+ var result = new List>(orthogonalMatrix.ColumnCount);
+ for (var i = 0; i < orthogonalMatrix.ColumnCount; i++)
+ {
+ result.Add(orthogonalMatrix.Column(i));
+
+ // Normalize the result vector
+ result[i].Multiply(1 / result[i].L2Norm(), result[i]);
+ }
+
+ return result;
+ }
+
+ ///
+ /// Create random vectors array
+ ///
+ /// Number of vectors
+ /// Size of each vector
+ /// Array of random vectors
+ static Vector[] CreateVectorArray(int arraySize, int vectorSize)
+ {
+ var result = new Vector[arraySize];
+ for (var i = 0; i < result.Length; i++)
+ {
+ result[i] = new DenseVector(vectorSize);
+ }
+
+ return result;
+ }
+
+ ///
+ /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax
+ ///
+ /// Source A.
+ /// Residual data.
+ /// x data.
+ /// b data.
+ static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+ {
+ // -Ax = residual
+ matrix.Multiply(x, residual);
+ residual.Multiply(-1, residual);
+
+ // residual + b
+ residual.Add(b, residual);
+ }
+
+ ///
+ /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
+ /// solution vector and x is the unknown vector.
+ ///
+ /// The coefficient matrix, A.
+ /// The solution vector, b
+ /// The result vector, x
+ /// The iterator to use to control when to stop iterating.
+ /// The preconditioner to use for approximations.
+ public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator, IPreconditioner preconditioner)
+ {
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException("Matrix must be square.", nameof(matrix));
+ }
+
+ if (input.Count != matrix.RowCount || result.Count != input.Count)
+ {
+ throw Matrix.DimensionsDontMatch(matrix, input, result);
+ }
+
+ if (iterator == null)
+ {
+ iterator = new Iterator();
+ }
+
+ if (preconditioner == null)
+ {
+ preconditioner = new UnitPreconditioner();
+ }
+
+ preconditioner.Initialize(matrix);
+
+ // Choose an initial guess x_0
+ // Take x_0 = 0
+ var xtemp = new DenseVector(input.Count);
+
+ // Choose k vectors q_1, q_2, ..., q_k
+ // Build a new set if:
+ // a) the stored set doesn't exist (i.e. == null)
+ // b) Is of an incorrect length (i.e. too long)
+ // c) The vectors are of an incorrect length (i.e. too long or too short)
+ var useOld = false;
+ if (_startingVectors != null)
+ {
+ // We don't accept collections with zero starting vectors so ...
+ if (_startingVectors.Count <= NumberOfStartingVectorsToCreate(_numberOfStartingVectors, input.Count))
+ {
+ // Only check the first vector for sizing. If that matches we assume the
+ // other vectors match too. If they don't the process will crash
+ if (_startingVectors[0].Count == input.Count)
+ {
+ useOld = true;
+ }
+ }
+ }
+
+ _startingVectors = useOld ? _startingVectors : CreateStartingVectors(_numberOfStartingVectors, input.Count);
+
+ // Store the number of starting vectors. Not really necessary but easier to type :)
+ var k = _startingVectors.Count;
+
+ // r_0 = b - Ax_0
+ // This is basically a SAXPY so it could be made a lot faster
+ var residuals = new DenseVector(matrix.RowCount);
+ CalculateTrueResidual(matrix, residuals, xtemp, input);
+
+ // Define the temporary values
+ var c = new Complex[k];
+
+ // Define the temporary vectors
+ var gtemp = new DenseVector(residuals.Count);
+
+ var u = new DenseVector(residuals.Count);
+ var utemp = new DenseVector(residuals.Count);
+ var temp = new DenseVector(residuals.Count);
+ var temp1 = new DenseVector(residuals.Count);
+ var temp2 = new DenseVector(residuals.Count);
+
+ var zd = new DenseVector(residuals.Count);
+ var zg = new DenseVector(residuals.Count);
+ var zw = new DenseVector(residuals.Count);
+
+ var d = CreateVectorArray(_startingVectors.Count, residuals.Count);
+
+ // g_0 = r_0
+ var g = CreateVectorArray(_startingVectors.Count, residuals.Count);
+ residuals.CopyTo(g[k - 1]);
+
+ var w = CreateVectorArray(_startingVectors.Count, residuals.Count);
+
+ // FOR (j = 0, 1, 2 ....)
+ var iterationNumber = 0;
+ while (iterator.DetermineStatus(iterationNumber, xtemp, input, residuals) == IterationStatus.Continue)
+ {
+ // SOLVE M g~_((j-1)k+k) = g_((j-1)k+k)
+ preconditioner.Approximate(g[k - 1], gtemp);
+
+ // w_((j-1)k+k) = A g~_((j-1)k+k)
+ matrix.Multiply(gtemp, w[k - 1]);
+
+ // c_((j-1)k+k) = q^T_1 w_((j-1)k+k)
+ c[k - 1] = _startingVectors[0].ConjugateDotProduct(w[k - 1]);
+ if (c[k - 1].Real.AlmostEqualNumbersBetween(0, 1) && c[k - 1].Imaginary.AlmostEqualNumbersBetween(0, 1))
+ {
+ throw new NumericalBreakdownException();
+ }
+
+ // alpha_(jk+1) = q^T_1 r_((j-1)k+k) / c_((j-1)k+k)
+ var alpha = _startingVectors[0].ConjugateDotProduct(residuals)/c[k - 1];
+
+ // u_(jk+1) = r_((j-1)k+k) - alpha_(jk+1) w_((j-1)k+k)
+ w[k - 1].Multiply(-alpha, temp);
+ residuals.Add(temp, u);
+
+ // SOLVE M u~_(jk+1) = u_(jk+1)
+ preconditioner.Approximate(u, temp1);
+ temp1.CopyTo(utemp);
+
+ // rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2
+ matrix.Multiply(temp1, temp);
+ var rho = temp.ConjugateDotProduct(temp);
+
+ // If rho is zero then temp is a zero vector and we're probably
+ // about to have zero residuals (i.e. an exact solution).
+ // So set rho to 1.0 because in the next step it will turn to zero.
+ if (rho.Real.AlmostEqualNumbersBetween(0, 1) && rho.Imaginary.AlmostEqualNumbersBetween(0, 1))
+ {
+ rho = 1.0;
+ }
+
+ rho = -u.ConjugateDotProduct(temp)/rho;
+
+ // r_(jk+1) = rho_(j+1) A u~_(jk+1) + u_(jk+1)
+ u.CopyTo(residuals);
+
+ // Reuse temp
+ temp.Multiply(rho, temp);
+ residuals.Add(temp, temp2);
+ temp2.CopyTo(residuals);
+
+ // x_(jk+1) = x_((j-1)k_k) - rho_(j+1) u~_(jk+1) + alpha_(jk+1) g~_((j-1)k+k)
+ utemp.Multiply(-rho, temp);
+ xtemp.Add(temp, temp2);
+ temp2.CopyTo(xtemp);
+
+ gtemp.Multiply(alpha, gtemp);
+ xtemp.Add(gtemp, temp2);
+ temp2.CopyTo(xtemp);
+
+ // Check convergence and stop if we are converged.
+ if (iterator.DetermineStatus(iterationNumber, xtemp, input, residuals) != IterationStatus.Continue)
+ {
+ // Calculate the true residual
+ CalculateTrueResidual(matrix, residuals, xtemp, input);
+
+ // Now recheck the convergence
+ if (iterator.DetermineStatus(iterationNumber, xtemp, input, residuals) != IterationStatus.Continue)
+ {
+ // We're all good now.
+ // Exit from the while loop.
+ break;
+ }
+ }
+
+ // FOR (i = 1,2, ...., k)
+ for (var i = 0; i < k; i++)
+ {
+ // z_d = u_(jk+1)
+ u.CopyTo(zd);
+
+ // z_g = r_(jk+i)
+ residuals.CopyTo(zg);
+
+ // z_w = 0
+ zw.Clear();
+
+ // FOR (s = i, ...., k-1) AND j >= 1
+ Complex beta;
+ if (iterationNumber >= 1)
+ {
+ for (var s = i; s < k - 1; s++)
+ {
+ // beta^(jk+i)_((j-1)k+s) = -q^t_(s+1) z_d / c_((j-1)k+s)
+ beta = -_startingVectors[s + 1].ConjugateDotProduct(zd)/c[s];
+
+ // z_d = z_d + beta^(jk+i)_((j-1)k+s) d_((j-1)k+s)
+ d[s].Multiply(beta, temp);
+ zd.Add(temp, temp2);
+ temp2.CopyTo(zd);
+
+ // z_g = z_g + beta^(jk+i)_((j-1)k+s) g_((j-1)k+s)
+ g[s].Multiply(beta, temp);
+ zg.Add(temp, temp2);
+ temp2.CopyTo(zg);
+
+ // z_w = z_w + beta^(jk+i)_((j-1)k+s) w_((j-1)k+s)
+ w[s].Multiply(beta, temp);
+ zw.Add(temp, temp2);
+ temp2.CopyTo(zw);
+ }
+ }
+
+ beta = rho*c[k - 1];
+ if (beta.Real.AlmostEqualNumbersBetween(0, 1) && beta.Imaginary.AlmostEqualNumbersBetween(0, 1))
+ {
+ throw new NumericalBreakdownException();
+ }
+
+ // beta^(jk+i)_((j-1)k+k) = -(q^T_1 (r_(jk+1) + rho_(j+1) z_w)) / (rho_(j+1) c_((j-1)k+k))
+ zw.Multiply(rho, temp2);
+ residuals.Add(temp2, temp);
+ beta = -_startingVectors[0].ConjugateDotProduct(temp)/beta;
+
+ // z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k)
+ g[k - 1].Multiply(beta, temp);
+ zg.Add(temp, temp2);
+ temp2.CopyTo(zg);
+
+ // z_w = rho_(j+1) (z_w + beta^(jk+i)_((j-1)k+k) w_((j-1)k+k))
+ w[k - 1].Multiply(beta, temp);
+ zw.Add(temp, temp2);
+ temp2.CopyTo(zw);
+ zw.Multiply(rho, zw);
+
+ // z_d = r_(jk+i) + z_w
+ residuals.Add(zw, zd);
+
+ // FOR (s = 1, ... i - 1)
+ for (var s = 0; s < i - 1; s++)
+ {
+ // beta^(jk+i)_(jk+s) = -q^T_s+1 z_d / c_(jk+s)
+ beta = -_startingVectors[s + 1].ConjugateDotProduct(zd)/c[s];
+
+ // z_d = z_d + beta^(jk+i)_(jk+s) * d_(jk+s)
+ d[s].Multiply(beta, temp);
+ zd.Add(temp, temp2);
+ temp2.CopyTo(zd);
+
+ // z_g = z_g + beta^(jk+i)_(jk+s) * g_(jk+s)
+ g[s].Multiply(beta, temp);
+ zg.Add(temp, temp2);
+ temp2.CopyTo(zg);
+ }
+
+ // d_(jk+i) = z_d - u_(jk+i)
+ zd.Subtract(u, d[i]);
+
+ // g_(jk+i) = z_g + z_w
+ zg.Add(zw, g[i]);
+
+ // IF (i < k - 1)
+ if (i < k - 1)
+ {
+ // c_(jk+1) = q^T_i+1 d_(jk+i)
+ c[i] = _startingVectors[i + 1].ConjugateDotProduct(d[i]);
+ if (c[i].Real.AlmostEqualNumbersBetween(0, 1) && c[i].Imaginary.AlmostEqualNumbersBetween(0, 1))
+ {
+ throw new NumericalBreakdownException();
+ }
+
+ // alpha_(jk+i+1) = q^T_(i+1) u_(jk+i) / c_(jk+i)
+ alpha = _startingVectors[i + 1].ConjugateDotProduct(u)/c[i];
+
+ // u_(jk+i+1) = u_(jk+i) - alpha_(jk+i+1) d_(jk+i)
+ d[i].Multiply(-alpha, temp);
+ u.Add(temp, temp2);
+ temp2.CopyTo(u);
+
+ // SOLVE M g~_(jk+i) = g_(jk+i)
+ preconditioner.Approximate(g[i], gtemp);
+
+ // x_(jk+i+1) = x_(jk+i) + rho_(j+1) alpha_(jk+i+1) g~_(jk+i)
+ gtemp.Multiply(rho*alpha, temp);
+ xtemp.Add(temp, temp2);
+ temp2.CopyTo(xtemp);
+
+ // w_(jk+i) = A g~_(jk+i)
+ matrix.Multiply(gtemp, w[i]);
+
+ // r_(jk+i+1) = r_(jk+i) - rho_(j+1) alpha_(jk+i+1) w_(jk+i)
+ w[i].Multiply(-rho*alpha, temp);
+ residuals.Add(temp, temp2);
+ temp2.CopyTo(residuals);
+
+ // We can check the residuals here if they're close
+ if (iterator.DetermineStatus(iterationNumber, xtemp, input, residuals) != IterationStatus.Continue)
+ {
+ // Recalculate the residuals and go round again. This is done to ensure that
+ // we have the proper residuals.
+ CalculateTrueResidual(matrix, residuals, xtemp, input);
+ }
+ }
+ } // END ITERATION OVER i
+
+ iterationNumber++;
+ }
+
+ // copy the temporary result to the real result vector
+ xtemp.CopyTo(result);
+ }
+ }
+}
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs.meta
new file mode 100644
index 00000000..edf605b1
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: 9d3c7bcc0c90dfd4b91702a199ae813e
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs
new file mode 100644
index 00000000..f88edfb4
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs
@@ -0,0 +1,273 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// 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;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// A Transpose Free Quasi-Minimal Residual (TFQMR) iterative matrix solver.
+ ///
+ ///
+ ///
+ /// The TFQMR algorithm was taken from:
+ /// Iterative methods for sparse linear systems.
+ ///
+ /// Yousef Saad
+ ///
+ /// Algorithm is described in Chapter 7, section 7.4.3, page 219
+ ///
+ ///
+ /// The example code below provides an indication of the possible use of the
+ /// solver.
+ ///
+ ///
+ public sealed class TFQMR : IIterativeSolver
+ {
+ ///
+ /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax
+ ///
+ /// Instance of the A.
+ /// Residual values in .
+ /// Instance of the x.
+ /// Instance of the b.
+ static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+ {
+ // -Ax = residual
+ matrix.Multiply(x, residual);
+ residual.Multiply(-1, residual);
+
+ // residual + b
+ residual.Add(b, residual);
+ }
+
+ ///
+ /// Is even?
+ ///
+ /// Number to check
+ /// true if even, otherwise false
+ static bool IsEven(int number)
+ {
+ return number % 2 == 0;
+ }
+
+ ///
+ /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
+ /// solution vector and x is the unknown vector.
+ ///
+ /// The coefficient matrix, A.
+ /// The solution vector, b
+ /// The result vector, x
+ /// The iterator to use to control when to stop iterating.
+ /// The preconditioner to use for approximations.
+ public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator, IPreconditioner preconditioner)
+ {
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException("Matrix must be square.", nameof(matrix));
+ }
+
+ if (input.Count != matrix.RowCount || result.Count != input.Count)
+ {
+ throw Matrix.DimensionsDontMatch(matrix, input, result);
+ }
+
+ if (iterator == null)
+ {
+ iterator = new Iterator();
+ }
+
+ if (preconditioner == null)
+ {
+ preconditioner = new UnitPreconditioner();
+ }
+
+ preconditioner.Initialize(matrix);
+
+ var d = new DenseVector(input.Count);
+ var r = DenseVector.OfVector(input);
+
+ var uodd = new DenseVector(input.Count);
+ var ueven = new DenseVector(input.Count);
+
+ var v = new DenseVector(input.Count);
+ var pseudoResiduals = DenseVector.OfVector(input);
+
+ var x = new DenseVector(input.Count);
+ var yodd = new DenseVector(input.Count);
+ var yeven = DenseVector.OfVector(input);
+
+ // Temp vectors
+ var temp = new DenseVector(input.Count);
+ var temp1 = new DenseVector(input.Count);
+ var temp2 = new DenseVector(input.Count);
+
+ // Define the scalars
+ Complex alpha = 0;
+ Complex eta = 0;
+ double theta = 0;
+
+ // Initialize
+ var tau = input.L2Norm();
+ Complex rho = tau*tau;
+
+ // Calculate the initial values for v
+ // M temp = yEven
+ preconditioner.Approximate(yeven, temp);
+
+ // v = A temp
+ matrix.Multiply(temp, v);
+
+ // Set uOdd
+ v.CopyTo(ueven);
+
+ // Start the iteration
+ var iterationNumber = 0;
+ while (iterator.DetermineStatus(iterationNumber, result, input, pseudoResiduals) == IterationStatus.Continue)
+ {
+ // First part of the step, the even bit
+ if (IsEven(iterationNumber))
+ {
+ // sigma = (v, r)
+ var sigma = r.ConjugateDotProduct(v);
+ if (sigma.Real.AlmostEqualNumbersBetween(0, 1) && sigma.Imaginary.AlmostEqualNumbersBetween(0, 1))
+ {
+ // FAIL HERE
+ iterator.Cancel();
+ break;
+ }
+
+ // alpha = rho / sigma
+ alpha = rho/sigma;
+
+ // yOdd = yEven - alpha * v
+ v.Multiply(-alpha, temp1);
+ yeven.Add(temp1, yodd);
+
+ // Solve M temp = yOdd
+ preconditioner.Approximate(yodd, temp);
+
+ // uOdd = A temp
+ matrix.Multiply(temp, uodd);
+ }
+
+ // The intermediate step which is equal for both even and
+ // odd iteration steps.
+ // Select the correct vector
+ var uinternal = IsEven(iterationNumber) ? ueven : uodd;
+ var yinternal = IsEven(iterationNumber) ? yeven : yodd;
+
+ // pseudoResiduals = pseudoResiduals - alpha * uOdd
+ uinternal.Multiply(-alpha, temp1);
+ pseudoResiduals.Add(temp1, temp2);
+ temp2.CopyTo(pseudoResiduals);
+
+ // d = yOdd + theta * theta * eta / alpha * d
+ d.Multiply(theta*theta*eta/alpha, temp);
+ yinternal.Add(temp, d);
+
+ // theta = ||pseudoResiduals||_2 / tau
+ theta = pseudoResiduals.L2Norm()/tau;
+ var c = 1/Math.Sqrt(1 + (theta*theta));
+
+ // tau = tau * theta * c
+ tau *= theta*c;
+
+ // eta = c^2 * alpha
+ eta = c*c*alpha;
+
+ // x = x + eta * d
+ d.Multiply(eta, temp1);
+ x.Add(temp1, temp2);
+ temp2.CopyTo(x);
+
+ // Check convergence and see if we can bail
+ if (iterator.DetermineStatus(iterationNumber, result, input, pseudoResiduals) != IterationStatus.Continue)
+ {
+ // Calculate the real values
+ preconditioner.Approximate(x, result);
+
+ // Calculate the true residual. Use the temp vector for that
+ // so that we don't pollute the pseudoResidual vector for no
+ // good reason.
+ CalculateTrueResidual(matrix, temp, result, input);
+
+ // Now recheck the convergence
+ if (iterator.DetermineStatus(iterationNumber, result, input, temp) != IterationStatus.Continue)
+ {
+ // We're all good now.
+ return;
+ }
+ }
+
+ // The odd step
+ if (!IsEven(iterationNumber))
+ {
+ if (rho.Real.AlmostEqualNumbersBetween(0, 1) && rho.Imaginary.AlmostEqualNumbersBetween(0, 1))
+ {
+ // FAIL HERE
+ iterator.Cancel();
+ break;
+ }
+
+ var rhoNew = r.ConjugateDotProduct(pseudoResiduals);
+ var beta = rhoNew/rho;
+
+ // Update rho for the next loop
+ rho = rhoNew;
+
+ // yOdd = pseudoResiduals + beta * yOdd
+ yodd.Multiply(beta, temp1);
+ pseudoResiduals.Add(temp1, yeven);
+
+ // Solve M temp = yOdd
+ preconditioner.Approximate(yeven, temp);
+
+ // uOdd = A temp
+ matrix.Multiply(temp, ueven);
+
+ // v = uEven + beta * (uOdd + beta * v)
+ v.Multiply(beta, temp1);
+ uodd.Add(temp1, temp);
+
+ temp.Multiply(beta, temp1);
+ ueven.Add(temp1, v);
+ }
+
+ // Calculate the real values
+ preconditioner.Approximate(x, result);
+
+ iterationNumber++;
+ }
+ }
+ }
+}
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs.meta
new file mode 100644
index 00000000..5bfd8b0d
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: 60a2d257264b54a44b3dc7d2e5971550
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
new file mode 100644
index 00000000..5ab1f838
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
@@ -0,0 +1,1543 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// 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 System.Collections.Generic;
+using System.Diagnostics;
+using System.Linq;
+using MathNet.Numerics.LinearAlgebra.Storage;
+using MathNet.Numerics.Providers.LinearAlgebra;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// A Matrix with sparse storage, intended for very large matrices where most of the cells are zero.
+ /// The underlying storage scheme is 3-array compressed-sparse-row (CSR) Format.
+ /// Wikipedia - CSR.
+ ///
+ [Serializable]
+ [DebuggerDisplay("SparseMatrix {RowCount}x{ColumnCount}-Complex {NonZerosCount}-NonZero")]
+ public class SparseMatrix : Matrix
+ {
+ readonly SparseCompressedRowMatrixStorage _storage;
+
+ ///
+ /// Gets the number of non zero elements in the matrix.
+ ///
+ /// The number of non zero elements.
+ public int NonZerosCount => _storage.ValueCount;
+
+ ///
+ /// Create a new sparse matrix straight from an initialized matrix storage instance.
+ /// The storage is used directly without copying.
+ /// Intended for advanced scenarios where you're working directly with
+ /// storage for performance or interop reasons.
+ ///
+ public SparseMatrix(SparseCompressedRowMatrixStorage storage)
+ : base(storage)
+ {
+ _storage = storage;
+ }
+
+ ///
+ /// Create a new square sparse matrix with the given number of rows and columns.
+ /// All cells of the matrix will be initialized to zero.
+ ///
+ /// If the order is less than one.
+ public SparseMatrix(int order)
+ : this(order, order)
+ {
+ }
+
+ ///
+ /// Create a new sparse matrix with the given number of rows and columns.
+ /// All cells of the matrix will be initialized to zero.
+ ///
+ /// If the row or column count is less than one.
+ public SparseMatrix(int rows, int columns)
+ : this(new SparseCompressedRowMatrixStorage(rows, columns))
+ {
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given other matrix.
+ /// This new matrix will be independent from the other matrix.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfMatrix(Matrix matrix)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfMatrix(matrix.Storage));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given two-dimensional array.
+ /// This new matrix will be independent from the provided array.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfArray(Complex[,] array)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfArray(array));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given indexed enumerable.
+ /// Keys must be provided at most once, zero is assumed if a key is omitted.
+ /// This new matrix will be independent from the enumerable.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfIndexed(int rows, int columns, IEnumerable> enumerable)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfIndexedEnumerable(rows, columns, enumerable));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given enumerable.
+ /// The enumerable is assumed to be in row-major order (row by row).
+ /// This new matrix will be independent from the enumerable.
+ /// A new memory block will be allocated for storing the vector.
+ ///
+ ///
+ public static SparseMatrix OfRowMajor(int rows, int columns, IEnumerable rowMajor)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfRowMajorEnumerable(rows, columns, rowMajor));
+ }
+
+ ///
+ /// Create a new sparse matrix with the given number of rows and columns as a copy of the given array.
+ /// The array is assumed to be in column-major order (column by column).
+ /// This new matrix will be independent from the provided array.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ ///
+ public static SparseMatrix OfColumnMajor(int rows, int columns, IList columnMajor)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfColumnMajorList(rows, columns, columnMajor));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given enumerable of enumerable columns.
+ /// Each enumerable in the master enumerable specifies a column.
+ /// This new matrix will be independent from the enumerables.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfColumns(IEnumerable> data)
+ {
+ return OfColumnArrays(data.Select(v => v.ToArray()).ToArray());
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given enumerable of enumerable columns.
+ /// Each enumerable in the master enumerable specifies a column.
+ /// This new matrix will be independent from the enumerables.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfColumns(int rows, int columns, IEnumerable> data)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfColumnEnumerables(rows, columns, data));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given column arrays.
+ /// This new matrix will be independent from the arrays.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfColumnArrays(params Complex[][] columns)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfColumnArrays(columns));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given column arrays.
+ /// This new matrix will be independent from the arrays.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfColumnArrays(IEnumerable columns)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfColumnArrays((columns as Complex[][]) ?? columns.ToArray()));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given column vectors.
+ /// This new matrix will be independent from the vectors.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfColumnVectors(params Vector[] columns)
+ {
+ var storage = new VectorStorage[columns.Length];
+ for (int i = 0; i < columns.Length; i++)
+ {
+ storage[i] = columns[i].Storage;
+ }
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfColumnVectors(storage));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given column vectors.
+ /// This new matrix will be independent from the vectors.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfColumnVectors(IEnumerable> columns)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfColumnVectors(columns.Select(c => c.Storage).ToArray()));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given enumerable of enumerable rows.
+ /// Each enumerable in the master enumerable specifies a row.
+ /// This new matrix will be independent from the enumerables.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfRows(IEnumerable> data)
+ {
+ return OfRowArrays(data.Select(v => v.ToArray()).ToArray());
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given enumerable of enumerable rows.
+ /// Each enumerable in the master enumerable specifies a row.
+ /// This new matrix will be independent from the enumerables.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfRows(int rows, int columns, IEnumerable> data)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfRowEnumerables(rows, columns, data));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given row arrays.
+ /// This new matrix will be independent from the arrays.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfRowArrays(params Complex[][] rows)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfRowArrays(rows));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given row arrays.
+ /// This new matrix will be independent from the arrays.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfRowArrays(IEnumerable rows)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfRowArrays((rows as Complex[][]) ?? rows.ToArray()));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given row vectors.
+ /// This new matrix will be independent from the vectors.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfRowVectors(params Vector[] rows)
+ {
+ var storage = new VectorStorage[rows.Length];
+ for (int i = 0; i < rows.Length; i++)
+ {
+ storage[i] = rows[i].Storage;
+ }
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfRowVectors(storage));
+ }
+
+ ///
+ /// Create a new sparse matrix as a copy of the given row vectors.
+ /// This new matrix will be independent from the vectors.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfRowVectors(IEnumerable> rows)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfRowVectors(rows.Select(r => r.Storage).ToArray()));
+ }
+
+ ///
+ /// Create a new sparse matrix with the diagonal as a copy of the given vector.
+ /// This new matrix will be independent from the vector.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfDiagonalVector(Vector diagonal)
+ {
+ var m = new SparseMatrix(diagonal.Count, diagonal.Count);
+ m.SetDiagonal(diagonal);
+ return m;
+ }
+
+ ///
+ /// Create a new sparse matrix with the diagonal as a copy of the given vector.
+ /// This new matrix will be independent from the vector.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfDiagonalVector(int rows, int columns, Vector diagonal)
+ {
+ var m = new SparseMatrix(rows, columns);
+ m.SetDiagonal(diagonal);
+ return m;
+ }
+
+ ///
+ /// Create a new sparse matrix with the diagonal as a copy of the given array.
+ /// This new matrix will be independent from the array.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfDiagonalArray(Complex[] diagonal)
+ {
+ var m = new SparseMatrix(diagonal.Length, diagonal.Length);
+ m.SetDiagonal(diagonal);
+ return m;
+ }
+
+ ///
+ /// Create a new sparse matrix with the diagonal as a copy of the given array.
+ /// This new matrix will be independent from the array.
+ /// A new memory block will be allocated for storing the matrix.
+ ///
+ public static SparseMatrix OfDiagonalArray(int rows, int columns, Complex[] diagonal)
+ {
+ var m = new SparseMatrix(rows, columns);
+ m.SetDiagonal(diagonal);
+ return m;
+ }
+
+ ///
+ /// Create a new sparse matrix and initialize each value to the same provided value.
+ ///
+ public static SparseMatrix Create(int rows, int columns, Complex value)
+ {
+ if (value == Complex.Zero) return new SparseMatrix(rows, columns);
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfValue(rows, columns, value));
+ }
+
+ ///
+ /// Create a new sparse matrix and initialize each value using the provided init function.
+ ///
+ public static SparseMatrix Create(int rows, int columns, Func init)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfInit(rows, columns, init));
+ }
+
+ ///
+ /// Create a new diagonal sparse matrix and initialize each diagonal value to the same provided value.
+ ///
+ public static SparseMatrix CreateDiagonal(int rows, int columns, Complex value)
+ {
+ if (value == Complex.Zero) return new SparseMatrix(rows, columns);
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfDiagonalInit(rows, columns, i => value));
+ }
+
+ ///
+ /// Create a new diagonal sparse matrix and initialize each diagonal value using the provided init function.
+ ///
+ public static SparseMatrix CreateDiagonal(int rows, int columns, Func init)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfDiagonalInit(rows, columns, init));
+ }
+
+ ///
+ /// Create a new square sparse identity matrix where each diagonal value is set to One.
+ ///
+ public static SparseMatrix CreateIdentity(int order)
+ {
+ return new SparseMatrix(SparseCompressedRowMatrixStorage.OfDiagonalInit(order, order, i => One));
+ }
+
+ ///
+ /// Returns a new matrix containing the lower triangle of this matrix.
+ ///
+ /// The lower triangle of this matrix.
+ public override Matrix LowerTriangle()
+ {
+ var result = Build.SameAs(this);
+ LowerTriangleImpl(result);
+ return result;
+ }
+
+ ///
+ /// Puts the lower triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ /// If is .
+ /// If the result matrix's dimensions are not the same as this matrix.
+ public override void LowerTriangle(Matrix result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException(nameof(result));
+ }
+
+ if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
+ {
+ throw DimensionsDontMatch(this, result);
+ }
+
+ if (ReferenceEquals(this, result))
+ {
+ var tmp = Build.SameAs(result);
+ LowerTriangle(tmp);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ result.Clear();
+ LowerTriangleImpl(result);
+ }
+ }
+
+ ///
+ /// Puts the lower triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ private void LowerTriangleImpl(Matrix result)
+ {
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var row = 0; row < result.RowCount; row++)
+ {
+ var endIndex = rowPointers[row + 1];
+ for (var j = rowPointers[row]; j < endIndex; j++)
+ {
+ if (row >= columnIndices[j])
+ {
+ result.At(row, columnIndices[j], values[j]);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Returns a new matrix containing the upper triangle of this matrix.
+ ///
+ /// The upper triangle of this matrix.
+ public override Matrix UpperTriangle()
+ {
+ var result = Build.SameAs(this);
+ UpperTriangleImpl(result);
+ return result;
+ }
+
+ ///
+ /// Puts the upper triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ /// If is .
+ /// If the result matrix's dimensions are not the same as this matrix.
+ public override void UpperTriangle(Matrix result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException(nameof(result));
+ }
+
+ if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
+ {
+ throw DimensionsDontMatch(this, result);
+ }
+
+ if (ReferenceEquals(this, result))
+ {
+ var tmp = Build.SameAs(result);
+ UpperTriangle(tmp);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ result.Clear();
+ UpperTriangleImpl(result);
+ }
+ }
+
+ ///
+ /// Puts the upper triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ private void UpperTriangleImpl(Matrix result)
+ {
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var row = 0; row < result.RowCount; row++)
+ {
+ var endIndex = rowPointers[row + 1];
+ for (var j = rowPointers[row]; j < endIndex; j++)
+ {
+ if (row <= columnIndices[j])
+ {
+ result.At(row, columnIndices[j], values[j]);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Returns a new matrix containing the lower triangle of this matrix. The new matrix
+ /// does not contain the diagonal elements of this matrix.
+ ///
+ /// The lower triangle of this matrix.
+ public override Matrix StrictlyLowerTriangle()
+ {
+ var result = Build.SameAs(this);
+ StrictlyLowerTriangleImpl(result);
+ return result;
+ }
+
+ ///
+ /// Puts the strictly lower triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ /// If is .
+ /// If the result matrix's dimensions are not the same as this matrix.
+ public override void StrictlyLowerTriangle(Matrix result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException(nameof(result));
+ }
+
+ if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
+ {
+ throw DimensionsDontMatch(this, result);
+ }
+
+ if (ReferenceEquals(this, result))
+ {
+ var tmp = Build.SameAs(result);
+ StrictlyLowerTriangle(tmp);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ result.Clear();
+ StrictlyLowerTriangleImpl(result);
+ }
+ }
+
+ ///
+ /// Puts the strictly lower triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ private void StrictlyLowerTriangleImpl(Matrix result)
+ {
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var row = 0; row < result.RowCount; row++)
+ {
+ var endIndex = rowPointers[row + 1];
+ for (var j = rowPointers[row]; j < endIndex; j++)
+ {
+ if (row > columnIndices[j])
+ {
+ result.At(row, columnIndices[j], values[j]);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Returns a new matrix containing the upper triangle of this matrix. The new matrix
+ /// does not contain the diagonal elements of this matrix.
+ ///
+ /// The upper triangle of this matrix.
+ public override Matrix StrictlyUpperTriangle()
+ {
+ var result = Build.SameAs(this);
+ StrictlyUpperTriangleImpl(result);
+ return result;
+ }
+
+ ///
+ /// Puts the strictly upper triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ /// If is .
+ /// If the result matrix's dimensions are not the same as this matrix.
+ public override void StrictlyUpperTriangle(Matrix result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException(nameof(result));
+ }
+
+ if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
+ {
+ throw DimensionsDontMatch(this, result);
+ }
+
+ if (ReferenceEquals(this, result))
+ {
+ var tmp = Build.SameAs(result);
+ StrictlyUpperTriangle(tmp);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ result.Clear();
+ StrictlyUpperTriangleImpl(result);
+ }
+ }
+
+ ///
+ /// Puts the strictly upper triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ private void StrictlyUpperTriangleImpl(Matrix result)
+ {
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var row = 0; row < result.RowCount; row++)
+ {
+ var endIndex = rowPointers[row + 1];
+ for (var j = rowPointers[row]; j < endIndex; j++)
+ {
+ if (row < columnIndices[j])
+ {
+ result.At(row, columnIndices[j], values[j]);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Negate each element of this matrix and place the results into the result matrix.
+ ///
+ /// The result of the negation.
+ protected override void DoNegate(Matrix result)
+ {
+ CopyTo(result);
+ DoMultiply(-1, result);
+ }
+
+ /// Calculates the induced infinity norm of this matrix.
+ /// The maximum absolute row sum of the matrix.
+ public override double InfinityNorm()
+ {
+ var rowPointers = _storage.RowPointers;
+ var values = _storage.Values;
+ var norm = 0d;
+ for (var i = 0; i < RowCount; i++)
+ {
+ var startIndex = rowPointers[i];
+ var endIndex = rowPointers[i + 1];
+
+ if (startIndex == endIndex)
+ {
+ continue;
+ }
+
+ var s = 0d;
+ for (var j = startIndex; j < endIndex; j++)
+ {
+ s += values[j].Magnitude;
+ }
+ norm = Math.Max(norm, s);
+ }
+ return norm;
+ }
+
+ /// Calculates the entry-wise Frobenius norm of this matrix.
+ /// The square root of the sum of the squared values.
+ public override double FrobeniusNorm()
+ {
+ var aat = (SparseCompressedRowMatrixStorage) (this*ConjugateTranspose()).Storage;
+ var norm = 0d;
+ for (var i = 0; i < aat.RowCount; i++)
+ {
+ var startIndex = aat.RowPointers[i];
+ var endIndex = aat.RowPointers[i + 1];
+
+ if (startIndex == endIndex)
+ {
+ continue;
+ }
+
+ for (var j = startIndex; j < endIndex; j++)
+ {
+ if (i == aat.ColumnIndices[j])
+ {
+ norm += aat.Values[j].Magnitude;
+ }
+ }
+ }
+ return Math.Sqrt(norm);
+ }
+
+ ///
+ /// Adds another matrix to this matrix.
+ ///
+ /// The matrix to add to this matrix.
+ /// The matrix to store the result of the addition.
+ /// If the other matrix is .
+ /// If the two matrices don't have the same dimensions.
+ protected override void DoAdd(Matrix other, Matrix result)
+ {
+ if (other is SparseMatrix sparseOther && result is SparseMatrix sparseResult)
+ {
+ if (ReferenceEquals(this, other))
+ {
+ if (!ReferenceEquals(this, result))
+ {
+ CopyTo(result);
+ }
+
+ LinearAlgebraControl.Provider.ScaleArray(2.0, sparseResult._storage.Values, sparseResult._storage.Values);
+ return;
+ }
+
+ SparseMatrix left;
+
+ if (ReferenceEquals(sparseOther, sparseResult))
+ {
+ left = this;
+ }
+ else if (ReferenceEquals(this, sparseResult))
+ {
+ left = sparseOther;
+ }
+ else
+ {
+ CopyTo(sparseResult);
+ left = sparseOther;
+ }
+
+ var leftStorage = left._storage;
+ for (var i = 0; i < leftStorage.RowCount; i++)
+ {
+ var endIndex = leftStorage.RowPointers[i + 1];
+ for (var j = leftStorage.RowPointers[i]; j < endIndex; j++)
+ {
+ var columnIndex = leftStorage.ColumnIndices[j];
+ var resVal = leftStorage.Values[j] + result.At(i, columnIndex);
+ result.At(i, columnIndex, resVal);
+ }
+ }
+ }
+ else
+ {
+ base.DoAdd(other, result);
+ }
+ }
+
+ ///
+ /// Subtracts another matrix from this matrix.
+ ///
+ /// The matrix to subtract to this matrix.
+ /// The matrix to store the result of subtraction.
+ /// If the other matrix is .
+ /// If the two matrices don't have the same dimensions.
+ protected override void DoSubtract(Matrix other, Matrix result)
+ {
+ var sparseOther = other as SparseMatrix;
+ var sparseResult = result as SparseMatrix;
+
+ if (sparseOther == null || sparseResult == null)
+ {
+ base.DoSubtract(other, result);
+ return;
+ }
+
+ if (ReferenceEquals(this, other))
+ {
+ result.Clear();
+ return;
+ }
+
+ var otherStorage = sparseOther._storage;
+
+ if (ReferenceEquals(this, sparseResult))
+ {
+ for (var i = 0; i < otherStorage.RowCount; i++)
+ {
+ var endIndex = otherStorage.RowPointers[i + 1];
+ for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
+ {
+ var columnIndex = otherStorage.ColumnIndices[j];
+ var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
+ result.At(i, columnIndex, resVal);
+ }
+ }
+ }
+ else
+ {
+ if (!ReferenceEquals(sparseOther, sparseResult))
+ {
+ sparseOther.CopyTo(sparseResult);
+ }
+
+ sparseResult.Negate(sparseResult);
+
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var i = 0; i < RowCount; i++)
+ {
+ var endIndex = rowPointers[i + 1];
+ for (var j = rowPointers[i]; j < endIndex; j++)
+ {
+ var columnIndex = columnIndices[j];
+ var resVal = sparseResult.At(i, columnIndex) + values[j];
+ result.At(i, columnIndex, resVal);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Multiplies each element of the matrix by a scalar and places results into the result matrix.
+ ///
+ /// The scalar to multiply the matrix with.
+ /// The matrix to store the result of the multiplication.
+ protected override void DoMultiply(Complex scalar, Matrix result)
+ {
+ if (scalar == 1.0)
+ {
+ CopyTo(result);
+ return;
+ }
+
+ if (scalar == 0.0 || NonZerosCount == 0)
+ {
+ result.Clear();
+ return;
+ }
+
+ if (result is SparseMatrix sparseResult)
+ {
+ if (!ReferenceEquals(this, result))
+ {
+ CopyTo(sparseResult);
+ }
+
+ LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
+ }
+ else
+ {
+ result.Clear();
+
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var row = 0; row < RowCount; row++)
+ {
+ var start = rowPointers[row];
+ var end = rowPointers[row + 1];
+
+ if (start == end)
+ {
+ continue;
+ }
+
+ for (var index = start; index < end; index++)
+ {
+ var column = columnIndices[index];
+ result.At(row, column, values[index] * scalar);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Multiplies this matrix with another matrix and places the results into the result matrix.
+ ///
+ /// The matrix to multiply with.
+ /// The result of the multiplication.
+ protected override void DoMultiply(Matrix other, Matrix result)
+ {
+ var sparseOther = other as SparseMatrix;
+ var sparseResult = result as SparseMatrix;
+ if (sparseOther != null && sparseResult != null)
+ {
+ DoMultiplySparse(sparseOther, sparseResult);
+ return;
+ }
+
+ if (other.Storage is DiagonalMatrixStorage diagonalOther && sparseResult != null)
+ {
+ var diagonal = diagonalOther.Data;
+ if (other.ColumnCount == other.RowCount)
+ {
+ Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], Zeros.AllowSkip, ExistingData.Clear);
+ }
+ else
+ {
+ result.Storage.Clear();
+ Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, Math.Min(ColumnCount, other.ColumnCount), Zeros.AllowSkip, ExistingData.AssumeZeros);
+ }
+ return;
+ }
+
+ result.Clear();
+
+
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+ if (other.Storage is DenseColumnMajorMatrixStorage denseOther)
+ {
+ // in this case we can directly address the underlying data-array
+ for (var row = 0; row < RowCount; row++)
+ {
+ var startIndex = rowPointers[row];
+ var endIndex = rowPointers[row + 1];
+
+ if (startIndex == endIndex)
+ {
+ continue;
+ }
+
+ for (var column = 0; column < other.ColumnCount; column++)
+ {
+ int otherColumnStartPosition = column * other.RowCount;
+ var sum = Complex.Zero;
+ for (var index = startIndex; index < endIndex; index++)
+ {
+ sum += values[index] * denseOther.Data[otherColumnStartPosition + columnIndices[index]];
+ }
+
+ result.At(row, column, sum);
+ }
+ }
+ return;
+ }
+
+ var columnVector = new DenseVector(other.RowCount);
+ for (var row = 0; row < RowCount; row++)
+ {
+ var startIndex = rowPointers[row];
+ var endIndex = rowPointers[row + 1];
+
+ if (startIndex == endIndex)
+ {
+ continue;
+ }
+
+ for (var column = 0; column < other.ColumnCount; column++)
+ {
+ // Multiply row of matrix A on column of matrix B
+ other.Column(column, columnVector);
+
+ var sum = Complex.Zero;
+ for (var index = startIndex; index < endIndex; index++)
+ {
+ sum += values[index] * columnVector[columnIndices[index]];
+ }
+
+ result.At(row, column, sum);
+ }
+ }
+ }
+
+ void DoMultiplySparse(SparseMatrix other, SparseMatrix result)
+ {
+ result.Clear();
+
+ var ax = _storage.Values;
+ var ap = _storage.RowPointers;
+ var ai = _storage.ColumnIndices;
+
+ var bx = other._storage.Values;
+ var bp = other._storage.RowPointers;
+ var bi = other._storage.ColumnIndices;
+
+ int rows = RowCount;
+ int cols = other.ColumnCount;
+
+ int[] cp = result._storage.RowPointers;
+
+ var marker = new int[cols];
+ for (int ib = 0; ib < cols; ib++)
+ {
+ marker[ib] = -1;
+ }
+
+ int count = 0;
+ for (int i = 0; i < rows; i++)
+ {
+ // For each row of A
+ for (int j = ap[i]; j < ap[i + 1]; j++)
+ {
+ // Row number to be added
+ int a = ai[j];
+ for (int k = bp[a]; k < bp[a + 1]; k++)
+ {
+ int b = bi[k];
+ if (marker[b] != i)
+ {
+ marker[b] = i;
+ count++;
+ }
+ }
+ }
+
+ // Record non-zero count.
+ cp[i + 1] = count;
+ }
+
+ var ci = new int[count];
+ var cx = new Complex[count];
+
+ for (int ib = 0; ib < cols; ib++)
+ {
+ marker[ib] = -1;
+ }
+
+ count = 0;
+ for (int i = 0; i < rows; i++)
+ {
+ int rowStart = cp[i];
+ for (int j = ap[i]; j < ap[i + 1]; j++)
+ {
+ int a = ai[j];
+ Complex aEntry = ax[j];
+ for (int k = bp[a]; k < bp[a + 1]; k++)
+ {
+ int b = bi[k];
+ Complex bEntry = bx[k];
+ if (marker[b] < rowStart)
+ {
+ marker[b] = count;
+ ci[marker[b]] = b;
+ cx[marker[b]] = aEntry * bEntry;
+ count++;
+ }
+ else
+ {
+ cx[marker[b]] += aEntry * bEntry;
+ }
+ }
+ }
+ }
+
+ result._storage.Values = cx;
+ result._storage.ColumnIndices = ci;
+ result._storage.Normalize();
+ }
+
+ ///
+ /// Multiplies this matrix with a vector and places the results into the result vector.
+ ///
+ /// The vector to multiply with.
+ /// The result of the multiplication.
+ protected override void DoMultiply(Vector rightSide, Vector result)
+ {
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var row = 0; row < RowCount; row++)
+ {
+ var startIndex = rowPointers[row];
+ var endIndex = rowPointers[row + 1];
+
+ if (startIndex == endIndex)
+ {
+ continue;
+ }
+
+ var sum = Complex.Zero;
+ for (var index = startIndex; index < endIndex; index++)
+ {
+ sum += values[index] * rightSide[columnIndices[index]];
+ }
+
+ result[row] = sum;
+ }
+ }
+
+ ///
+ /// Multiplies this matrix with transpose of another matrix and places the results into the result matrix.
+ ///
+ /// The matrix to multiply with.
+ /// The result of the multiplication.
+ protected override void DoTransposeAndMultiply(Matrix other, Matrix result)
+ {
+ if (other is SparseMatrix otherSparse && result is SparseMatrix resultSparse)
+ {
+ resultSparse.Clear();
+
+ var rowPointers = _storage.RowPointers;
+ var values = _storage.Values;
+
+ var otherStorage = otherSparse._storage;
+
+ for (var j = 0; j < RowCount; j++)
+ {
+ var startIndexOther = otherStorage.RowPointers[j];
+ var endIndexOther = otherStorage.RowPointers[j + 1];
+
+ if (startIndexOther == endIndexOther)
+ {
+ continue;
+ }
+
+ for (var i = 0; i < RowCount; i++)
+ {
+ // Multiply row of matrix A on row of matrix B
+
+ var startIndexThis = rowPointers[i];
+ var endIndexThis = rowPointers[i + 1];
+
+ if (startIndexThis == endIndexThis)
+ {
+ continue;
+ }
+
+ var sum = Complex.Zero;
+ for (var index = startIndexOther; index < endIndexOther; index++)
+ {
+ var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
+ if (ind >= 0)
+ {
+ sum += otherStorage.Values[index] * values[ind];
+ }
+ }
+
+ resultSparse._storage.At(i, j, sum + result.At(i, j));
+ }
+ }
+ }
+ else
+ {
+ base.DoTransposeAndMultiply(other, result);
+ }
+ }
+
+ ///
+ /// Multiplies the transpose of this matrix with a vector and places the results into the result vector.
+ ///
+ /// The vector to multiply with.
+ /// The result of the multiplication.
+ protected override void DoTransposeThisAndMultiply(Vector rightSide, Vector result)
+ {
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var row = 0; row < RowCount; row++)
+ {
+ var startIndex = rowPointers[row];
+ var endIndex = rowPointers[row + 1];
+
+ if (startIndex == endIndex)
+ {
+ continue;
+ }
+
+ var rightSideValue = rightSide[row];
+ for (var index = startIndex; index < endIndex; index++)
+ {
+ result[columnIndices[index]] += values[index] * rightSideValue;
+ }
+ }
+ }
+
+ ///
+ /// Pointwise multiplies this matrix with another matrix and stores the result into the result matrix.
+ ///
+ /// The matrix to pointwise multiply with this one.
+ /// The matrix to store the result of the pointwise multiplication.
+ protected override void DoPointwiseMultiply(Matrix other, Matrix result)
+ {
+ result.Clear();
+
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var i = 0; i < RowCount; i++)
+ {
+ var endIndex = rowPointers[i + 1];
+ for (var j = rowPointers[i]; j < endIndex; j++)
+ {
+ var resVal = values[j]*other.At(i, columnIndices[j]);
+ if (!resVal.IsZero())
+ {
+ result.At(i, columnIndices[j], resVal);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Pointwise divide this matrix by another matrix and stores the result into the result matrix.
+ ///
+ /// The matrix to pointwise divide this one by.
+ /// The matrix to store the result of the pointwise division.
+ protected override void DoPointwiseDivide(Matrix divisor, Matrix result)
+ {
+ result.Clear();
+
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var i = 0; i < RowCount; i++)
+ {
+ var endIndex = rowPointers[i + 1];
+ for (var j = rowPointers[i]; j < endIndex; j++)
+ {
+ if (!values[j].IsZero())
+ {
+ result.At(i, columnIndices[j], values[j]/divisor.At(i, columnIndices[j]));
+ }
+ }
+ }
+ }
+
+ public override void KroneckerProduct(Matrix other, Matrix result)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException(nameof(other));
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException(nameof(result));
+ }
+
+ if (result.RowCount != (RowCount*other.RowCount) || result.ColumnCount != (ColumnCount*other.ColumnCount))
+ {
+ throw DimensionsDontMatch(this, other, result);
+ }
+
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var i = 0; i < RowCount; i++)
+ {
+ var endIndex = rowPointers[i + 1];
+ for (var j = rowPointers[i]; j < endIndex; j++)
+ {
+ if (!values[j].IsZero())
+ {
+ result.SetSubMatrix(i*other.RowCount, other.RowCount, columnIndices[j]*other.ColumnCount, other.ColumnCount, values[j]*other);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Evaluates whether this matrix is symmetric.
+ ///
+ public override bool IsSymmetric()
+ {
+ if (RowCount != ColumnCount)
+ {
+ return false;
+ }
+
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var row = 0; row < RowCount; row++)
+ {
+ var start = rowPointers[row];
+ var end = rowPointers[row + 1];
+
+ if (start == end)
+ {
+ continue;
+ }
+
+ for (var index = start; index < end; index++)
+ {
+ var column = columnIndices[index];
+ if (!values[index].Equals(At(column, row)))
+ {
+ return false;
+ }
+ }
+ }
+
+ return true;
+ }
+
+ ///
+ /// Evaluates whether this matrix is Hermitian (conjugate symmetric).
+ ///
+ public override bool IsHermitian()
+ {
+ if (RowCount != ColumnCount)
+ {
+ return false;
+ }
+
+ var rowPointers = _storage.RowPointers;
+ var columnIndices = _storage.ColumnIndices;
+ var values = _storage.Values;
+
+ for (var row = 0; row < RowCount; row++)
+ {
+ var start = rowPointers[row];
+ var end = rowPointers[row + 1];
+
+ if (start == end)
+ {
+ continue;
+ }
+
+ for (var index = start; index < end; index++)
+ {
+ var column = columnIndices[index];
+ if (!values[index].Equals(At(column, row).Conjugate()))
+ {
+ return false;
+ }
+ }
+ }
+
+ return true;
+ }
+
+ ///
+ /// Adds two matrices together and returns the results.
+ ///
+ /// This operator will allocate new memory for the result. It will
+ /// choose the representation of either or depending on which
+ /// is denser.
+ /// The left matrix to add.
+ /// The right matrix to add.
+ /// The result of the addition.
+ /// If and don't have the same dimensions.
+ /// If or is .
+ public static SparseMatrix operator +(SparseMatrix leftSide, SparseMatrix rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException(nameof(rightSide));
+ }
+
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ if (leftSide.RowCount != rightSide.RowCount || leftSide.ColumnCount != rightSide.ColumnCount)
+ {
+ throw DimensionsDontMatch(leftSide, rightSide);
+ }
+
+ return (SparseMatrix)leftSide.Add(rightSide);
+ }
+
+ ///
+ /// Returns a Matrix containing the same values of .
+ ///
+ /// The matrix to get the values from.
+ /// A matrix containing a the same values as .
+ /// If is .
+ public static SparseMatrix operator +(SparseMatrix rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException(nameof(rightSide));
+ }
+
+ return (SparseMatrix)rightSide.Clone();
+ }
+
+ ///
+ /// Subtracts two matrices together and returns the results.
+ ///
+ /// This operator will allocate new memory for the result. It will
+ /// choose the representation of either or depending on which
+ /// is denser.
+ /// The left matrix to subtract.
+ /// The right matrix to subtract.
+ /// The result of the addition.
+ /// If and don't have the same dimensions.
+ /// If or is .
+ public static SparseMatrix operator -(SparseMatrix leftSide, SparseMatrix rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException(nameof(rightSide));
+ }
+
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ if (leftSide.RowCount != rightSide.RowCount || leftSide.ColumnCount != rightSide.ColumnCount)
+ {
+ throw DimensionsDontMatch(leftSide, rightSide);
+ }
+
+ return (SparseMatrix)leftSide.Subtract(rightSide);
+ }
+
+ ///
+ /// Negates each element of the matrix.
+ ///
+ /// The matrix to negate.
+ /// A matrix containing the negated values.
+ /// If is .
+ public static SparseMatrix operator -(SparseMatrix rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException(nameof(rightSide));
+ }
+
+ return (SparseMatrix)rightSide.Negate();
+ }
+
+ ///
+ /// Multiplies a Matrix by a constant and returns the result.
+ ///
+ /// The matrix to multiply.
+ /// The constant to multiply the matrix by.
+ /// The result of the multiplication.
+ /// If is .
+ public static SparseMatrix operator *(SparseMatrix leftSide, Complex rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ return (SparseMatrix)leftSide.Multiply(rightSide);
+ }
+
+ ///
+ /// Multiplies a Matrix by a constant and returns the result.
+ ///
+ /// The matrix to multiply.
+ /// The constant to multiply the matrix by.
+ /// The result of the multiplication.
+ /// If is .
+ public static SparseMatrix operator *(Complex leftSide, SparseMatrix rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException(nameof(rightSide));
+ }
+
+ return (SparseMatrix)rightSide.Multiply(leftSide);
+ }
+
+ ///
+ /// Multiplies two matrices.
+ ///
+ /// This operator will allocate new memory for the result. It will
+ /// choose the representation of either or depending on which
+ /// is denser.
+ /// The left matrix to multiply.
+ /// The right matrix to multiply.
+ /// The result of multiplication.
+ /// If or is .
+ /// If the dimensions of or don't conform.
+ public static SparseMatrix operator *(SparseMatrix leftSide, SparseMatrix rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException(nameof(rightSide));
+ }
+
+ if (leftSide.ColumnCount != rightSide.RowCount)
+ {
+ throw DimensionsDontMatch(leftSide, rightSide);
+ }
+
+ return (SparseMatrix)leftSide.Multiply(rightSide);
+ }
+
+ ///
+ /// Multiplies a Matrix and a Vector.
+ ///
+ /// The matrix to multiply.
+ /// The vector to multiply.
+ /// The result of multiplication.
+ /// If or is .
+ public static SparseVector operator *(SparseMatrix leftSide, SparseVector rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ return (SparseVector)leftSide.Multiply(rightSide);
+ }
+
+ ///
+ /// Multiplies a Vector and a Matrix.
+ ///
+ /// The vector to multiply.
+ /// The matrix to multiply.
+ /// The result of multiplication.
+ /// If or is .
+ public static SparseVector operator *(SparseVector leftSide, SparseMatrix rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException(nameof(rightSide));
+ }
+
+ return (SparseVector)rightSide.LeftMultiply(leftSide);
+ }
+
+ ///
+ /// Multiplies a Matrix by a constant and returns the result.
+ ///
+ /// The matrix to multiply.
+ /// The constant to multiply the matrix by.
+ /// The result of the multiplication.
+ /// If is .
+ public static SparseMatrix operator %(SparseMatrix leftSide, Complex rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ return (SparseMatrix)leftSide.Remainder(rightSide);
+ }
+
+ public override string ToTypeString()
+ {
+ return FormattableString.Invariant($"SparseMatrix {RowCount}x{ColumnCount}-Complex {NonZerosCount / (RowCount * (double) ColumnCount):P2} Filled");
+ }
+ }
+ }
+
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseMatrix.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseMatrix.cs.meta
new file mode 100644
index 00000000..6c15965e
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseMatrix.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: d5190d79b37cfeb4abe6a0321fc026e7
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseVector.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseVector.cs
new file mode 100644
index 00000000..aaf3d691
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseVector.cs
@@ -0,0 +1,921 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2015 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 System.Collections.Generic;
+using System.Diagnostics;
+using MathNet.Numerics.LinearAlgebra.Storage;
+using MathNet.Numerics.Providers.LinearAlgebra;
+using MathNet.Numerics.Threading;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// A vector with sparse storage, intended for very large vectors where most of the cells are zero.
+ ///
+ /// The sparse vector is not thread safe.
+ [Serializable]
+ [DebuggerDisplay("SparseVector {Count}-Complex {NonZerosCount}-NonZero")]
+ public class SparseVector : Vector
+ {
+ readonly SparseVectorStorage _storage;
+
+ ///
+ /// Gets the number of non zero elements in the vector.
+ ///
+ /// The number of non zero elements.
+ public int NonZerosCount => _storage.ValueCount;
+
+ ///
+ /// Create a new sparse vector straight from an initialized vector storage instance.
+ /// The storage is used directly without copying.
+ /// Intended for advanced scenarios where you're working directly with
+ /// storage for performance or interop reasons.
+ ///
+ public SparseVector(SparseVectorStorage storage)
+ : base(storage)
+ {
+ _storage = storage;
+ }
+
+ ///
+ /// Create a new sparse vector with the given length.
+ /// All cells of the vector will be initialized to zero.
+ ///
+ /// If length is less than one.
+ public SparseVector(int length)
+ : this(new SparseVectorStorage(length))
+ {
+ }
+
+ ///
+ /// Create a new sparse vector as a copy of the given other vector.
+ /// This new vector will be independent from the other vector.
+ /// A new memory block will be allocated for storing the vector.
+ ///
+ public static SparseVector OfVector(Vector vector)
+ {
+ return new SparseVector(SparseVectorStorage.OfVector(vector.Storage));
+ }
+
+ ///
+ /// Create a new sparse vector as a copy of the given enumerable.
+ /// This new vector will be independent from the enumerable.
+ /// A new memory block will be allocated for storing the vector.
+ ///
+ public static SparseVector OfEnumerable(IEnumerable enumerable)
+ {
+ return new SparseVector(SparseVectorStorage.OfEnumerable(enumerable));
+ }
+
+ ///
+ /// Create a new sparse vector as a copy of the given indexed enumerable.
+ /// Keys must be provided at most once, zero is assumed if a key is omitted.
+ /// This new vector will be independent from the enumerable.
+ /// A new memory block will be allocated for storing the vector.
+ ///
+ public static SparseVector OfIndexedEnumerable(int length, IEnumerable> enumerable)
+ {
+ return new SparseVector(SparseVectorStorage.OfIndexedEnumerable(length, enumerable));
+ }
+
+ ///
+ /// Create a new sparse vector and initialize each value using the provided value.
+ ///
+ public static SparseVector Create(int length, Complex value)
+ {
+ return new SparseVector(SparseVectorStorage.OfValue(length, value));
+ }
+
+ ///
+ /// Create a new sparse vector and initialize each value using the provided init function.
+ ///
+ public static SparseVector Create(int length, Func init)
+ {
+ return new SparseVector(SparseVectorStorage.OfInit(length, init));
+ }
+
+ ///
+ /// Adds a scalar to each element of the vector and stores the result in the result vector.
+ /// Warning, the new 'sparse vector' with a non-zero scalar added to it will be a 100% filled
+ /// sparse vector and very inefficient. Would be better to work with a dense vector instead.
+ ///
+ ///
+ /// The scalar to add.
+ ///
+ ///
+ /// The vector to store the result of the addition.
+ ///
+ protected override void DoAdd(Complex scalar, Vector result)
+ {
+ if (scalar == Complex.Zero)
+ {
+ if (!ReferenceEquals(this, result))
+ {
+ CopyTo(result);
+ }
+
+ return;
+ }
+
+ if (ReferenceEquals(this, result))
+ {
+ //populate a new vector with the scalar
+ var vnonZeroValues = new Complex[Count];
+ var vnonZeroIndices = new int[Count];
+ for (int index = 0; index < Count; index++)
+ {
+ vnonZeroIndices[index] = index;
+ vnonZeroValues[index] = scalar;
+ }
+
+ //populate the non zero values from this
+ var indices = _storage.Indices;
+ var values = _storage.Values;
+ for (int j = 0; j < _storage.ValueCount; j++)
+ {
+ vnonZeroValues[indices[j]] = values[j] + scalar;
+ }
+
+ //assign this vectors array to the new arrays.
+ _storage.Values = vnonZeroValues;
+ _storage.Indices = vnonZeroIndices;
+ _storage.ValueCount = Count;
+ }
+ else
+ {
+ for (var index = 0; index < Count; index++)
+ {
+ result.At(index, At(index) + scalar);
+ }
+ }
+ }
+
+ ///
+ /// Adds another vector to this vector and stores the result into the result vector.
+ ///
+ ///
+ /// The vector to add to this one.
+ ///
+ ///
+ /// The vector to store the result of the addition.
+ ///
+ protected override void DoAdd(Vector other, Vector result)
+ {
+ if (other is SparseVector otherSparse && result is SparseVector resultSparse)
+ {
+ // TODO (ruegg, 2011-10-11): Options to optimize?
+ var otherStorage = otherSparse._storage;
+ if (ReferenceEquals(this, resultSparse))
+ {
+ int i = 0, j = 0;
+ while (j < otherStorage.ValueCount)
+ {
+ if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
+ {
+ var otherValue = otherStorage.Values[j];
+ if (!Complex.Zero.Equals(otherValue))
+ {
+ _storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], otherValue);
+ }
+
+ j++;
+ }
+ else if (_storage.Indices[i] == otherStorage.Indices[j])
+ {
+ // TODO: result can be zero, remove?
+ _storage.Values[i++] += otherStorage.Values[j++];
+ }
+ else
+ {
+ i++;
+ }
+ }
+ }
+ else
+ {
+ result.Clear();
+ int i = 0, j = 0, last = -1;
+ while (i < _storage.ValueCount || j < otherStorage.ValueCount)
+ {
+ if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
+ {
+ var next = _storage.Indices[i];
+ if (next != last)
+ {
+ last = next;
+ result.At(next, _storage.Values[i] + otherSparse.At(next));
+ }
+
+ i++;
+ }
+ else
+ {
+ var next = otherStorage.Indices[j];
+ if (next != last)
+ {
+ last = next;
+ result.At(next, At(next) + otherStorage.Values[j]);
+ }
+
+ j++;
+ }
+ }
+ }
+ }
+ else
+ {
+ base.DoAdd(other, result);
+ }
+ }
+
+ ///
+ /// Subtracts a scalar from each element of the vector and stores the result in the result vector.
+ ///
+ ///
+ /// The scalar to subtract.
+ ///
+ ///
+ /// The vector to store the result of the subtraction.
+ ///
+ protected override void DoSubtract(Complex scalar, Vector result)
+ {
+ DoAdd(-scalar, result);
+ }
+
+ ///
+ /// Subtracts another vector to this vector and stores the result into the result vector.
+ ///
+ ///
+ /// The vector to subtract from this one.
+ ///
+ ///
+ /// The vector to store the result of the subtraction.
+ ///
+ protected override void DoSubtract(Vector other, Vector result)
+ {
+ if (ReferenceEquals(this, other))
+ {
+ result.Clear();
+ return;
+ }
+
+ if (other is SparseVector otherSparse && result is SparseVector resultSparse)
+ {
+ // TODO (ruegg, 2011-10-11): Options to optimize?
+ var otherStorage = otherSparse._storage;
+ if (ReferenceEquals(this, resultSparse))
+ {
+ int i = 0, j = 0;
+ while (j < otherStorage.ValueCount)
+ {
+ if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
+ {
+ var otherValue = otherStorage.Values[j];
+ if (!Complex.Zero.Equals(otherValue))
+ {
+ _storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], -otherValue);
+ }
+
+ j++;
+ }
+ else if (_storage.Indices[i] == otherStorage.Indices[j])
+ {
+ // TODO: result can be zero, remove?
+ _storage.Values[i++] -= otherStorage.Values[j++];
+ }
+ else
+ {
+ i++;
+ }
+ }
+ }
+ else
+ {
+ result.Clear();
+ int i = 0, j = 0, last = -1;
+ while (i < _storage.ValueCount || j < otherStorage.ValueCount)
+ {
+ if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
+ {
+ var next = _storage.Indices[i];
+ if (next != last)
+ {
+ last = next;
+ result.At(next, _storage.Values[i] - otherSparse.At(next));
+ }
+
+ i++;
+ }
+ else
+ {
+ var next = otherStorage.Indices[j];
+ if (next != last)
+ {
+ last = next;
+ result.At(next, At(next) - otherStorage.Values[j]);
+ }
+
+ j++;
+ }
+ }
+ }
+ }
+ else
+ {
+ base.DoSubtract(other, result);
+ }
+ }
+
+ ///
+ /// Negates vector and saves result to
+ ///
+ /// Target vector
+ protected override void DoNegate(Vector result)
+ {
+ if (result is SparseVector sparseResult)
+ {
+ if (!ReferenceEquals(this, result))
+ {
+ sparseResult._storage.ValueCount = _storage.ValueCount;
+ sparseResult._storage.Indices = new int[_storage.ValueCount];
+ Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt);
+ sparseResult._storage.Values = new Complex[_storage.ValueCount];
+ Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
+ }
+
+ LinearAlgebraControl.Provider.ScaleArray(-Complex.One, sparseResult._storage.Values, sparseResult._storage.Values);
+ }
+ else
+ {
+ result.Clear();
+ for (var index = 0; index < _storage.ValueCount; index++)
+ {
+ result.At(_storage.Indices[index], -_storage.Values[index]);
+ }
+ }
+ }
+
+ ///
+ /// Conjugates vector and save result to
+ ///
+ /// Target vector
+ protected override void DoConjugate(Vector result)
+ {
+ if (result is SparseVector sparseResult)
+ {
+ if (!ReferenceEquals(this, result))
+ {
+ sparseResult._storage.ValueCount = _storage.ValueCount;
+ sparseResult._storage.Indices = new int[_storage.ValueCount];
+ Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount*Constants.SizeOfInt);
+ sparseResult._storage.Values = new Complex[_storage.ValueCount];
+ Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
+ }
+
+ LinearAlgebraControl.Provider.ConjugateArray(sparseResult._storage.Values, sparseResult._storage.Values);
+ return;
+ }
+
+ result.Clear();
+ for (var index = 0; index < _storage.ValueCount; index++)
+ {
+ result.At(_storage.Indices[index], _storage.Values[index].Conjugate());
+ }
+ }
+
+ ///
+ /// Multiplies a scalar to each element of the vector and stores the result in the result vector.
+ ///
+ ///
+ /// The scalar to multiply.
+ ///
+ ///
+ /// The vector to store the result of the multiplication.
+ ///
+ protected override void DoMultiply(Complex scalar, Vector result)
+ {
+ if (result is SparseVector sparseResult)
+ {
+ if (!ReferenceEquals(this, result))
+ {
+ sparseResult._storage.ValueCount = _storage.ValueCount;
+ sparseResult._storage.Indices = new int[_storage.ValueCount];
+ Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt);
+ sparseResult._storage.Values = new Complex[_storage.ValueCount];
+ Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
+ }
+
+ LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
+ }
+ else
+ {
+ result.Clear();
+ for (var index = 0; index < _storage.ValueCount; index++)
+ {
+ result.At(_storage.Indices[index], scalar * _storage.Values[index]);
+ }
+ }
+ }
+
+ ///
+ /// Computes the dot product between this vector and another vector.
+ ///
+ /// The other vector.
+ /// The sum of a[i]*b[i] for all i.
+ protected override Complex DoDotProduct(Vector other)
+ {
+ var result = Complex.Zero;
+ if (ReferenceEquals(this, other))
+ {
+ for (var i = 0; i < _storage.ValueCount; i++)
+ {
+ result += _storage.Values[i] * _storage.Values[i];
+ }
+ }
+ else
+ {
+ for (var i = 0; i < _storage.ValueCount; i++)
+ {
+ result += _storage.Values[i] * other.At(_storage.Indices[i]);
+ }
+ }
+ return result;
+ }
+
+ ///
+ /// Computes the dot product between the conjugate of this vector and another vector.
+ ///
+ /// The other vector.
+ /// The sum of conj(a[i])*b[i] for all i.
+ protected override Complex DoConjugateDotProduct(Vector other)
+ {
+ var result = Complex.Zero;
+ if (ReferenceEquals(this, other))
+ {
+ for (var i = 0; i < _storage.ValueCount; i++)
+ {
+ result += _storage.Values[i].Conjugate() * _storage.Values[i];
+ }
+ }
+ else
+ {
+ for (var i = 0; i < _storage.ValueCount; i++)
+ {
+ result += _storage.Values[i].Conjugate() * other.At(_storage.Indices[i]);
+ }
+ }
+ return result;
+ }
+
+ ///
+ /// Adds two Vectors together and returns the results.
+ ///
+ /// One of the vectors to add.
+ /// The other vector to add.
+ /// The result of the addition.
+ /// If and are not the same size.
+ /// If or is .
+ public static SparseVector operator +(SparseVector leftSide, SparseVector rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ return (SparseVector)leftSide.Add(rightSide);
+ }
+
+ ///
+ /// Returns a Vector containing the negated values of .
+ ///
+ /// The vector to get the values from.
+ /// A vector containing the negated values as .
+ /// If is .
+ public static SparseVector operator -(SparseVector rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException(nameof(rightSide));
+ }
+
+ return (SparseVector)rightSide.Negate();
+ }
+
+ ///
+ /// Subtracts two Vectors and returns the results.
+ ///
+ /// The vector to subtract from.
+ /// The vector to subtract.
+ /// The result of the subtraction.
+ /// If and are not the same size.
+ /// If or is .
+ public static SparseVector operator -(SparseVector leftSide, SparseVector rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ return (SparseVector)leftSide.Subtract(rightSide);
+ }
+
+ ///
+ /// Multiplies a vector with a complex.
+ ///
+ /// The vector to scale.
+ /// The complex value.
+ /// The result of the multiplication.
+ /// If is .
+ public static SparseVector operator *(SparseVector leftSide, Complex rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ return (SparseVector)leftSide.Multiply(rightSide);
+ }
+
+ ///
+ /// Multiplies a vector with a complex.
+ ///
+ /// The complex value.
+ /// The vector to scale.
+ /// The result of the multiplication.
+ /// If is .
+ public static SparseVector operator *(Complex leftSide, SparseVector rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException(nameof(rightSide));
+ }
+
+ return (SparseVector)rightSide.Multiply(leftSide);
+ }
+
+ ///
+ /// Computes the dot product between two Vectors.
+ ///
+ /// The left row vector.
+ /// The right column vector.
+ /// The dot product between the two vectors.
+ /// If and are not the same size.
+ /// If or is .
+ public static Complex operator *(SparseVector leftSide, SparseVector rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ return leftSide.DotProduct(rightSide);
+ }
+
+ ///
+ /// Divides a vector with a complex.
+ ///
+ /// The vector to divide.
+ /// The complex value.
+ /// The result of the division.
+ /// If is .
+ public static SparseVector operator /(SparseVector leftSide, Complex rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ return (SparseVector)leftSide.Divide(rightSide);
+ }
+
+ ///
+ /// Computes the modulus of each element of the vector of the given divisor.
+ ///
+ /// The vector whose elements we want to compute the modulus of.
+ /// The divisor to use,
+ /// The result of the calculation
+ /// If is .
+ public static SparseVector operator %(SparseVector leftSide, Complex rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException(nameof(leftSide));
+ }
+
+ return (SparseVector)leftSide.Modulus(rightSide);
+ }
+
+ ///
+ /// Returns the index of the absolute minimum element.
+ ///
+ /// The index of absolute minimum element.
+ public override int AbsoluteMinimumIndex()
+ {
+ if (_storage.ValueCount == 0)
+ {
+ // No non-zero elements. Return 0
+ return 0;
+ }
+
+ var index = 0;
+ var min = _storage.Values[index].Magnitude;
+ for (var i = 1; i < _storage.ValueCount; i++)
+ {
+ var test = _storage.Values[i].Magnitude;
+ if (test < min)
+ {
+ index = i;
+ min = test;
+ }
+ }
+
+ return _storage.Indices[index];
+ }
+
+ ///
+ /// Returns the index of the absolute maximum element.
+ ///
+ /// The index of absolute maximum element.
+ public override int AbsoluteMaximumIndex()
+ {
+ if (_storage.ValueCount == 0)
+ {
+ // No non-zero elements. Return 0
+ return 0;
+ }
+
+ var index = 0;
+ var max = _storage.Values[index].Magnitude;
+ for (var i = 1; i < _storage.ValueCount; i++)
+ {
+ var test = _storage.Values[i].Magnitude;
+ if (test > max)
+ {
+ index = i;
+ max = test;
+ }
+ }
+
+ return _storage.Indices[index];
+ }
+
+ ///
+ /// Computes the sum of the vector's elements.
+ ///
+ /// The sum of the vector's elements.
+ public override Complex Sum()
+ {
+ var result = Complex.Zero;
+ for (var i = 0; i < _storage.ValueCount; i++)
+ {
+ result += _storage.Values[i];
+ }
+ return result;
+ }
+
+ ///
+ /// Calculates the L1 norm of the vector, also known as Manhattan norm.
+ ///
+ /// The sum of the absolute values.
+ public override double L1Norm()
+ {
+ double result = 0d;
+ for (var i = 0; i < _storage.ValueCount; i++)
+ {
+ result += _storage.Values[i].Magnitude;
+ }
+ return result;
+ }
+
+ ///
+ /// Calculates the infinity norm of the vector.
+ ///
+ /// The maximum absolute value.
+ public override double InfinityNorm()
+ {
+ return CommonParallel.Aggregate(0, _storage.ValueCount, i => _storage.Values[i].Magnitude, Math.Max, 0d);
+ }
+
+ ///
+ /// Computes the p-Norm.
+ ///
+ /// The p value.
+ /// Scalar ret = ( ∑|this[i]|^p )^(1/p)
+ public override double Norm(double p)
+ {
+ if (p < 0d) throw new ArgumentOutOfRangeException(nameof(p));
+
+ if (_storage.ValueCount == 0)
+ {
+ return 0d;
+ }
+
+ if (p == 1d) return L1Norm();
+ if (p == 2d) return L2Norm();
+ if (double.IsPositiveInfinity(p)) return InfinityNorm();
+
+ double sum = 0d;
+ for (var index = 0; index < _storage.ValueCount; index++)
+ {
+ sum += Math.Pow(_storage.Values[index].Magnitude, p);
+ }
+ return Math.Pow(sum, 1.0 / p);
+ }
+
+ ///
+ /// Pointwise multiplies this vector with another vector and stores the result into the result vector.
+ ///
+ /// The vector to pointwise multiply with this one.
+ /// The vector to store the result of the pointwise multiplication.
+ protected override void DoPointwiseMultiply(Vector other, Vector result)
+ {
+ if (ReferenceEquals(this, other) && ReferenceEquals(this, result))
+ {
+ for (var i = 0; i < _storage.ValueCount; i++)
+ {
+ _storage.Values[i] *= _storage.Values[i];
+ }
+ }
+ else
+ {
+ base.DoPointwiseMultiply(other, result);
+ }
+ }
+
+ #region Parse Functions
+
+ ///
+ /// Creates a double sparse vector based on a string. The string can be in the following formats (without the
+ /// quotes): 'n', 'n;n;..', '(n;n;..)', '[n;n;...]', where n is a Complex.
+ ///
+ ///
+ /// A double sparse vector containing the values specified by the given string.
+ ///
+ ///
+ /// the string to parse.
+ ///
+ ///
+ /// An that supplies culture-specific formatting information.
+ ///
+ public static SparseVector Parse(string value, IFormatProvider formatProvider = null)
+ {
+ if (value == null)
+ {
+ throw new ArgumentNullException(nameof(value));
+ }
+
+ value = value.Trim();
+ if (value.Length == 0)
+ {
+ throw new FormatException();
+ }
+
+ // strip out parens
+ if (value.StartsWith("(", StringComparison.Ordinal))
+ {
+ if (!value.EndsWith(")", StringComparison.Ordinal))
+ {
+ throw new FormatException();
+ }
+
+ value = value.Substring(1, value.Length - 2).Trim();
+ }
+
+ if (value.StartsWith("[", StringComparison.Ordinal))
+ {
+ if (!value.EndsWith("]", StringComparison.Ordinal))
+ {
+ throw new FormatException();
+ }
+
+ value = value.Substring(1, value.Length - 2).Trim();
+ }
+
+ // parsing
+ var strongTokens = value.Split(new[] { formatProvider.GetTextInfo().ListSeparator }, StringSplitOptions.RemoveEmptyEntries);
+ var data = new List();
+ foreach (string strongToken in strongTokens)
+ {
+ var weakTokens = strongToken.Split(new[] { " ", "\t" }, StringSplitOptions.RemoveEmptyEntries);
+ string current = string.Empty;
+ for (int i = 0; i < weakTokens.Length; i++)
+ {
+ current += weakTokens[i];
+ if (current.EndsWith("+") || current.EndsWith("-") || current.StartsWith("(") && !current.EndsWith(")"))
+ {
+ continue;
+ }
+ var ahead = i < weakTokens.Length - 1 ? weakTokens[i + 1] : string.Empty;
+ if (ahead.StartsWith("+") || ahead.StartsWith("-"))
+ {
+ continue;
+ }
+ data.Add(current.ToComplex(formatProvider));
+ current = string.Empty;
+ }
+ if (current != string.Empty)
+ {
+ throw new FormatException();
+ }
+ }
+ if (data.Count == 0)
+ {
+ throw new FormatException();
+ }
+ return OfEnumerable(data);
+ }
+
+ ///
+ /// Converts the string representation of a complex sparse vector to double-precision sparse vector equivalent.
+ /// A return value indicates whether the conversion succeeded or failed.
+ ///
+ ///
+ /// A string containing a complex vector to convert.
+ ///
+ ///
+ /// The parsed value.
+ ///
+ ///
+ /// If the conversion succeeds, the result will contain a complex number equivalent to value.
+ /// Otherwise the result will be null.
+ ///
+ public static bool TryParse(string value, out SparseVector result)
+ {
+ return TryParse(value, null, out result);
+ }
+
+ ///
+ /// Converts the string representation of a complex sparse vector to double-precision sparse vector equivalent.
+ /// A return value indicates whether the conversion succeeded or failed.
+ ///
+ ///
+ /// A string containing a complex vector to convert.
+ ///
+ ///
+ /// An that supplies culture-specific formatting information about value.
+ ///
+ ///
+ /// The parsed value.
+ ///
+ ///
+ /// If the conversion succeeds, the result will contain a complex number equivalent to value.
+ /// Otherwise the result will be null.
+ ///
+ public static bool TryParse(string value, IFormatProvider formatProvider, out SparseVector result)
+ {
+ bool ret;
+ try
+ {
+ result = Parse(value, formatProvider);
+ ret = true;
+ }
+ catch (ArgumentNullException)
+ {
+ result = null;
+ ret = false;
+ }
+ catch (FormatException)
+ {
+ result = null;
+ ret = false;
+ }
+
+ return ret;
+ }
+ #endregion
+
+ public override string ToTypeString()
+ {
+ return FormattableString.Invariant($"SparseVector {Count}-Complex {NonZerosCount / (double) Count:P2} Filled");
+ }
+ }
+}
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseVector.cs.meta b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseVector.cs.meta
new file mode 100644
index 00000000..9b0468c4
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/SparseVector.cs.meta
@@ -0,0 +1,11 @@
+fileFormatVersion: 2
+guid: 35f129457993fad4ea51c9731320b863
+MonoImporter:
+ externalObjects: {}
+ serializedVersion: 2
+ defaultReferences: []
+ executionOrder: 0
+ icon: {instanceID: 0}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Vector.cs b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Vector.cs
new file mode 100644
index 00000000..018ddb1c
--- /dev/null
+++ b/Assets/Plugins/UnityMeshDecimaton/3rdParty/Math.Net/Numerics/LinearAlgebra/Complex/Vector.cs
@@ -0,0 +1,618 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2015 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.Storage;
+using MathNet.Numerics.Threading;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex
+{
+ using Complex = System.Numerics.Complex;
+
+ ///
+ /// Complex version of the class.
+ ///
+ [Serializable]
+ public abstract class Vector : Vector
+ {
+ ///
+ /// Initializes a new instance of the Vector class.
+ ///
+ protected Vector(VectorStorage storage)
+ : base(storage)
+ {
+ }
+
+ ///
+ /// Set all values whose absolute value is smaller than the threshold to zero.
+ ///
+ public override void CoerceZero(double threshold)
+ {
+ MapInplace(x => x.Magnitude < threshold ? Complex.Zero : x, Zeros.AllowSkip);
+ }
+
+ ///
+ /// Conjugates vector and save result to
+ ///
+ /// Target vector
+ protected override void DoConjugate(Vector result)
+ {
+ Map(Complex.Conjugate, result, Zeros.AllowSkip);
+ }
+
+ ///
+ /// Negates vector and saves result to
+ ///
+ /// Target vector
+ protected override void DoNegate(Vector result)
+ {
+ Map(Complex.Negate, result, Zeros.AllowSkip);
+ }
+
+ ///
+ /// Adds a scalar to each element of the vector and stores the result in the result vector.
+ ///
+ ///
+ /// The scalar to add.
+ ///
+ ///
+ /// The vector to store the result of the addition.
+ ///
+ protected override void DoAdd(Complex scalar, Vector result)
+ {
+ Map(x => x + scalar, result, Zeros.Include);
+ }
+
+ ///
+ /// Adds another vector to this vector and stores the result into the result vector.
+ ///
+ ///
+ /// The vector to add to this one.
+ ///
+ ///
+ /// The vector to store the result of the addition.
+ ///
+ protected override void DoAdd(Vector other, Vector result)
+ {
+ Map2(Complex.Add, other, result, Zeros.AllowSkip);
+ }
+
+ ///
+ /// Subtracts a scalar from each element of the vector and stores the result in the result vector.
+ ///
+ ///
+ /// The scalar to subtract.
+ ///
+ ///
+ /// The vector to store the result of the subtraction.
+ ///
+ protected override void DoSubtract(Complex scalar, Vector result)
+ {
+ Map(x => x - scalar, result, Zeros.Include);
+ }
+
+ ///
+ /// Subtracts another vector to this vector and stores the result into the result vector.
+ ///
+ ///
+ /// The vector to subtract from this one.
+ ///
+ ///
+ /// The vector to store the result of the subtraction.
+ ///
+ protected override void DoSubtract(Vector other, Vector result)
+ {
+ Map2(Complex.Subtract, other, result, Zeros.AllowSkip);
+ }
+
+ ///
+ /// Multiplies a scalar to each element of the vector and stores the result in the result vector.
+ ///
+ ///
+ /// The scalar to multiply.
+ ///
+ ///
+ /// The vector to store the result of the multiplication.
+ ///
+ protected override void DoMultiply(Complex scalar, Vector result)
+ {
+ Map(x => x*scalar, result, Zeros.AllowSkip);
+ }
+
+ ///
+ /// Divides each element of the vector by a scalar and stores the result in the result vector.
+ ///
+ ///
+ /// The scalar to divide with.
+ ///
+ ///
+ /// The vector to store the result of the division.
+ ///
+ protected override void DoDivide(Complex divisor, Vector result)
+ {
+ Map(x => x/divisor, result, divisor.IsZero() ? Zeros.Include : Zeros.AllowSkip);
+ }
+
+ ///
+ /// Divides a scalar by each element of the vector and stores the result in the result vector.
+ ///
+ /// The scalar to divide.
+ /// The vector to store the result of the division.
+ protected override void DoDivideByThis(Complex dividend, Vector result)
+ {
+ Map(x => dividend/x, result, Zeros.Include);
+ }
+
+ ///
+ /// Pointwise multiplies this vector with another vector and stores the result into the result vector.
+ ///
+ /// The vector to pointwise multiply with this one.
+ /// The vector to store the result of the pointwise multiplication.
+ protected override void DoPointwiseMultiply(Vector other, Vector result)
+ {
+ Map2(Complex.Multiply, other, result, Zeros.AllowSkip);
+ }
+
+ ///
+ /// Pointwise divide this vector with another vector and stores the result into the result vector.
+ ///
+ /// The vector to pointwise divide this one by.
+ /// The vector to store the result of the pointwise division.
+ protected override void DoPointwiseDivide(Vector divisor, Vector result)
+ {
+ Map2(Complex.Divide, divisor, result, Zeros.Include);
+ }
+
+ ///
+ /// Pointwise raise this vector to an exponent and store the result into the result vector.
+ ///
+ /// The exponent to raise this vector values to.
+ /// The vector to store the result of the pointwise power.
+ protected override void DoPointwisePower(Complex exponent, Vector result)
+ {
+ Map(x => x.Power(exponent), result, Zeros.Include);
+ }
+
+ ///
+ /// Pointwise raise this vector to an exponent vector and store the result into the result vector.
+ ///
+ /// The exponent vector to raise this vector values to.
+ /// The vector to store the result of the pointwise power.
+ protected override void DoPointwisePower(Vector exponent, Vector result)
+ {
+ Map2(Complex.Pow, exponent, result, Zeros.Include);
+ }
+
+ ///
+ /// Pointwise canonical modulus, where the result has the sign of the divisor,
+ /// of this vector with another vector and stores the result into the result vector.
+ ///
+ /// The pointwise denominator vector to use.
+ /// The result of the modulus.
+ protected sealed override void DoPointwiseModulus(Vector divisor, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ ///
+ /// Pointwise remainder (% operator), where the result has the sign of the dividend,
+ /// of this vector with another vector and stores the result into the result vector.
+ ///
+ /// The pointwise denominator vector to use.
+ /// The result of the modulus.
+ protected sealed override void DoPointwiseRemainder(Vector divisor, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ ///
+ /// Pointwise applies the exponential function to each value and stores the result into the result vector.
+ ///
+ /// The vector to store the result.
+ protected override void DoPointwiseExp(Vector result)
+ {
+ Map(Complex.Exp, result, Zeros.Include);
+ }
+
+ ///
+ /// Pointwise applies the natural logarithm function to each value and stores the result into the result vector.
+ ///
+ /// The vector to store the result.
+ protected override void DoPointwiseLog(Vector result)
+ {
+ Map(Complex.Log, result, Zeros.Include);
+ }
+
+ protected override void DoPointwiseAbs(Vector result)
+ {
+ Map(x => (Complex)Complex.Abs(x), result, Zeros.AllowSkip);
+ }
+ protected override void DoPointwiseAcos(Vector result)
+ {
+ Map(Complex.Acos, result, Zeros.Include);
+ }
+ protected override void DoPointwiseAsin(Vector result)
+ {
+ Map(Complex.Asin, result, Zeros.AllowSkip);
+ }
+ protected override void DoPointwiseAtan(Vector result)
+ {
+ Map(Complex.Atan, result, Zeros.AllowSkip);
+ }
+ protected override void DoPointwiseAtan2(Vector other, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+ protected override void DoPointwiseAtan2(Complex scalar, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+ protected override void DoPointwiseCeiling(Vector result)
+ {
+ throw new NotSupportedException();
+ }
+ protected override void DoPointwiseCos(Vector result)
+ {
+ Map(Complex.Cos, result, Zeros.Include);
+ }
+ protected override void DoPointwiseCosh(Vector result)
+ {
+ Map(Complex.Cosh, result, Zeros.Include);
+ }
+ protected override void DoPointwiseFloor(Vector result)
+ {
+ throw new NotSupportedException();
+ }
+ protected override void DoPointwiseLog10(Vector result)
+ {
+ Map(Complex.Log10, result, Zeros.Include);
+ }
+ protected override void DoPointwiseRound(Vector result)
+ {
+ throw new NotSupportedException();
+ }
+ protected override void DoPointwiseSign(Vector result)
+ {
+ throw new NotSupportedException();
+ }
+ protected override void DoPointwiseSin(Vector result)
+ {
+ Map(Complex.Sin, result, Zeros.AllowSkip);
+ }
+ protected override void DoPointwiseSinh(Vector result)
+ {
+ Map(Complex.Sinh, result, Zeros.AllowSkip);
+ }
+ protected override void DoPointwiseSqrt(Vector result)
+ {
+ Map(Complex.Sqrt, result, Zeros.AllowSkip);
+ }
+ protected override void DoPointwiseTan(Vector result)
+ {
+ Map(Complex.Tan, result, Zeros.AllowSkip);
+ }
+ protected override void DoPointwiseTanh(Vector result)
+ {
+ Map(Complex.Tanh, result, Zeros.AllowSkip);
+ }
+
+ ///
+ /// Computes the dot product between this vector and another vector.
+ ///
+ /// The other vector.
+ /// The sum of a[i]*b[i] for all i.
+ protected override Complex DoDotProduct(Vector other)
+ {
+ var dot = Complex.Zero;
+ for (var i = 0; i < Count; i++)
+ {
+ dot += At(i) * other.At(i);
+ }
+ return dot;
+ }
+
+ ///
+ /// Computes the dot product between the conjugate of this vector and another vector.
+ ///
+ /// The other vector.
+ /// The sum of conj(a[i])*b[i] for all i.
+ protected override Complex DoConjugateDotProduct(Vector other)
+ {
+ var dot = Complex.Zero;
+ for (var i = 0; i < Count; i++)
+ {
+ dot += At(i).Conjugate() * other.At(i);
+ }
+ return dot;
+ }
+
+ ///
+ /// Computes the canonical modulus, where the result has the sign of the divisor,
+ /// for each element of the vector for the given divisor.
+ ///
+ /// The scalar denominator to use.
+ /// A vector to store the results in.
+ protected sealed override void DoModulus(Complex divisor, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ ///
+ /// Computes the canonical modulus, where the result has the sign of the divisor,
+ /// for the given dividend for each element of the vector.
+ ///
+ /// The scalar numerator to use.
+ /// A vector to store the results in.
+ protected sealed override void DoModulusByThis(Complex dividend, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ ///
+ /// Computes the canonical modulus, where the result has the sign of the divisor,
+ /// for each element of the vector for the given divisor.
+ ///
+ /// The scalar denominator to use.
+ /// A vector to store the results in.
+ protected sealed override void DoRemainder(Complex divisor, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ ///
+ /// Computes the canonical modulus, where the result has the sign of the divisor,
+ /// for the given dividend for each element of the vector.
+ ///
+ /// The scalar numerator to use.
+ /// A vector to store the results in.
+ protected sealed override void DoRemainderByThis(Complex dividend, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ protected override void DoPointwiseMinimum(Complex scalar, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ protected override void DoPointwiseMaximum(Complex scalar, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ protected override void DoPointwiseAbsoluteMinimum(Complex scalar, Vector result)
+ {
+ double absolute = scalar.Magnitude;
+ Map(x => Math.Min(absolute, x.Magnitude), result, Zeros.AllowSkip);
+ }
+
+ protected override void DoPointwiseAbsoluteMaximum(Complex scalar, Vector result)
+ {
+ double absolute = scalar.Magnitude;
+ Map(x => Math.Max(absolute, x.Magnitude), result, Zeros.Include);
+ }
+
+ protected override void DoPointwiseMinimum(Vector other, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ protected override void DoPointwiseMaximum(Vector other, Vector result)
+ {
+ throw new NotSupportedException();
+ }
+
+ protected override void DoPointwiseAbsoluteMinimum(Vector other, Vector result)
+ {
+ Map2((x, y) => Math.Min(x.Magnitude, y.Magnitude), other, result, Zeros.AllowSkip);
+ }
+
+ protected override void DoPointwiseAbsoluteMaximum(Vector other, Vector result)
+ {
+ Map2((x, y) => Math.Max(x.Magnitude, y.Magnitude), other, result, Zeros.AllowSkip);
+ }
+
+ ///
+ /// Returns the value of the absolute minimum element.
+ ///
+ /// The value of the absolute minimum element.
+ public sealed override Complex AbsoluteMinimum()
+ {
+ return At(AbsoluteMinimumIndex()).Magnitude;
+ }
+
+ ///
+ /// Returns the index of the absolute minimum element.
+ ///
+ /// The index of absolute minimum element.
+ public override int AbsoluteMinimumIndex()
+ {
+ var index = 0;
+ var min = At(index).Magnitude;
+ for (var i = 1; i < Count; i++)
+ {
+ var test = At(i).Magnitude;
+ if (test < min)
+ {
+ index = i;
+ min = test;
+ }
+ }
+
+ return index;
+ }
+
+ ///
+ /// Returns the value of the absolute maximum element.
+ ///
+ /// The value of the absolute maximum element.
+ public override Complex AbsoluteMaximum()
+ {
+ return At(AbsoluteMaximumIndex()).Magnitude;
+ }
+
+ ///
+ /// Returns the index of the absolute maximum element.
+ ///
+ /// The index of absolute maximum element.
+ public override int AbsoluteMaximumIndex()
+ {
+ var index = 0;
+ var max = At(index).Magnitude;
+ for (var i = 1; i < Count; i++)
+ {
+ var test = At(i).Magnitude;
+ if (test > max)
+ {
+ index = i;
+ max = test;
+ }
+ }
+
+ return index;
+ }
+
+ ///
+ /// Computes the sum of the vector's elements.
+ ///
+ /// The sum of the vector's elements.
+ public override Complex Sum()
+ {
+ var sum = Complex.Zero;
+ for (var i = 0; i < Count; i++)
+ {
+ sum += At(i);
+ }
+ return sum;
+ }
+
+ ///
+ /// Calculates the L1 norm of the vector, also known as Manhattan norm.
+ ///
+ /// The sum of the absolute values.
+ public override double L1Norm()
+ {
+ double sum = 0d;
+ for (var i = 0; i < Count; i++)
+ {
+ sum += At(i).Magnitude;
+ }
+ return sum;
+ }
+
+ ///
+ /// Calculates the L2 norm of the vector, also known as Euclidean norm.
+ ///
+ /// The square root of the sum of the squared values.
+ public override double L2Norm()
+ {
+ return DoConjugateDotProduct(this).SquareRoot().Real;
+ }
+
+ ///
+ /// Calculates the infinity norm of the vector.
+ ///
+ /// The maximum absolute value.
+ public override double InfinityNorm()
+ {
+ return CommonParallel.Aggregate(0, Count, i => At(i).Magnitude, Math.Max, 0d);
+ }
+
+ ///
+ /// Computes the p-Norm.
+ ///
+ ///
+ /// The p value.
+ ///
+ ///
+ ///