|
|
|
@ -4,7 +4,7 @@ |
|
|
|
// http://github.com/mathnet/mathnet-numerics
|
|
|
|
// http://mathnetnumerics.codeplex.com
|
|
|
|
//
|
|
|
|
// Copyright (c) 2009-2010 Math.NET
|
|
|
|
// Copyright (c) 2009-2013 Math.NET
|
|
|
|
//
|
|
|
|
// Permission is hereby granted, free of charge, to any person
|
|
|
|
// obtaining a copy of this software and associated documentation
|
|
|
|
@ -38,10 +38,8 @@ using MathNet.Numerics.LinearAlgebra.Solvers; |
|
|
|
using MathNet.Numerics.LinearAlgebra.Solvers.Status; |
|
|
|
using MathNet.Numerics.Properties; |
|
|
|
|
|
|
|
namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers |
|
|
|
{ |
|
|
|
using Numerics; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// A Multiple-Lanczos Bi-Conjugate Gradient stabilized iterative matrix solver.
|
|
|
|
/// </summary>
|
|
|
|
@ -65,7 +63,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// solver.
|
|
|
|
/// </para>
|
|
|
|
/// </remarks>
|
|
|
|
public sealed class MlkBiCgStab : IIterativeSolver<Complex32> |
|
|
|
public sealed class MlkBiCgStab : IIterativeSolver<Numerics.Complex32> |
|
|
|
{ |
|
|
|
/// <summary>
|
|
|
|
/// The default number of starting vectors.
|
|
|
|
@ -82,17 +80,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// The preconditioner that will be used. Can be set to <see langword="null" />, in which case the default
|
|
|
|
/// pre-conditioner will be used.
|
|
|
|
/// </summary>
|
|
|
|
IPreConditioner<Complex32> _preconditioner; |
|
|
|
IPreConditioner<Numerics.Complex32> _preconditioner; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// The iterative process controller.
|
|
|
|
/// </summary>
|
|
|
|
IIterator<Complex32> _iterator; |
|
|
|
IIterator<Numerics.Complex32> _iterator; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// The collection of starting vectors which are used as the basis for the Krylov sub-space.
|
|
|
|
/// </summary>
|
|
|
|
IList<Vector<Complex32>> _startingVectors; |
|
|
|
IList<Vector<Numerics.Complex32>> _startingVectors; |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// The number of starting vectors used by the algorithm
|
|
|
|
@ -135,7 +133,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// </para>
|
|
|
|
/// </remarks>
|
|
|
|
/// <param name="iterator">The <see cref="IIterator{T}"/> that will be used to monitor the iterative process.</param>
|
|
|
|
public MlkBiCgStab(IIterator<Complex32> iterator) |
|
|
|
public MlkBiCgStab(IIterator<Numerics.Complex32> iterator) |
|
|
|
: this(null, iterator) |
|
|
|
{ |
|
|
|
} |
|
|
|
@ -148,7 +146,7 @@ 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<Complex32> preconditioner) |
|
|
|
public MlkBiCgStab(IPreConditioner<Numerics.Complex32> preconditioner) |
|
|
|
: this(preconditioner, null) |
|
|
|
{ |
|
|
|
} |
|
|
|
@ -170,7 +168,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// </remarks>
|
|
|
|
/// <param name="preconditioner">The <see cref="IPreConditioner"/> that will be used to precondition the matrix equation.</param>
|
|
|
|
/// <param name="iterator">The <see cref="IIterator{T}"/> that will be used to monitor the iterative process.</param>
|
|
|
|
public MlkBiCgStab(IPreConditioner<Complex32> preconditioner, IIterator<Complex32> iterator) |
|
|
|
public MlkBiCgStab(IPreConditioner<Numerics.Complex32> preconditioner, IIterator<Numerics.Complex32> iterator) |
|
|
|
{ |
|
|
|
_iterator = iterator; |
|
|
|
_preconditioner = preconditioner; |
|
|
|
@ -212,7 +210,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// Sets the <see cref="IPreConditioner"/> that will be used to precondition the iterative process.
|
|
|
|
/// </summary>
|
|
|
|
/// <param name="preconditioner">The preconditioner.</param>
|
|
|
|
public void SetPreconditioner(IPreConditioner<Complex32> preconditioner) |
|
|
|
public void SetPreconditioner(IPreConditioner<Numerics.Complex32> preconditioner) |
|
|
|
{ |
|
|
|
_preconditioner = preconditioner; |
|
|
|
} |
|
|
|
@ -221,7 +219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// Sets the <see cref="IIterator{T}"/> that will be used to track the iterative process.
|
|
|
|
/// </summary>
|
|
|
|
/// <param name="iterator">The iterator.</param>
|
|
|
|
public void SetIterator(IIterator<Complex32> iterator) |
|
|
|
public void SetIterator(IIterator<Numerics.Complex32> iterator) |
|
|
|
{ |
|
|
|
_iterator = iterator; |
|
|
|
} |
|
|
|
@ -230,7 +228,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// Gets or sets a series of orthonormal vectors which will be used as basis for the
|
|
|
|
/// Krylov sub-space.
|
|
|
|
/// </summary>
|
|
|
|
public IList<Vector<Complex32>> StartingVectors |
|
|
|
public IList<Vector<Numerics.Complex32>> StartingVectors |
|
|
|
{ |
|
|
|
[DebuggerStepThrough] |
|
|
|
get { return _startingVectors; } |
|
|
|
@ -276,7 +274,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
|
|
|
|
/// <param name="vector">The solution vector, <c>b</c>.</param>
|
|
|
|
/// <returns>The result vector, <c>x</c>.</returns>
|
|
|
|
public Vector<Complex32> Solve(Matrix<Complex32> matrix, Vector<Complex32> vector) |
|
|
|
public Vector<Numerics.Complex32> Solve(Matrix<Numerics.Complex32> matrix, Vector<Numerics.Complex32> vector) |
|
|
|
{ |
|
|
|
if (vector == null) |
|
|
|
{ |
|
|
|
@ -295,7 +293,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
|
|
|
|
/// <param name="input">The solution vector, <c>b</c></param>
|
|
|
|
/// <param name="result">The result vector, <c>x</c></param>
|
|
|
|
public void Solve(Matrix<Complex32> matrix, Vector<Complex32> input, Vector<Complex32> result) |
|
|
|
public void Solve(Matrix<Numerics.Complex32> matrix, Vector<Numerics.Complex32> input, Vector<Numerics.Complex32> result) |
|
|
|
{ |
|
|
|
// If we were stopped before, we are no longer
|
|
|
|
// We're doing this at the start of the method to ensure
|
|
|
|
@ -382,7 +380,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
CalculateTrueResidual(matrix, residuals, xtemp, input); |
|
|
|
|
|
|
|
// Define the temporary values
|
|
|
|
var c = new Complex32[k]; |
|
|
|
var c = new Numerics.Complex32[k]; |
|
|
|
|
|
|
|
// Define the temporary vectors
|
|
|
|
var gtemp = new DenseVector(residuals.Count); |
|
|
|
@ -492,7 +490,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
zw.Clear(); |
|
|
|
|
|
|
|
// FOR (s = i, ...., k-1) AND j >= 1
|
|
|
|
Complex32 beta; |
|
|
|
Numerics.Complex32 beta; |
|
|
|
if (iterationNumber >= 1) |
|
|
|
{ |
|
|
|
for (var s = i; s < k - 1; s++) |
|
|
|
@ -639,7 +637,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// the <paramref name="numberOfVariables"/> is smaller than
|
|
|
|
/// the <paramref name="maximumNumberOfStartingVectors"/>.
|
|
|
|
/// </returns>
|
|
|
|
static IList<Vector<Complex32>> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) |
|
|
|
static IList<Vector<Numerics.Complex32>> 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
|
|
|
|
@ -653,12 +651,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
var matrix = new DenseMatrix(numberOfVariables, count); |
|
|
|
for (var i = 0; i < matrix.ColumnCount; i++) |
|
|
|
{ |
|
|
|
var samples = new Complex32[matrix.RowCount]; |
|
|
|
var samples = new Numerics.Complex32[matrix.RowCount]; |
|
|
|
var samplesRe = distribution.Samples().Take(matrix.RowCount).ToArray(); |
|
|
|
var samplesIm = distribution.Samples().Take(matrix.RowCount).ToArray(); |
|
|
|
for (int j = 0; j < matrix.RowCount; j++) |
|
|
|
{ |
|
|
|
samples[j] = new Complex32((float) samplesRe[j], (float) samplesIm[j]); |
|
|
|
samples[j] = new Numerics.Complex32((float) samplesRe[j], (float) samplesIm[j]); |
|
|
|
} |
|
|
|
|
|
|
|
// Set the column
|
|
|
|
@ -670,7 +668,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
var orthogonalMatrix = gs.Q; |
|
|
|
|
|
|
|
// Now transfer this to vectors
|
|
|
|
var result = new List<Vector<Complex32>>(); |
|
|
|
var result = new List<Vector<Numerics.Complex32>>(); |
|
|
|
for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) |
|
|
|
{ |
|
|
|
result.Add(orthogonalMatrix.Column(i)); |
|
|
|
@ -688,9 +686,9 @@ 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>
|
|
|
|
static Vector<Complex32>[] CreateVectorArray(int arraySize, int vectorSize) |
|
|
|
static Vector<Numerics.Complex32>[] CreateVectorArray(int arraySize, int vectorSize) |
|
|
|
{ |
|
|
|
var result = new Vector<Complex32>[arraySize]; |
|
|
|
var result = new Vector<Numerics.Complex32>[arraySize]; |
|
|
|
for (var i = 0; i < result.Length; i++) |
|
|
|
{ |
|
|
|
result[i] = new DenseVector(vectorSize); |
|
|
|
@ -706,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>
|
|
|
|
static void CalculateTrueResidual(Matrix<Complex32> matrix, Vector<Complex32> residual, Vector<Complex32> x, Vector<Complex32> b) |
|
|
|
static void CalculateTrueResidual(Matrix<Numerics.Complex32> matrix, Vector<Numerics.Complex32> residual, Vector<Numerics.Complex32> x, Vector<Numerics.Complex32> b) |
|
|
|
{ |
|
|
|
// -Ax = residual
|
|
|
|
matrix.Multiply(x, residual); |
|
|
|
@ -724,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>
|
|
|
|
bool ShouldContinue(int iterationNumber, Vector<Complex32> result, Vector<Complex32> source, Vector<Complex32> residuals) |
|
|
|
bool ShouldContinue(int iterationNumber, Vector<Numerics.Complex32> result, Vector<Numerics.Complex32> source, Vector<Numerics.Complex32> residuals) |
|
|
|
{ |
|
|
|
if (_hasBeenStopped) |
|
|
|
{ |
|
|
|
@ -748,7 +746,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
|
|
|
|
/// <param name="input">The solution matrix, <c>B</c>.</param>
|
|
|
|
/// <returns>The result matrix, <c>X</c>.</returns>
|
|
|
|
public Matrix<Complex32> Solve(Matrix<Complex32> matrix, Matrix<Complex32> input) |
|
|
|
public Matrix<Numerics.Complex32> Solve(Matrix<Numerics.Complex32> matrix, Matrix<Numerics.Complex32> input) |
|
|
|
{ |
|
|
|
if (matrix == null) |
|
|
|
{ |
|
|
|
@ -772,7 +770,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative |
|
|
|
/// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
|
|
|
|
/// <param name="input">The solution matrix, <c>B</c>.</param>
|
|
|
|
/// <param name="result">The result matrix, <c>X</c></param>
|
|
|
|
public void Solve(Matrix<Complex32> matrix, Matrix<Complex32> input, Matrix<Complex32> result) |
|
|
|
public void Solve(Matrix<Numerics.Complex32> matrix, Matrix<Numerics.Complex32> input, Matrix<Numerics.Complex32> result) |
|
|
|
{ |
|
|
|
if (matrix == null) |
|
|
|
{ |
Christoph Ruegg
13 years ago