|
|
|
@ -68,39 +68,39 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// <summary>
|
|
|
|
/// The default number of starting vectors.
|
|
|
|
/// </summary>
|
|
|
|
private const int DefaultNumberOfStartingVectors = 50; |
|
|
|
|
|
|
|
const int DefaultNumberOfStartingVectors = 50; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// The status used if there is no status, i.e. the solver hasn't run yet and there is no
|
|
|
|
/// iterator.
|
|
|
|
/// </summary>
|
|
|
|
private static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); |
|
|
|
static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// The preconditioner that will be used. Can be set to <see langword="null" />, in which case the default
|
|
|
|
/// pre-conditioner will be used.
|
|
|
|
/// </summary>
|
|
|
|
private IPreConditioner _preconditioner; |
|
|
|
IPreConditioner _preconditioner; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// The iterative process controller.
|
|
|
|
/// </summary>
|
|
|
|
private IIterator _iterator; |
|
|
|
IIterator _iterator; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// The collection of starting vectors which are used as the basis for the Krylov sub-space.
|
|
|
|
/// </summary>
|
|
|
|
private IList<Vector> _startingVectors; |
|
|
|
IList<Vector> _startingVectors; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// The number of starting vectors used by the algorithm
|
|
|
|
/// </summary>
|
|
|
|
private int _numberOfStartingVectors = DefaultNumberOfStartingVectors; |
|
|
|
int _numberOfStartingVectors = DefaultNumberOfStartingVectors; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Indicates if the user has stopped the solver.
|
|
|
|
/// </summary>
|
|
|
|
private bool _hasBeenStopped; |
|
|
|
bool _hasBeenStopped; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Initializes a new instance of the <see cref="MlkBiCgStab"/> class.
|
|
|
|
@ -109,7 +109,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// When using this constructor the solver will use the <see cref="IIterator"/> with
|
|
|
|
/// the standard settings and a default preconditioner.
|
|
|
|
/// </remarks>
|
|
|
|
public MlkBiCgStab() : this(null, null) |
|
|
|
public MlkBiCgStab() |
|
|
|
: this(null, null) |
|
|
|
{ |
|
|
|
} |
|
|
|
|
|
|
|
@ -132,7 +133,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// </para>
|
|
|
|
/// </remarks>
|
|
|
|
/// <param name="iterator">The <see cref="IIterator"/> that will be used to monitor the iterative process.</param>
|
|
|
|
public MlkBiCgStab(IIterator iterator) : this(null, iterator) |
|
|
|
public MlkBiCgStab(IIterator iterator) |
|
|
|
: this(null, iterator) |
|
|
|
{ |
|
|
|
} |
|
|
|
|
|
|
|
@ -144,7 +146,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// the standard settings.
|
|
|
|
/// </remarks>
|
|
|
|
/// <param name="preconditioner">The <see cref="IPreConditioner"/> that will be used to precondition the matrix equation.</param>
|
|
|
|
public MlkBiCgStab(IPreConditioner preconditioner) : this(preconditioner, null) |
|
|
|
public MlkBiCgStab(IPreConditioner preconditioner) |
|
|
|
: this(preconditioner, null) |
|
|
|
{ |
|
|
|
} |
|
|
|
|
|
|
|
@ -181,10 +184,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
public int NumberOfStartingVectors |
|
|
|
{ |
|
|
|
[DebuggerStepThrough] |
|
|
|
get |
|
|
|
{ |
|
|
|
return _numberOfStartingVectors; |
|
|
|
} |
|
|
|
get { return _numberOfStartingVectors; } |
|
|
|
|
|
|
|
[DebuggerStepThrough] |
|
|
|
set |
|
|
|
@ -231,10 +231,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
public IList<Vector> StartingVectors |
|
|
|
{ |
|
|
|
[DebuggerStepThrough] |
|
|
|
get |
|
|
|
{ |
|
|
|
return _startingVectors; |
|
|
|
} |
|
|
|
get { return _startingVectors; } |
|
|
|
|
|
|
|
[DebuggerStepThrough] |
|
|
|
set |
|
|
|
@ -256,10 +253,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
public ICalculationStatus IterationResult |
|
|
|
{ |
|
|
|
[DebuggerStepThrough] |
|
|
|
get |
|
|
|
{ |
|
|
|
return (_iterator != null) ? _iterator.Status : DefaultStatus; |
|
|
|
} |
|
|
|
get { return (_iterator != null) ? _iterator.Status : DefaultStatus; } |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
@ -348,7 +342,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
{ |
|
|
|
_preconditioner = new UnitPreconditioner(); |
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
_preconditioner.Initialize(matrix); |
|
|
|
|
|
|
|
// Choose an initial guess x_0
|
|
|
|
@ -402,7 +396,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
Vector 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]); |
|
|
|
@ -420,14 +414,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
matrix.Multiply(gtemp, w[k - 1]); |
|
|
|
|
|
|
|
// c_((j-1)k+k) = q^T_1 w_((j-1)k+k)
|
|
|
|
c[k - 1] = _startingVectors[0].DotProduct(w[k - 1]); |
|
|
|
c[k - 1] = _startingVectors[0].ConjugateDotProduct(w[k - 1]); |
|
|
|
if (c[k - 1].Real.AlmostEqual(0, 1) && c[k - 1].Imaginary.AlmostEqual(0, 1)) |
|
|
|
{ |
|
|
|
throw new Exception("Iterative solver experience a numerical break down"); |
|
|
|
} |
|
|
|
|
|
|
|
// alpha_(jk+1) = q^T_1 r_((j-1)k+k) / c_((j-1)k+k)
|
|
|
|
var alpha = _startingVectors[0].DotProduct(residuals) / c[k - 1]; |
|
|
|
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); |
|
|
|
@ -439,7 +433,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
|
|
|
|
// rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2
|
|
|
|
matrix.Multiply(temp1, temp); |
|
|
|
var rho = temp.DotProduct(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).
|
|
|
|
@ -449,7 +443,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
rho = 1.0f; |
|
|
|
} |
|
|
|
|
|
|
|
rho = -u.DotProduct(temp) / rho; |
|
|
|
rho = -u.ConjugateDotProduct(temp)/rho; |
|
|
|
|
|
|
|
// r_(jk+1) = rho_(j+1) A u~_(jk+1) + u_(jk+1)
|
|
|
|
u.CopyTo(residuals); |
|
|
|
@ -502,7 +496,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
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].DotProduct(zd) / c[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); |
|
|
|
@ -521,7 +515,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
beta = rho * c[k - 1]; |
|
|
|
beta = rho*c[k - 1]; |
|
|
|
if (beta.Real.AlmostEqual(0, 1) && beta.Imaginary.AlmostEqual(0, 1)) |
|
|
|
{ |
|
|
|
throw new Exception("Iterative solver experience a numerical break down"); |
|
|
|
@ -530,7 +524,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.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, temp2); |
|
|
|
residuals.Add(temp2, temp); |
|
|
|
beta = -_startingVectors[0].DotProduct(temp) / beta; |
|
|
|
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); |
|
|
|
@ -550,7 +544,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
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].DotProduct(zd) / c[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); |
|
|
|
@ -573,14 +567,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
if (i < k - 1) |
|
|
|
{ |
|
|
|
// c_(jk+1) = q^T_i+1 d_(jk+i)
|
|
|
|
c[i] = _startingVectors[i + 1].DotProduct(d[i]); |
|
|
|
c[i] = _startingVectors[i + 1].ConjugateDotProduct(d[i]); |
|
|
|
if (c[i].Real.AlmostEqual(0, 1) && c[i].Imaginary.AlmostEqual(0, 1)) |
|
|
|
{ |
|
|
|
throw new Exception("Iterative solver experience a numerical break down"); |
|
|
|
} |
|
|
|
|
|
|
|
// alpha_(jk+i+1) = q^T_(i+1) u_(jk+i) / c_(jk+i)
|
|
|
|
alpha = _startingVectors[i + 1].DotProduct(u) / c[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); |
|
|
|
@ -591,7 +585,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
_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); |
|
|
|
gtemp.Multiply(rho*alpha, temp); |
|
|
|
xtemp.Add(temp, temp2); |
|
|
|
temp2.CopyTo(xtemp); |
|
|
|
|
|
|
|
@ -599,7 +593,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
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); |
|
|
|
w[i].Multiply(-rho*alpha, temp); |
|
|
|
residuals.Add(temp, temp2); |
|
|
|
temp2.CopyTo(residuals); |
|
|
|
|
|
|
|
@ -626,7 +620,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// <param name="maximumNumberOfStartingVectors">Maximum number</param>
|
|
|
|
/// <param name="numberOfVariables">Number of variables</param>
|
|
|
|
/// <returns>Number of starting vectors to create</returns>
|
|
|
|
private static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables) |
|
|
|
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)); |
|
|
|
@ -643,7 +637,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// the <paramref name="numberOfVariables"/> is smaller than
|
|
|
|
/// the <paramref name="maximumNumberOfStartingVectors"/>.
|
|
|
|
/// </returns>
|
|
|
|
private static IList<Vector> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) |
|
|
|
static IList<Vector> 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
|
|
|
|
@ -662,7 +656,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
var samplesIm = distribution.Samples().Take(matrix.RowCount).ToArray(); |
|
|
|
for (int j = 0; j < matrix.RowCount; j++) |
|
|
|
{ |
|
|
|
samples[j] = new Complex32((float)samplesRe[j], (float)samplesIm[j]); |
|
|
|
samples[j] = new Complex32((float) samplesRe[j], (float) samplesIm[j]); |
|
|
|
} |
|
|
|
|
|
|
|
// Set the column
|
|
|
|
@ -677,10 +671,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
var result = new List<Vector>(); |
|
|
|
for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) |
|
|
|
{ |
|
|
|
result.Add((Vector)orthogonalMatrix.Column(i)); |
|
|
|
|
|
|
|
result.Add((Vector) orthogonalMatrix.Column(i)); |
|
|
|
|
|
|
|
// Normalize the result vector
|
|
|
|
result[i].Multiply(1 / result[i].L2Norm().Real, result[i]); |
|
|
|
result[i].Multiply(1/result[i].L2Norm().Real, result[i]); |
|
|
|
} |
|
|
|
|
|
|
|
return result; |
|
|
|
@ -692,7 +686,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// <param name="arraySize">Number of vectors</param>
|
|
|
|
/// <param name="vectorSize">Size of each vector</param>
|
|
|
|
/// <returns>Array of random vectors</returns>
|
|
|
|
private static Vector[] CreateVectorArray(int arraySize, int vectorSize) |
|
|
|
static Vector[] CreateVectorArray(int arraySize, int vectorSize) |
|
|
|
{ |
|
|
|
var result = new Vector[arraySize]; |
|
|
|
for (var i = 0; i < result.Length; i++) |
|
|
|
@ -710,7 +704,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// <param name="residual">Residual <see cref="Vector"/> data.</param>
|
|
|
|
/// <param name="x">x <see cref="Vector"/> data.</param>
|
|
|
|
/// <param name="b">b <see cref="Vector"/> data.</param>
|
|
|
|
private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) |
|
|
|
static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) |
|
|
|
{ |
|
|
|
// -Ax = residual
|
|
|
|
matrix.Multiply(x, residual); |
|
|
|
@ -728,7 +722,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// <param name="source">Source <see cref="Vector"/>.</param>
|
|
|
|
/// <param name="residuals">Residual <see cref="Vector"/>.</param>
|
|
|
|
/// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
|
|
|
|
private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) |
|
|
|
bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) |
|
|
|
{ |
|
|
|
if (_hasBeenStopped) |
|
|
|
{ |
|
|
|
@ -764,7 +758,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
throw new ArgumentNullException("input"); |
|
|
|
} |
|
|
|
|
|
|
|
var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount); |
|
|
|
var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount); |
|
|
|
Solve(matrix, input, result); |
|
|
|
return result; |
|
|
|
} |
|
|
|
@ -800,7 +794,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
|
|
|
|
for (var column = 0; column < input.ColumnCount; column++) |
|
|
|
{ |
|
|
|
var solution = Solve(matrix, (Vector)input.Column(column)); |
|
|
|
var solution = Solve(matrix, (Vector) input.Column(column)); |
|
|
|
foreach (var element in solution.GetIndexedEnumerator()) |
|
|
|
{ |
|
|
|
result.At(element.Item1, column, element.Item2); |
|
|
|
|
Christoph Ruegg
13 years ago