tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

native: added LU solve 

api: removed the LU solve methods with the transpose option from LA provider interface for consistency with other solvers

added SecuritySafeCritical and SecurityCritical attributes native provider and p/invoke methods
pull/36/head
Marcus Cuda 16 years ago
parent
4628fc4de6
commit
60f4ee5cc5
16 changed files with 739 additions and 1037 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 4
      src/MathNet.Numerics.5.1.ReSharper
  2. 1091
      src/NativeWrappers/MKL/lapack.cpp
  3. 25
      src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
  4. 120
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
  5. 120
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
  6. 110
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
  7. 107
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
  8. 8
      src/Numerics/Algorithms/LinearAlgebra/native.common.include
  9. 93
      src/Numerics/Algorithms/LinearAlgebra/native.generic.include
  10. 1
      src/Numerics/Algorithms/LinearAlgebra/native.header.include
  11. 19
      src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include
  12. 3
      src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs
  13. 3
      src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs
  14. 1
      src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs
  15. 3
      src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs
  16. 68
      src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs

4
src/MathNet.Numerics.5.1.ReSharper
View File

@ -83,7 +83,9 @@ kronecker
blocksize
ipiv
dll
Silverlight</UserWords>
Silverlight
namespace
da</UserWords>
</CustomDictionary>
</Dictionaries>
</CustomDictionaries>

1091
src/NativeWrappers/MKL/lapack.cpp
View File

File diff suppressed because it is too large

25
src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
View File

@ -3,7 +3,7 @@
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
// Copyright (c) 2009-2010 Math.NET
// Copyright (c) 2009-2011 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
@ -283,29 +283,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
void LUSolveFactored(int columnsOfB, T[] a, int order, int[] ipiv, T[] b);
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
void LUSolve(Transpose transposeA, int columnsOfB, T[] a, int order, T[] b);
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
void LUSolveFactored(Transpose transposeA, int columnsOfB, T[] a, int order, int[] ipiv, T[] b);
/// <summary>
/// Computes the Cholesky factorization of A.
/// </summary>

120
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
View File

@ -1026,7 +1026,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
public virtual void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b)
{
@ -1051,8 +1051,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
var ipiv = new int[order];
LUFactor(a, order, ipiv);
LUSolveFactored(columnsOfB, a, order, ipiv, b);
var clone = new Complex[a.Length];
Array.Copy(a, 0, clone, 0, a.Length);
LUFactor(clone, order, ipiv);
LUSolveFactored(columnsOfB, clone, order, ipiv, b);
}
/// <summary>
@ -1062,7 +1064,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
public virtual void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b)
{
@ -1151,116 +1153,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
public virtual void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, int order, Complex[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
var ipiv = new int[order];
LUFactor(a, order, ipiv);
LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b);
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
public virtual void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if (transposeA == Transpose.Transpose)
{
var aT = new Complex[a.Length];
for (var i = 0; i < order; i++)
{
for (var j = 0; j < order; j++)
{
aT[(j * order) + i] = a[(i * order) + j];
}
}
LUSolveFactored(columnsOfB, aT, order, ipiv, b);
}
else if (transposeA == Transpose.ConjugateTranspose)
{
var acT = new Complex[a.Length];
for (var i = 0; i < order; i++)
{
for (var j = 0; j < order; j++)
{
acT[(j * order) + i] = a[(i * order) + j].Conjugate();
}
}
LUSolveFactored(columnsOfB, acT, order, ipiv, b);
}
else
{
LUSolveFactored(columnsOfB, a, order, ipiv, b);
}
}
/// <summary>
/// Computes the Cholesky factorization of A.
/// </summary>

120
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
View File

@ -1026,7 +1026,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
public virtual void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b)
{
@ -1051,8 +1051,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
var ipiv = new int[order];
LUFactor(a, order, ipiv);
LUSolveFactored(columnsOfB, a, order, ipiv, b);
var clone = new Complex32[a.Length];
Array.Copy(a, 0, clone, 0, a.Length);
LUFactor(clone, order, ipiv);
LUSolveFactored(columnsOfB, clone, order, ipiv, b);
}
/// <summary>
@ -1062,7 +1064,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
public virtual void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b)
{
@ -1151,116 +1153,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
public virtual void LUSolve(Transpose transposeA, int columnsOfB, Complex32[] a, int order, Complex32[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
var ipiv = new int[order];
LUFactor(a, order, ipiv);
LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b);
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
public virtual void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if (transposeA == Transpose.Transpose)
{
var aT = new Complex32[a.Length];
for (var i = 0; i < order; i++)
{
for (var j = 0; j < order; j++)
{
aT[(j * order) + i] = a[(i * order) + j];
}
}
LUSolveFactored(columnsOfB, aT, order, ipiv, b);
}
else if (transposeA == Transpose.ConjugateTranspose)
{
var acT = new Complex32[a.Length];
for (var i = 0; i < order; i++)
{
for (var j = 0; j < order; j++)
{
acT[(j * order) + i] = a[(i * order) + j].Conjugate();
}
}
LUSolveFactored(columnsOfB, acT, order, ipiv, b);
}
else
{
LUSolveFactored(columnsOfB, a, order, ipiv, b);
}
}
/// <summary>
/// Computes the Cholesky factorization of A.
/// </summary>

110
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
View File

@ -1020,7 +1020,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
public virtual void LUSolve(int columnsOfB, double[] a, int order, double[] b)
{
@ -1045,8 +1045,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
var ipiv = new int[order];
LUFactor(a, order, ipiv);
LUSolveFactored(columnsOfB, a, order, ipiv, b);
var clone = new double[a.Length];
Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble);
LUFactor(clone, order, ipiv);
LUSolveFactored(columnsOfB, clone, order, ipiv, b);
}
/// <summary>
@ -1056,7 +1058,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
public virtual void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b)
{
@ -1145,103 +1147,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
public virtual void LUSolve(Transpose transposeA, int columnsOfB, double[] a, int order, double[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
var ipiv = new int[order];
LUFactor(a, order, ipiv);
LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b);
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
public virtual void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int order, int[] ipiv, double[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if ((transposeA == Transpose.Transpose) || (transposeA == Transpose.ConjugateTranspose))
{
var aT = new double[a.Length];
for (var i = 0; i < order; i++)
{
for (var j = 0; j < order; j++)
{
aT[(j * order) + i] = a[(i * order) + j];
}
}
LUSolveFactored(columnsOfB, aT, order, ipiv, b);
}
else
{
LUSolveFactored(columnsOfB, a, order, ipiv, b);
}
}
/// <summary>
/// Computes the Cholesky factorization of A.
/// </summary>
@ -1829,7 +1734,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
var im = i * rowsR;
var sum = 0.0;
for (var k = 0; k < rowsR; k++ )
for (var k = 0; k < rowsR; k++)
{
sum += q[im + k] * column[k];
}
@ -2021,7 +1926,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
{
// Compute the transformation for the l-th column and
// place the l-th diagonal in vector s[l].
var sum = 0.0;
for (var i1 = l; i1 < rowsA; i1++)
{

107
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
View File

@ -1021,7 +1021,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
public virtual void LUSolve(int columnsOfB, float[] a, int order, float[] b)
{
@ -1046,8 +1046,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
var ipiv = new int[order];
LUFactor(a, order, ipiv);
LUSolveFactored(columnsOfB, a, order, ipiv, b);
var clone = new float[a.Length];
Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat);
LUFactor(clone, order, ipiv);
LUSolveFactored(columnsOfB, clone, order, ipiv, b);
}
/// <summary>
@ -1057,7 +1059,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
public virtual void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b)
{
@ -1146,103 +1148,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
public virtual void LUSolve(Transpose transposeA, int columnsOfB, float[] a, int order, float[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
var ipiv = new int[order];
LUFactor(a, order, ipiv);
LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b);
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
public virtual void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int order, int[] ipiv, float[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if ((transposeA == Transpose.Transpose) || (transposeA == Transpose.ConjugateTranspose))
{
var aT = new float[a.Length];
for (var i = 0; i < order; i++)
{
for (var j = 0; j < order; j++)
{
aT[(j * order) + i] = a[(i * order) + j];
}
}
LUSolveFactored(columnsOfB, aT, order, ipiv, b);
}
else
{
LUSolveFactored(columnsOfB, a, order, ipiv, b);
}
}
/// <summary>
/// Computes the Cholesky factorization of A.
/// </summary>

8
src/Numerics/Algorithms/LinearAlgebra/native.common.include
View File

@ -8,6 +8,7 @@
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[SecuritySafeCritical]
public override float MatrixNorm(Norm norm, int rows, int columns, float[] matrix)
{
if (matrix == null)
@ -46,6 +47,7 @@
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[SecuritySafeCritical]
public override float MatrixNorm(Norm norm, int rows, int columns, float[] matrix, float[] work)
{
if (matrix == null)
@ -86,6 +88,7 @@
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[SecuritySafeCritical]
public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix)
{
if (matrix == null)
@ -124,6 +127,7 @@
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[SecuritySafeCritical]
public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix, double[] work)
{
if (matrix == null)
@ -164,6 +168,7 @@
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[SecuritySafeCritical]
public override Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix)
{
if (matrix == null)
@ -202,6 +207,7 @@
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[SecuritySafeCritical]
public override Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix, float[] work)
{
if (matrix == null)
@ -242,6 +248,7 @@
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[SecuritySafeCritical]
public override Complex MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix)
{
if (matrix == null)
@ -280,6 +287,7 @@
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[SecuritySafeCritical]
public override Complex MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix, double[] work)
{
if (matrix == null)

93
src/Numerics/Algorithms/LinearAlgebra/native.generic.include
View File

@ -6,6 +6,7 @@
/// <param name="x">The vector to add to <paramref name="y"/>.</param>
/// <param name="result">The result of the addition.</param>
/// <remarks>This is similar to the AXPY BLAS routine.</remarks>
[SecuritySafeCritical]
public override void AddVectorToScaledVector(<#=dataType#>[] y, <#=dataType#> alpha, <#=dataType#>[] x, <#=dataType#>[] result)
{
if (y == null)
@ -43,6 +44,7 @@
/// <param name="x">The values to scale.</param>
/// <param name="result">This result of the scaling.</param>
/// <remarks>This is similar to the SCAL BLAS routine.</remarks>
[SecuritySafeCritical]
public override void ScaleArray(<#=dataType#> alpha, <#=dataType#>[] x, <#=dataType#>[] result)
{
if (x == null)
@ -70,6 +72,7 @@
/// <param name="y">The vector y.</param>
/// <returns>The dot product of x and y.</returns>
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
[SecuritySafeCritical]
public override <#=dataType#> DotProduct(<#=dataType#>[] x, <#=dataType#>[] y)
{
if (y == null)
@ -121,6 +124,7 @@
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
[SecuritySafeCritical]
public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, <#=dataType#> alpha, <#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] b, int rowsB, int columnsB, <#=dataType#> beta, <#=dataType#>[] c)
{
if (a == null)
@ -164,6 +168,7 @@
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void LUFactor(<#=dataType#>[] data, int order, int[] ipiv)
{
if (data == null)
@ -195,6 +200,7 @@
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[SecuritySafeCritical]
public override void LUInverse(<#=dataType#>[] a, int order)
{
if (a == null)
@ -218,6 +224,7 @@
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void LUInverseFactored(<#=dataType#>[] a, int order, int[] ipiv)
{
if (a == null)
@ -253,6 +260,7 @@
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[SecuritySafeCritical]
public override void LUInverse(<#=dataType#>[] a, int order, <#=dataType#>[] work)
{
if (a == null)
@ -288,6 +296,7 @@
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void LUInverseFactored(<#=dataType#>[] a, int order, int[] ipiv, <#=dataType#>[] work)
{
if (a == null)
@ -331,9 +340,25 @@
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
[SecuritySafeCritical]
public override void LUSolve(int columnsOfB, <#=dataType#>[] a, int order, <#=dataType#>[] b)
{
throw new NotImplementedException();
if (a == null)
{
throw new ArgumentNullException("a");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (b.Length != columnsOfB * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
SafeNativeMethods.<#=prefix#>_lu_solve(order, columnsOfB, a, b);
}
/// <summary>
@ -345,38 +370,35 @@
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void LUSolveFactored(int columnsOfB, <#=dataType#>[] a, int order, int[] ipiv, <#=dataType#>[] b)
{
throw new NotImplementedException();
}
if (a == null)
{
throw new ArgumentNullException("a");
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
public override void LUSolve(Transpose transposeA, int columnsOfB, <#=dataType#>[] a, int order, <#=dataType#>[] b)
{
throw new NotImplementedException();
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
public override void LUSolveFactored(Transpose transposeA, int columnsOfB, <#=dataType#>[] a, int order, int[] ipiv, <#=dataType#>[] b)
{
throw new NotImplementedException();
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (b.Length != columnsOfB * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
SafeNativeMethods.<#=prefix#>_lu_solve_factored(order, columnsOfB, a, ipiv, b);
}
/// <summary>
@ -386,6 +408,7 @@
/// the Cholesky factorization.</param>
/// <param name="order">The number of rows or columns in the matrix.</param>
/// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void CholeskyFactor(<#=dataType#>[] a, int order)
{
if (a == null)
@ -416,6 +439,7 @@
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.
/// </remarks>
[SecuritySafeCritical]
public override void CholeskySolve(<#=dataType#>[] a, int orderA, <#=dataType#>[] b, int rowsB, int columnsB)
{
throw new NotImplementedException();
@ -430,6 +454,7 @@
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void CholeskySolveFactored(<#=dataType#>[] a, int orderA, <#=dataType#>[] b, int rowsB, int columnsB)
{
throw new NotImplementedException();
@ -445,6 +470,7 @@
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[SecuritySafeCritical]
public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q)
{
throw new NotImplementedException();
@ -463,6 +489,7 @@
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[SecuritySafeCritical]
public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] work)
{
throw new NotImplementedException();
@ -480,6 +507,7 @@
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
[SecuritySafeCritical]
public override void QRSolve(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x)
{
throw new NotImplementedException();
@ -500,6 +528,7 @@
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
[SecuritySafeCritical]
public override void QRSolve(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x, <#=dataType#>[] work)
{
throw new NotImplementedException();
@ -515,6 +544,7 @@
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
[SecuritySafeCritical]
public override void QRSolveFactored(<#=dataType#>[] q, <#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x)
{
throw new NotImplementedException();
@ -533,6 +563,7 @@
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, <#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt)
{
throw new NotImplementedException();
@ -554,6 +585,7 @@
/// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
/// On exit, work[0] contains the optimal work size value.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, <#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] work)
{
throw new NotImplementedException();
@ -571,6 +603,7 @@
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
[SecuritySafeCritical]
public override void SvdSolve(<#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x)
{
throw new NotImplementedException();
@ -591,6 +624,7 @@
/// <param name="work">The work array. For real matrices, the work array should be at least
/// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
/// On exit, work[0] contains the optimal work size value.</param>
[SecuritySafeCritical]
public override void SvdSolve(<#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x, <#=dataType#>[] work)
{
throw new NotImplementedException();
@ -607,6 +641,7 @@
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
[SecuritySafeCritical]
public override void SvdSolveFactored(int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x)
{
throw new NotImplementedException();

1
src/Numerics/Algorithms/LinearAlgebra/native.header.include
View File

@ -36,6 +36,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#>
{
using System;
using System.Numerics;
using System.Security;
using Properties;
/// <summary>

19
src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include
View File

@ -40,6 +40,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #>
/// P/Invoke methods to the native math libraries.
/// </summary>
[SuppressUnmanagedCodeSecurity]
[SecurityCritical]
internal static class SafeNativeMethods
{
/// <summary>
@ -161,27 +162,27 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #>
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void z_lu_inverse_factored(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] work, int lwork);
/* [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void s_lu_solve_factored(int n, int nrhs, float a[], int ipiv[], float b[])
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void s_lu_solve_factored(int n, int nrhs, float[] a, [In, Out]int[] ipiv, [In, Out] float[] b);
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void d_lu_solve_factored(int n, int nrhs, double a[], int ipiv[], double b[])
internal static extern void d_lu_solve_factored(int n, int nrhs, double[] a, [In, Out] int[] ipiv, [In, Out] double[] b);
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void c_lu_solve_factored(int n, int nrhs, MKL_Complex8 a[], int ipiv[], MKL_Complex8 b[])
internal static extern void c_lu_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] b);
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void z_lu_solve_factored(int n, int nrhs, MKL_Complex16 a[], int ipiv[], MKL_Complex16 b[])
internal static extern void z_lu_solve_factored(int n, int nrhs, Complex[] a, [In, Out]int[] ipiv, [In, Out] Complex[] b);
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void s_lu_solve(int n, int nrhs, float a[], float b[])
internal static extern void s_lu_solve(int n, int nrhs, float[] a, [In, Out] float[] b);
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void d_lu_solve(int n, int nrhs, double a[], double b[])
internal static extern void d_lu_solve(int n, int nrhs, double[] a, [In, Out] double[] b);
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void c_lu_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[])
internal static extern void c_lu_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b);
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void z_lu_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[])*/
internal static extern void z_lu_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b);
#endregion LAPACK

3
src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs
View File

@ -32,6 +32,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
{
using System;
using System.Numerics;
using Algorithms.LinearAlgebra;
using Generic;
using Properties;
using Threading;
@ -125,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
// LU solve by overwriting result.
var dfactors = (DenseMatrix)Factors;
Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data);
Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data);
}
/// <summary>

3
src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs
View File

@ -31,6 +31,7 @@
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
{
using System;
using Algorithms.LinearAlgebra;
using Generic;
using Numerics;
using Properties;
@ -125,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
// LU solve by overwriting result.
var dfactors = (DenseMatrix)Factors;
Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data);
Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data);
}
/// <summary>

1
src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs
View File

@ -31,6 +31,7 @@
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
{
using System;
using Algorithms.LinearAlgebra;
using Generic;
using Properties;

3
src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs
View File

@ -31,6 +31,7 @@
namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
{
using System;
using Algorithms.LinearAlgebra;
using Generic;
using Properties;
@ -123,7 +124,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
// LU solve by overwriting result.
var dfactors = (DenseMatrix)Factors;
Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data);
Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data);
}
/// <summary>

68
src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs
View File

@ -30,6 +30,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
using System.Collections.Generic;
using Algorithms.LinearAlgebra;
using LinearAlgebra.Double;
using MathNet.Numerics.LinearAlgebra.Generic;
using NUnit.Framework;
/// <summary>
@ -548,5 +551,70 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqual(a[8], -0.113636363636364, 14);
}
/// <summary>
/// Can solve Ax=b using LU factorization.
/// </summary>
[Test]
public void CanSolveUsingLU()
{
var matrix = _matrices["Square3x3"];
var a = new double[matrix.RowCount * matrix.RowCount];
Array.Copy(matrix.Data, a, a.Length);
var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
Provider.LUSolve(2, a, matrix.RowCount, b);
AssertHelpers.AlmostEqual(b[0], -1.477272727272726, 14);
AssertHelpers.AlmostEqual(b[1], -4.318181818181815, 14);
AssertHelpers.AlmostEqual(b[2], 3.068181818181816, 14);
AssertHelpers.AlmostEqual(b[3], -4.204545454545451, 14);
AssertHelpers.AlmostEqual(b[4], -12.499999999999989, 14);
AssertHelpers.AlmostEqual(b[5], 8.522727272727266, 14);
NotModified(matrix.RowCount, matrix.ColumnCount, a, matrix);
}
/// <summary>
/// Can solve Ax=b using LU factorization using a factored matrix.
/// </summary>
[Test]
public void CanSolveUsingLUOnFactoredMatrix()
{
var matrix = _matrices["Square3x3"];
var a = new double[matrix.RowCount * matrix.RowCount];
Array.Copy(matrix.Data, a, a.Length);
var ipiv = new int[matrix.RowCount];
Provider.LUFactor(a, matrix.RowCount, ipiv);
var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
Provider.LUSolveFactored(2, a, matrix.RowCount, ipiv, b);
AssertHelpers.AlmostEqual(b[0], -1.477272727272726, 14);
AssertHelpers.AlmostEqual(b[1], -4.318181818181815, 14);
AssertHelpers.AlmostEqual(b[2], 3.068181818181816, 14);
AssertHelpers.AlmostEqual(b[3], -4.204545454545451, 14);
AssertHelpers.AlmostEqual(b[4], -12.499999999999989, 14);
AssertHelpers.AlmostEqual(b[5], 8.522727272727266, 14);
}
/// <summary>
/// Checks to see if a matrix and array contain the same values.
/// </summary>
/// <param name="rows">number of rows.</param>
/// <param name="columns">number of columns.</param>
/// <param name="array">array to check.</param>
/// <param name="matrix">matrix to check against.</param>
private static void NotModified(int rows, int columns, IList<double> array, Matrix<double> matrix)
{
var index = 0;
for (var col = 0; col < columns; col++)
{
for (var row = 0; row < rows; row++)
{
Assert.AreEqual(array[index++], matrix[row, col]);
}
}
}
}
}

Write Preview
Loading…
Cancel
Save