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

merged andriy's changes

pull/36/head
Marcus Cuda 16 years ago
parent
1a8cce73c6
commit
b17bbf47cc
15 changed files with 442 additions and 275 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. 129
      src/Numerics/LinearAlgebra/Double/DenseVector.cs
  2. 4
      src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs
  3. 2
      src/Numerics/LinearAlgebra/Double/Factorization/QR.cs
  4. 31
      src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs
  5. 51
      src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs
  6. 57
      src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs
  7. 25
      src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs
  8. 4
      src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/Ilutp.cs
  9. 54
      src/Numerics/LinearAlgebra/Double/SparseVector.cs
  10. 51
      src/Numerics/LinearAlgebra/Double/Vector.cs
  11. 4
      src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs
  12. 76
      src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs
  13. 76
      src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs
  14. 76
      src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs
  15. 77
      src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs

129
src/Numerics/LinearAlgebra/Double/DenseVector.cs
View File

@ -29,6 +29,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
using System;
using System.Collections.Generic;
using System.Globalization;
using System.Linq;
using Distributions;
using NumberTheory;
using Properties;
@ -303,10 +304,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 0.0)
{
return this.Clone();
return Clone();
}
var copy = (DenseVector) this.Clone();
var copy = (DenseVector)Clone();
CommonParallel.For(
0,
Data.Length,
@ -343,7 +344,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
CommonParallel.For(
0,
Data.Length,
index => dense.Data[index] = this.Data[index] + scalar);
index => dense.Data[index] = Data[index] + scalar);
}
}
@ -372,12 +373,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
return base.Add(other);
}
else
{
var copy = (DenseVector)this.Clone();
Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, 1.0, denseVector.Data);
return copy;
}
var copy = (DenseVector)Clone();
Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, 1.0, denseVector.Data);
return copy;
}
/// <summary>
@ -408,7 +407,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = this.Add(other);
var tmp = Add(other);
tmp.CopyTo(result);
}
else
@ -425,7 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
CommonParallel.For(
0,
Data.Length,
index => result[index] = this.Data[index] + other[index]);
index => result[index] = Data[index] + other[index]);
}
}
}
@ -484,10 +483,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 0.0)
{
return this.Clone();
return Clone();
}
var copy = (DenseVector)this.Clone();
var copy = (DenseVector)Clone();
CommonParallel.For(
0,
Data.Length,
@ -517,14 +516,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var dense = result as DenseVector;
if (dense == null)
{
base.Add(scalar, result);
base.Subtract(scalar, result);
}
else
{
CommonParallel.For(
0,
Data.Length,
index => dense.Data[index] = this.Data[index] - scalar);
index => dense.Data[index] = Data[index] - scalar);
}
}
@ -553,12 +552,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
return base.Subtract(other);
}
else
{
var copy = (DenseVector)this.Clone();
Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, -1.0, denseVector.Data);
return copy;
}
var copy = (DenseVector)Clone();
Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, -1.0, denseVector.Data);
return copy;
}
/// <summary>
@ -589,8 +586,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = result.CreateVector(result.Count);
Subtract(other, tmp);
var tmp = Subtract(other);
tmp.CopyTo(result);
}
else
@ -607,7 +603,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
CommonParallel.For(
0,
Data.Length,
index => result[index] = this.Data[index] - other[index]);
index => result[index] = Data[index] - other[index]);
}
}
}
@ -681,10 +677,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 1.0)
{
return this.Clone();
return Clone();
}
var copy = (DenseVector)this.Clone();
var copy = (DenseVector)Clone();
Control.LinearAlgebraProvider.ScaleArray(scalar, copy.Data);
return copy;
}
@ -732,7 +728,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentNullException("leftSide");
}
return (DenseVector) leftSide.Multiply(rightSide);
return (DenseVector)leftSide.Multiply(rightSide);
}
/// <summary>
@ -749,7 +745,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentNullException("rightSide");
}
return (DenseVector) rightSide.Multiply(leftSide);
return (DenseVector)rightSide.Multiply(leftSide);
}
/// <summary>
@ -794,7 +790,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentNullException("leftSide");
}
return (DenseVector) leftSide.Multiply(1.0 / rightSide);
return (DenseVector)leftSide.Multiply(1.0 / rightSide);
}
/// <summary>
@ -1013,15 +1009,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
return base.PointwiseMultiply(other);
}
else
{
var copy = (DenseVector)this.Clone();
CommonParallel.For(
0,
Count,
index => copy[index] *= other[index]);
return copy;
}
var copy = (DenseVector)Clone();
CommonParallel.For(
0,
Count,
index => copy[index] *= other[index]);
return copy;
}
/// <summary>
@ -1057,8 +1051,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = result.CreateVector(result.Count);
PointwiseMultiply(other, tmp);
var tmp = PointwiseMultiply(other);
tmp.CopyTo(result);
}
else
@ -1073,7 +1066,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
CommonParallel.For(
0,
Data.Length,
index => dense.Data[index] = this.Data[index] * other[index]);
index => dense.Data[index] = Data[index] * other[index]);
}
}
}
@ -1103,15 +1096,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
return base.PointwiseMultiply(other);
}
else
{
var copy = (DenseVector)this.Clone();
CommonParallel.For(
0,
Count,
index => copy[index] /= other[index]);
return copy;
}
var copy = (DenseVector)Clone();
CommonParallel.For(
0,
Count,
index => copy[index] /= other[index]);
return copy;
}
/// <summary>
@ -1147,7 +1138,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = this.PointwiseDivide(other);
var tmp = PointwiseDivide(other);
tmp.CopyTo(result);
}
else
@ -1162,7 +1153,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
CommonParallel.For(
0,
Data.Length,
index => dense.Data[index] = this.Data[index] / other[index]);
index => dense.Data[index] = Data[index] / other[index]);
}
}
}
@ -1279,41 +1270,35 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
throw new ArgumentOutOfRangeException("p");
}
else if (1.0 == p)
if (1.0 == p)
{
return CommonParallel.Aggregate(
0,
Count,
index => Math.Abs(Data[index]));
}
else if (2.0 == p)
{
var sum = 0.0;
for (var i = 0; i < Data.Length; i++)
{
sum = SpecialFunctions.Hypotenuse(sum, Data[i]);
}
return sum;
if (2.0 == p)
{
return Data.Aggregate(0.0, SpecialFunctions.Hypotenuse);
}
else if (Double.IsPositiveInfinity(p))
if (Double.IsPositiveInfinity(p))
{
return CommonParallel.Select(
0,
Count,
(index, localData) => localData = Math.Max(localData, Math.Abs(Data[index])),
(index, localData) => Math.Max(localData, Math.Abs(Data[index])),
Math.Max);
}
else
{
var sum = CommonParallel.Aggregate(
0,
Count,
index => Math.Pow(Math.Abs(Data[index]), p));
return Math.Pow(sum, 1.0 / p);
}
var sum = CommonParallel.Aggregate(
0,
Count,
index => Math.Pow(Math.Abs(Data[index]), p));
return Math.Pow(sum, 1.0 / p);
}
#endregion

4
src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs
View File

@ -165,13 +165,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
var dinput = input as DenseVector;
if (dinput == null)
{
throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment.");
throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment.");
}
var dresult = result as DenseVector;
if (dresult == null)
{
throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment.");
throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment.");
}
Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, 1, dresult.Data);

2
src/Numerics/LinearAlgebra/Double/Factorization/QR.cs
View File

@ -97,7 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
}
/// <summary>
/// Gets the determinant of the matrix for which the QR matrix was computed.
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
/// </summary>
public virtual double Determinant
{

31
src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs
View File

@ -290,8 +290,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
Vector vecSdash = new DenseVector(residuals.Count);
Vector t = new DenseVector(residuals.Count);
Vector mult = new DenseVector(residuals.Count);
// create some temporary double variables that are needed
// to hold values in between iterations
double currentRho = 0;
@ -319,11 +317,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
var beta = (currentRho / oldRho) * (alpha / omega);
// p_i = r_(i-1) + beta_(i-1)(p_(i-1) - omega_(i-1) * nu_(i-1))
nu.Multiply(-omega, mult);
vecP.Add(mult);
vecP.Multiply(beta);
vecP.Add(residuals);
vecP.Add(nu.Multiply(-omega), vecP);
vecP.Multiply(beta, vecP);
vecP.Add(residuals, vecP);
}
else
{
@ -341,8 +338,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
alpha = currentRho * 1 / tempResiduals.DotProduct(nu);
// s = r_(i-1) - alpha_i nu_i
nu.Multiply(-alpha, mult);
residuals.Add(mult, vecS);
residuals.Add(nu.Multiply(-alpha), vecS);
// Check if we're converged. If so then stop. Otherwise continue;
// Calculate the temporary result.
@ -350,8 +346,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
// temp. Others will be used in the calculation later on.
// x_i = x_(i-1) + alpha_i * p^_i + s^_i
vecPdash.Multiply(alpha, temp);
temp.Add(vecSdash);
temp.Add(result);
temp.Add(vecSdash, temp);
temp.Add(result, temp);
// Check convergence and stop if we are converged.
if (!ShouldContinue(iterationNumber, temp, input, vecS))
@ -383,14 +379,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
omega = t.DotProduct(vecS) / t.DotProduct(t);
// x_i = x_(i-1) + alpha_i p^ + omega_i s^
vecSdash.Multiply(omega, mult);
result.Add(mult);
result.Add(vecSdash.Multiply(omega), result);
vecPdash.Multiply(alpha, mult);
result.Add(mult);
//vecPdash.Multiply(alpha, vecPdash);
result.Add(vecPdash.Multiply(alpha), result);
t.Multiply(-omega, residuals);
residuals.Add(vecS);
residuals.Add(vecS, residuals);
// for continuation it is necessary that omega_i != 0.0
// If omega is only 1 ULP from zero then we fail.
@ -427,10 +422,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
matrix.Multiply(x, residual);
// Do not use residual = residual.Negate() because it creates another object
residual.Multiply(-1);
residual.Multiply(-1, residual);
// residual + b
residual.Add(b);
residual.Add(b, residual);
}
/// <summary>

51
src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs
View File

@ -359,7 +359,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
Vector z = new DenseVector(residuals.Count);
Vector temp = new DenseVector(residuals.Count);
Vector mult = new DenseVector(residuals.Count);
// for (k = 0, 1, .... )
var iterationNumber = 0;
@ -368,8 +367,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
// p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1))
p.Subtract(u, temp);
temp.Multiply(beta, mult);
residuals.Add(mult, p);
residuals.Add(temp.Multiply(beta), p);
// Solve M b_k = p_k
_preconditioner.Approximate(p, temp);
@ -384,15 +382,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
s.Subtract(w, temp);
t.Subtract(residuals, y);
temp.Multiply(alpha, mult);
y.Add(mult);
y.Add(temp.Multiply(alpha), y);
// Store the old value of t in t0
t.CopyTo(t0);
// t_k = r_k - alpha_k s_k
s.Multiply(-alpha, mult);
residuals.Add(mult, t);
residuals.Add(s.Multiply(-alpha), t);
// Solve M d_k = t_k
_preconditioner.Approximate(t, temp);
@ -450,48 +446,39 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
}
// u_k = sigma_k s_k + eta_k (t_(k-1) - r_k + beta_(k-1) u_(k-1))
u.Multiply(beta, mult);
t0.Add(mult, temp);
temp.Subtract(residuals);
temp.Multiply(eta);
t0.Add(u.Multiply(beta), temp);
s.Multiply(sigma, mult);
temp.Add(mult, u);
temp.Subtract(residuals, temp);
temp.Multiply(eta, temp);
// z_k = sigma_k r_k +_ eta_k z_(k-1) - alpha_k u_k
z.Multiply(eta);
temp.Add(s.Multiply(sigma), u);
u.Multiply(-alpha, mult);
z.Add(mult);
// z_k = sigma_k r_k +_ eta_k z_(k-1) - alpha_k u_k
z.Multiply(eta, z);
z.Add(u.Multiply(-alpha), z);
residuals.Multiply(sigma, mult);
z.Add(mult);
z.Add(residuals.Multiply(sigma), z);
// x_(k+1) = x_k + alpha_k p_k + z_k
p.Multiply(alpha, mult);
xtemp.Add(mult);
xtemp.Add(z);
xtemp.Add(p.Multiply(alpha), xtemp);
xtemp.Add(z, 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, mult);
t.Add(mult, residuals);
t.Add(y.Multiply(-eta), residuals);
c.Multiply(-sigma, mult);
residuals.Add(mult);
residuals.Add(c.Multiply(-sigma), 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.AlmostEqual(0, 1)) ? alpha / sigma * rdash.DotProduct(residuals) / rdash.DotProduct(t0) : 0;
// w_k = c_k + beta_k s_k
s.Multiply(beta, mult);
c.Add(mult, w);
c.Add(s.Multiply(beta), w);
// Get the real value
_preconditioner.Approximate(xtemp, result);
@ -520,10 +507,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
{
// -Ax = residual
matrix.Multiply(x, residual);
residual.Multiply(-1);
residual.Multiply(-1, residual);
// residual + b
residual.Add(b);
residual.Add(b, residual);
}
/// <summary>

57
src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs
View File

@ -395,8 +395,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
Vector zg = new DenseVector(residuals.Count);
Vector zw = new DenseVector(residuals.Count);
Vector mult = new DenseVector(residuals.Count);
var d = CreateVectorArray(_startingVectors.Count, residuals.Count);
// g_0 = r_0
@ -426,8 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
var alpha = _startingVectors[0].DotProduct(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, mult);
residuals.Add(mult, u);
residuals.Add(w[k - 1].Multiply(-alpha), u);
// SOLVE M u~_(jk+1) = u_(jk+1)
_preconditioner.Approximate(u, temp1);
@ -448,18 +445,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
rho = -u.DotProduct(temp) / rho;
// 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, mult);
xtemp.Add(mult);
xtemp.Add(utemp.Multiply(-rho), xtemp);
gtemp.Multiply(alpha);
xtemp.Add(gtemp);
gtemp.Multiply(alpha, gtemp);
xtemp.Add(gtemp, xtemp);
// r_(jk+1) = rho_(j+1) A u~_(jk+1) + u_(jk+1)
u.CopyTo(residuals);
// Reuse temp
temp.Multiply(rho);
residuals.Add(temp);
temp.Multiply(rho, temp);
residuals.Add(temp, residuals);
// Check convergence and stop if we are converged.
if (!ShouldContinue(iterationNumber, xtemp, input, residuals))
@ -498,16 +494,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
beta = -_startingVectors[s + 1].DotProduct(zd) / c[s];
// z_d = z_d + beta^(jk+i)_((j-1)k+s) d_((j-1)k+s)
d[s].Multiply(beta, mult);
zd.Add(mult);
zd.Add(d[s].Multiply(beta), zd);
// z_g = z_g + beta^(jk+i)_((j-1)k+s) g_((j-1)k+s)
g[s].Multiply(beta, mult);
zg.Add(mult);
zg.Add(g[s].Multiply(beta), zg);
// z_w = z_w + beta^(jk+i)_((j-1)k+s) w_((j-1)k+s)
w[s].Multiply(beta, mult);
zw.Add(mult);
zw.Add(w[s].Multiply(beta), zw);
}
}
@ -518,18 +511,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
}
// 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, mult);
residuals.Add(mult, temp);
residuals.Add(zw.Multiply(rho), temp);
beta = -_startingVectors[0].DotProduct(temp) / beta;
// z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k)
g[k - 1].Multiply(beta, mult);
zg.Add(mult);
zg.Add(g[k - 1].Multiply(beta), 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, mult);
zw.Add(mult);
zw.Multiply(rho);
zw.Add(w[k - 1].Multiply(beta), zw);
zw.Multiply(rho, zw);
// z_d = r_(jk+i) + z_w
residuals.Add(zw, zd);
@ -541,12 +531,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
beta = -_startingVectors[s + 1].DotProduct(zd) / c[s];
// z_d = z_d + beta^(jk+i)_(jk+s) * d_(jk+s)
d[s].Multiply(beta, mult);
zd.Add(mult);
zd.Add(d[s].Multiply(beta), zd);
// z_g = z_g + beta^(jk+i)_(jk+s) * g_(jk+s)
g[s].Multiply(beta, mult);
zg.Add(mult);
zg.Add(g[s].Multiply(beta), zg);
}
// d_(jk+i) = z_d - u_(jk+i)
@ -569,22 +557,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
alpha = _startingVectors[i + 1].DotProduct(u) / c[i];
// u_(jk+i+1) = u_(jk+i) - alpha_(jk+i+1) d_(jk+i)
d[i].Multiply(-alpha, mult);
u.Add(mult);
u.Add(d[i].Multiply(-alpha), 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, mult);
xtemp.Add(mult);
xtemp.Add(gtemp.Multiply(rho * alpha), 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, mult);
residuals.Add(mult);
residuals.Add(w[i].Multiply(-rho * alpha), residuals);
// We can check the residuals here if they're close
if (!ShouldContinue(iterationNumber, xtemp, input, residuals))
@ -657,7 +642,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
result.Add(orthogonalMatrix.Column(i));
// Normalize the result vector
result[i].Multiply(1 / result[i].Norm(2));
result[i].Multiply(1 / result[i].Norm(2), result[i]);
}
return result;
@ -691,10 +676,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
{
// -Ax = residual
matrix.Multiply(x, residual);
residual.Multiply(-1);
residual.Multiply(-1, residual);
// residual + b
residual.Add(b);
residual.Add(b, residual);
}
/// <summary>

25
src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs
View File

@ -273,8 +273,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
// Temp vectors
var temp = new DenseVector(input.Count);
var mult = new DenseVector(input.Count);
// Initialize
var startNorm = input.Norm(2);
@ -316,8 +315,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
alpha = rho / sigma;
// yOdd = yEven - alpha * v
v.Multiply(-alpha, mult);
yeven.Add(mult, yodd);
yeven.Add(v.Multiply(-alpha), yodd);
// Solve M temp = yOdd
_preconditioner.Approximate(yodd, temp);
@ -333,8 +331,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
var yinternal = IsEven(iterationNumber) ? yeven : yodd;
// pseudoResiduals = pseudoResiduals - alpha * uOdd
uinternal.Multiply(-alpha, mult);
pseudoResiduals.Add(mult);
pseudoResiduals.Add(uinternal.Multiply(-alpha), pseudoResiduals);
// d = yOdd + theta * theta * eta / alpha * d
d.Multiply(theta * theta * eta / alpha, temp);
@ -351,8 +348,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
eta = c * c * alpha;
// x = x + eta * d
d.Multiply(eta, mult);
x.Add(mult);
x.Add(d.Multiply(eta), x);
// Check convergence and see if we can bail
if (!ShouldContinue(iterationNumber, result, input, pseudoResiduals))
@ -390,8 +386,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
rho = rhoNew;
// yOdd = pseudoResiduals + beta * yOdd
yodd.Multiply(beta, mult);
pseudoResiduals.Add(mult, yeven);
pseudoResiduals.Add(yodd.Multiply(beta), yeven);
// Solve M temp = yOdd
_preconditioner.Approximate(yeven, temp);
@ -400,11 +395,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
matrix.Multiply(temp, ueven);
// v = uEven + beta * (uOdd + beta * v)
v.Multiply(beta, mult);
uodd.Add(mult, temp);
uodd.Add(v.Multiply(beta), temp);
temp.Multiply(beta, mult);
ueven.Add(mult, v);
ueven.Add(temp.Multiply(beta), v);
}
// Calculate the real values
@ -425,10 +418,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
{
// -Ax = residual
matrix.Multiply(x, residual);
residual.Multiply(-1);
residual.Multiply(-1, residual);
// residual + b
residual.Add(b);
residual.Add(b, residual);
}
/// <summary>

4
src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/Ilutp.cs
View File

@ -426,8 +426,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
rowVector[k] = 0.0;
}
rowVector.Multiply(workVector[j]);
workVector.Subtract(rowVector);
rowVector.Multiply(workVector[j], rowVector);
workVector.Subtract(rowVector, workVector);
}
}
}

54
src/Numerics/LinearAlgebra/Double/SparseVector.cs
View File

@ -351,14 +351,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 0.0)
{
return this.Clone();
return Clone();
}
var copy = (SparseVector)this.Clone();
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] += scalar;
}
return copy;
}
@ -409,12 +410,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
return base.Add(other);
}
else
{
var copy = (SparseVector)this.Clone();
copy.AddScaledSparseVector(1.0, sparseVector);
return copy;
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(1.0, sparseVector);
return copy;
}
/// <summary>
@ -514,7 +513,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = this.Add(other);
var tmp = Add(other);
tmp.CopyTo(result);
}
else
@ -593,14 +592,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 0.0)
{
return this.Clone();
return Clone();
}
var copy = (SparseVector)this.Clone();
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] -= scalar;
}
return copy;
}
@ -651,12 +651,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
return base.Subtract(other);
}
else
{
var copy = (SparseVector)this.Clone();
copy.AddScaledSparseVector(-1.0, sparseVector);
return copy;
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(-1.0, sparseVector);
return copy;
}
/// <summary>
@ -687,7 +685,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = this.Subtract(other);
var tmp = Subtract(other);
tmp.CopyTo(result);
}
else
@ -789,18 +787,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 1.0)
{
return this.Clone();
return Clone();
}
else if (scalar == 0)
if (scalar == 0)
{
var copy = this.Clone();
var copy = Clone();
copy.Clear(); // Set array empty
return copy;
}
else
{
var copy = (SparseVector)this.Clone();
var copy = (SparseVector)Clone();
Control.LinearAlgebraProvider.ScaleArray(scalar, copy._nonZeroValues);
return copy;
}
@ -912,7 +910,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentNullException("leftSide");
}
return (SparseVector) leftSide.Multiply(1.0 / rightSide);
return (SparseVector)leftSide.Multiply(1.0 / rightSide);
}
/// <summary>
@ -1103,12 +1101,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
}
var copy = new SparseVector(Count);
for (int i = 0; i < this._nonZeroIndices.Length; i++)
for (int i = 0; i < _nonZeroIndices.Length; i++)
{
var d = this._nonZeroValues[i] * other[this._nonZeroIndices[i]];
var d = _nonZeroValues[i] * other[_nonZeroIndices[i]];
if (d != 0.0)
{
copy[this._nonZeroIndices[i]] = d;
copy[_nonZeroIndices[i]] = d;
}
}
@ -1238,7 +1236,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (Double.IsPositiveInfinity(p))
{
return CommonParallel.Select(0, NonZerosCount, (index, localData) => localData = Math.Max(localData, Math.Abs(_nonZeroValues[index])), Math.Max);
return CommonParallel.Select(0, NonZerosCount, (index, localData) => Math.Max(localData, Math.Abs(_nonZeroValues[index])), Math.Max);
}
var sum = CommonParallel.Aggregate(

51
src/Numerics/LinearAlgebra/Double/Vector.cs
View File

@ -125,10 +125,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 0.0)
{
return this.Clone();
return Clone();
}
var copy = this.Clone();
var copy = Clone();
CommonParallel.For(
0,
Count,
@ -185,7 +185,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </remarks>
public virtual Vector Plus()
{
return this.Clone();
return Clone();
}
/// <summary>
@ -213,7 +213,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var copy = this.Clone();
var copy = Clone();
CommonParallel.For(
0,
Count,
@ -256,7 +256,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = this.Add(other);
var tmp = Add(other);
tmp.CopyTo(result);
}
else
@ -279,10 +279,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 0.0)
{
return this.Clone();
return Clone();
}
var copy = this.Clone();
var copy = Clone();
CommonParallel.For(
0,
Count,
@ -367,7 +367,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var copy = this.Clone();
var copy = Clone();
CommonParallel.For(
0,
Count,
@ -434,10 +434,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 1.0)
{
return this.Clone();
return Clone();
}
var copy = this.Clone();
var copy = Clone();
CommonParallel.For(
0,
Count,
@ -530,7 +530,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
if (scalar == 1.0)
{
return this.Clone();
return Clone();
}
return Multiply(1.0 / scalar);
@ -593,7 +593,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var copy = this.Clone();
var copy = Clone();
CommonParallel.For(
0,
Count,
@ -634,7 +634,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = this.PointwiseMultiply(other);
var tmp = PointwiseMultiply(other);
tmp.CopyTo(result);
}
else
@ -665,7 +665,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var copy = this.Clone();
var copy = Clone();
CommonParallel.For(
0,
Count,
@ -706,7 +706,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = this.PointwiseDivide(other);
var tmp = PointwiseDivide(other);
tmp.CopyTo(result);
}
else
@ -1156,23 +1156,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
throw new ArgumentOutOfRangeException("p");
}
else if (Double.IsPositiveInfinity(p))
if (Double.IsPositiveInfinity(p))
{
return CommonParallel.Select(
0,
Count,
(index, localData) => localData = Math.Max(localData, Math.Abs(this[index])),
(index, localData) => Math.Max(localData, Math.Abs(this[index])),
Math.Max);
}
else
{
var sum = CommonParallel.Aggregate(
0,
Count,
index => Math.Pow(Math.Abs(this[index]), p));
return Math.Pow(sum, 1.0 / p);
}
var sum = CommonParallel.Aggregate(
0,
Count,
index => Math.Pow(Math.Abs(this[index]), p));
return Math.Pow(sum, 1.0 / p);
}
/// <summary>
@ -1192,7 +1191,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
}
var norm = Norm(p);
var clone = this.Clone();
var clone = Clone();
if (norm == 0.0)
{
return clone;

4
src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs
View File

@ -695,7 +695,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
[Row(-4, ExpectedException = typeof(ArgumentOutOfRangeException))]
public void NormalizeColumns(int pValue)
{
var matrix = this.TestMatrices["Singular3x3"];
var matrix = this.TestMatrices["Square4x4"];
var result = matrix.NormalizeColumns(pValue);
for (var j = 0; j < result.ColumnCount; j++)
{
@ -710,7 +710,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
[Row(-3, ExpectedException = typeof(ArgumentOutOfRangeException))]
public void NormalizeRows(int pValue)
{
var matrix = this.TestMatrices["Singular3x3"].NormalizeRows(pValue);
var matrix = this.TestMatrices["Square4x4"].NormalizeRows(pValue);
for (var i = 0; i < matrix.RowCount; i++)
{
var row = matrix.Row(i);

76
src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs
View File

@ -40,11 +40,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
public void SolveUnitMatrixAndBackMultiply()
{
// Create the identity matrix
Matrix matrix = new SparseMatrix(100);
for (var i = 0; i < matrix.RowCount; i++)
{
matrix[i, i] = 1.0;
}
Matrix matrix = SparseMatrix.Identity(100);
// Create the y vector
Vector y = new DenseVector(matrix.RowCount, 1);
@ -84,11 +80,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
public void SolveScaledUnitMatrixAndBackMultiply()
{
// Create the identity matrix
Matrix matrix = new SparseMatrix(100);
for (var i = 0; i < matrix.RowCount; i++)
{
matrix[i, i] = 1.0;
}
Matrix matrix = SparseMatrix.Identity(100);
// Scale it with a funny number
matrix.Multiply(System.Math.PI);
@ -201,5 +193,69 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
Assert.IsTrue(System.Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#05-" + i);
}
}
[Test]
[Row(4)]
[Row(8)]
[Row(10)]
[MultipleAsserts]
public void CanSolveForRandomVector(int order)
{
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var vectorb = MatrixLoader.GenerateRandomDenseVector(order);
var monitor = new Iterator(new IIterationStopCriterium[]
{
new IterationCountStopCriterium(1000),
new ResidualStopCriterium(1e-10),
});
var solver = new BiCgStab(monitor);
var resultx = solver.Solve(matrixA, vectorb);
Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
var bReconstruct = matrixA * resultx;
// Check the reconstruction.
for (var i = 0; i < order; i++)
{
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-7);
}
}
[Test]
[Row(4)]
[Row(8)]
[Row(10)]
[MultipleAsserts]
public void CanSolveForRandomMatrix(int order)
{
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var monitor = new Iterator(new IIterationStopCriterium[]
{
new IterationCountStopCriterium(1000),
new ResidualStopCriterium(1e-10)
});
var solver = new BiCgStab(monitor);
var matrixX = solver.Solve(matrixA, matrixB);
// The solution X row dimension is equal to the column dimension of A
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
// The solution X has the same number of columns as B
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
var matrixBReconstruct = matrixA * matrixX;
// Check the reconstruction.
for (var i = 0; i < matrixB.RowCount; i++)
{
for (var j = 0; j < matrixB.ColumnCount; j++)
{
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-7);
}
}
}
}
}

76
src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs
View File

@ -40,11 +40,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
public void SolveUnitMatrixAndBackMultiply()
{
// Create the identity matrix
Matrix matrix = new SparseMatrix(100);
for (var i = 0; i < matrix.RowCount; i++)
{
matrix[i, i] = 1.0;
}
Matrix matrix = SparseMatrix.Identity(100);
// Create the y vector
Vector y = new DenseVector(matrix.RowCount, 1);
@ -84,11 +80,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
public void SolveScaledUnitMatrixAndBackMultiply()
{
// Create the identity matrix
Matrix matrix = new SparseMatrix(100);
for (var i = 0; i < matrix.RowCount; i++)
{
matrix[i, i] = 1.0;
}
Matrix matrix = SparseMatrix.Identity(100);
// Scale it with a funny number
matrix.Multiply(System.Math.PI);
@ -203,5 +195,69 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
Assert.IsTrue(System.Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#05-" + i);
}
}
[Test]
[Row(4)]
[Row(8)]
[Row(10)]
[MultipleAsserts]
public void CanSolveForRandomVector(int order)
{
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var vectorb = MatrixLoader.GenerateRandomDenseVector(order);
var monitor = new Iterator(new IIterationStopCriterium[]
{
new IterationCountStopCriterium(1000),
new ResidualStopCriterium(1e-10),
});
var solver = new GpBiCg(monitor);
var resultx = solver.Solve(matrixA, vectorb);
Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
var bReconstruct = matrixA * resultx;
// Check the reconstruction.
for (var i = 0; i < order; i++)
{
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-7);
}
}
[Test]
[Row(4)]
[Row(8)]
[Row(10)]
[MultipleAsserts]
public void CanSolveForRandomMatrix(int order)
{
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var monitor = new Iterator(new IIterationStopCriterium[]
{
new IterationCountStopCriterium(1000),
new ResidualStopCriterium(1e-10)
});
var solver = new GpBiCg(monitor);
var matrixX = solver.Solve(matrixA, matrixB);
// The solution X row dimension is equal to the column dimension of A
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
// The solution X has the same number of columns as B
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
var matrixBReconstruct = matrixA * matrixX;
// Check the reconstruction.
for (var i = 0; i < matrixB.RowCount; i++)
{
for (var j = 0; j < matrixB.ColumnCount; j++)
{
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-7);
}
}
}
}
}

76
src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs
View File

@ -40,11 +40,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
public void SolveUnitMatrixAndBackMultiply()
{
// Create the identity matrix
Matrix matrix = new SparseMatrix(100);
for (var i = 0; i < matrix.RowCount; i++)
{
matrix[i, i] = 1.0;
}
Matrix matrix = SparseMatrix.Identity(100);
// Create the y vector
Vector y = new DenseVector(matrix.RowCount, 1);
@ -85,11 +81,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
public void SolveScaledUnitMatrixAndBackMultiply()
{
// Create the identity matrix
Matrix matrix = new SparseMatrix(100);
for (var i = 0; i < matrix.RowCount; i++)
{
matrix[i, i] = 1.0;
}
Matrix matrix = SparseMatrix.Identity(100);
// Scale it with a funny number
matrix.Multiply(System.Math.PI);
@ -201,5 +193,69 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
Assert.IsTrue(System.Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#05-" + i);
}
}
[Test]
[Row(4)]
[Row(8)]
[Row(10)]
[MultipleAsserts]
public void CanSolveForRandomVector(int order)
{
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var vectorb = MatrixLoader.GenerateRandomDenseVector(order);
var monitor = new Iterator(new IIterationStopCriterium[]
{
new IterationCountStopCriterium(1000),
new ResidualStopCriterium(1e-10),
});
var solver = new MlkBiCgStab(monitor);
var resultx = solver.Solve(matrixA, vectorb);
Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
var bReconstruct = matrixA * resultx;
// Check the reconstruction.
for (var i = 0; i < order; i++)
{
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-7);
}
}
[Test]
[Row(4)]
[Row(8)]
[Row(10)]
[MultipleAsserts]
public void CanSolveForRandomMatrix(int order)
{
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var monitor = new Iterator(new IIterationStopCriterium[]
{
new IterationCountStopCriterium(1000),
new ResidualStopCriterium(1e-10)
});
var solver = new MlkBiCgStab(monitor);
var matrixX = solver.Solve(matrixA, matrixB);
// The solution X row dimension is equal to the column dimension of A
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
// The solution X has the same number of columns as B
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
var matrixBReconstruct = matrixA * matrixX;
// Check the reconstruction.
for (var i = 0; i < matrixB.RowCount; i++)
{
for (var j = 0; j < matrixB.ColumnCount; j++)
{
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-7);
}
}
}
}
}

77
src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs
View File

@ -40,11 +40,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
public void SolveUnitMatrixAndBackMultiply()
{
// Create the identity matrix
Matrix matrix = new SparseMatrix(100);
for (var i = 0; i < matrix.RowCount; i++)
{
matrix[i, i] = 1.0;
}
Matrix matrix = SparseMatrix.Identity(100);
// Create the y vector
Vector y = new DenseVector(matrix.RowCount, 1);
@ -85,11 +81,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
public void SolveScaledUnitMatrixAndBackMultiply()
{
// Create the identity matrix
Matrix matrix = new SparseMatrix(100);
for (var i = 0; i < matrix.RowCount; i++)
{
matrix[i, i] = 1.0;
}
Matrix matrix = SparseMatrix.Identity(100);
// Scale it with a funny number
matrix.Multiply(System.Math.PI);
@ -200,5 +193,69 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative
Assert.IsTrue(System.Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#05-" + i);
}
}
[Test]
[Row(4)]
[Row(8)]
[Row(10)]
[MultipleAsserts]
public void CanSolveForRandomVector(int order)
{
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var vectorb = MatrixLoader.GenerateRandomDenseVector(order);
var monitor = new Iterator(new IIterationStopCriterium[]
{
new IterationCountStopCriterium(1000),
new ResidualStopCriterium(1e-10),
});
var solver = new TFQMR(monitor);
var resultx = solver.Solve(matrixA, vectorb);
Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
var bReconstruct = matrixA * resultx;
// Check the reconstruction.
for (var i = 0; i < order; i++)
{
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-7);
}
}
[Test]
[Row(4)]
[Row(8)]
[Row(10)]
[MultipleAsserts]
public void CanSolveForRandomMatrix(int order)
{
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order);
var monitor = new Iterator(new IIterationStopCriterium[]
{
new IterationCountStopCriterium(1000),
new ResidualStopCriterium(1e-10)
});
var solver = new TFQMR(monitor);
var matrixX = solver.Solve(matrixA, matrixB);
// The solution X row dimension is equal to the column dimension of A
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
// The solution X has the same number of columns as B
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
var matrixBReconstruct = matrixA * matrixX;
// Check the reconstruction.
for (var i = 0; i < matrixB.RowCount; i++)
{
for (var j = 0; j < matrixB.ColumnCount; j++)
{
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-7);
}
}
}
}
}

Write Preview
Loading…
Cancel
Save