diff --git a/src/Examples/LinearAlgebra/IterativeSolvers/CompositeSolverExample.cs b/src/Examples/LinearAlgebra/IterativeSolvers/CompositeSolverExample.cs index aaf68536..fdd1e23f 100644 --- a/src/Examples/LinearAlgebra/IterativeSolvers/CompositeSolverExample.cs +++ b/src/Examples/LinearAlgebra/IterativeSolvers/CompositeSolverExample.cs @@ -109,10 +109,7 @@ namespace Examples.LinearAlgebra.IterativeSolversExamples // Load all suitable solvers from current assembly. Below in this example, there is user-defined solver // "class UserBiCgStab : IIterativeSolverSetup" which uses regular BiCgStab solver. But user may create any other solver // and solver setup classes which implement IIterativeSolverSetup and pass assembly to next function: - CompositeSolver.LoadSolverInformationFromAssembly(Assembly.GetExecutingAssembly()); - - // Create composite solver - var solver = new CompositeSolver(); + var solver = new CompositeSolver(SolverSetup.LoadFromAssembly(Assembly.GetExecutingAssembly())); // 1. Solve the matrix equation var resultX = solver.Solve(matrixA, vectorB, monitor); @@ -178,11 +175,16 @@ namespace Examples.LinearAlgebra.IterativeSolversExamples /// given by this setup. /// /// A new . - public IIterativeSolver CreateNew() + public IIterativeSolver CreateSolver() { return new BiCgStab(); } + public IPreConditioner CreatePreconditioner() + { + return null; + } + /// /// Gets the relative speed of the solver. /// diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/BiCgStab.cs index 95dda13f..fc5d2f58 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/BiCgStab.cs @@ -72,12 +72,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// public sealed class BiCgStab : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates if the user has stopped the solver. /// @@ -94,36 +88,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public BiCgStab(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -142,10 +106,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient , A. /// The solution , b. /// The result , x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -156,7 +120,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient , A. /// The solution , b. /// The result , x. - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -201,12 +165,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // Compute r_0 = b - Ax_0 for some initial guess x_0 // In this case we take x_0 = vector @@ -270,7 +234,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } // SOLVE Mp~ = p_i // M = preconditioner - _preconditioner.Approximate(vecP, vecPdash); + preconditioner.Approximate(vecP, vecPdash); // nu_i = Ap~ matrix.Multiply(vecPdash, nu); @@ -314,7 +278,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } // SOLVE Ms~ = s - _preconditioner.Approximate(vecS, vecSdash); + preconditioner.Approximate(vecS, vecSdash); // temp = As~ matrix.Multiply(vecSdash, temp); @@ -407,10 +371,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient , A. /// The solution , B. /// The result , X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -421,7 +385,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient , A. /// The solution , B. /// The result , X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -433,9 +397,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs index b04a03fb..598e1048 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs @@ -30,9 +30,7 @@ using System; using System.Collections.Generic; -using System.IO; using System.Linq; -using System.Reflection; using MathNet.Numerics.LinearAlgebra.Solvers; using MathNet.Numerics.Properties; @@ -43,7 +41,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; - #endif /// @@ -63,276 +60,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// public sealed class CompositeSolver : IIterativeSolver { - #region Internal class - DoubleComparer - - /// - /// An IComparer used to compare double precision floating points. - /// - /// NOTE: The instance of this class is used only in . If C# suppports interface inheritence - /// NOTE: and methods in anonymous types, then this class should be deleted and anonymous type implemented with IComaprer support - /// NOTE: in constructor - public sealed class DoubleComparer : IComparer - { - /// - /// Compares two double values based on the selected comparison method. - /// - /// The first double to compare. - /// The second double to compare. - /// - /// A 32-bit signed integer that indicates the relative order of the objects being compared. The return - /// value has the following meanings: - /// Value Meaning Less than zero This object is less than the other parameter. - /// Zero This object is equal to other. - /// Greater than zero This object is greater than other. - /// - public int Compare(double x, double y) - { - return x.CompareTo(y, 1); - } - } - - #endregion - -#if PORTABLE - private static readonly Dictionary>> SolverSetups = new Dictionary>>(); -#else - /// - /// The collection of iterative solver setups. Stored based on the - /// ratio between the relative speed and relative accuracy. - /// - static readonly SortedList>> SolverSetups = new SortedList>>(new DoubleComparer()); -#endif - - #region Solver information loading methods - - /// - /// Loads all the available objects from the MathNet.Numerics assembly. - /// - public static void LoadSolverInformation() - { - LoadSolverInformation(new Type[0]); - } - - /// - /// Loads the available objects from the MathNet.Numerics assembly. - /// - /// The types that should not be loaded. - public static void LoadSolverInformation(Type[] typesToExclude) - { - LoadSolverInformationFromAssembly(Assembly.GetExecutingAssembly(), typesToExclude); - } - -#if !PORTABLE - /// - /// Loads the available objects from the assembly specified by the file location. - /// - /// The fully qualified path to the assembly. - public static void LoadSolverInformationFromAssembly(string assemblyLocation) - { - LoadSolverInformationFromAssembly(assemblyLocation, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the file location. - /// - /// The fully qualified path to the assembly. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(string assemblyLocation, params Type[] typesToExclude) - { - if (assemblyLocation == null) - { - throw new ArgumentNullException("assemblyLocation"); - } - - if (assemblyLocation.Length == 0) - { - throw new ArgumentException(); - } - - if (!File.Exists(assemblyLocation)) - { - throw new FileNotFoundException(); - } - - // Get the assembly name - var assemblyFileName = Path.GetFileNameWithoutExtension(assemblyLocation); - - // Now load the assembly with an AssemblyName - var assemblyName = new AssemblyName(assemblyFileName); - var assembly = Assembly.Load(assemblyName.FullName); - - // Can't get this because we checked that the file exists. - // FileNotFoundException --> Can't get this because we checked that the file exists. - // FileLoadException - // BadImageFormatException - // Now we can load the solver information. - LoadSolverInformationFromAssembly(assembly, typesToExclude); - } -#endif - - /// - /// Loads the available objects from the assembly specified by the assembly name. - /// - /// The of the assembly that should be searched for setup objects. - public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName) - { - LoadSolverInformationFromAssembly(assemblyName, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the assembly name. - /// - /// The of the assembly that should be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName, params Type[] typesToExclude) - { - if (assemblyName == null) - { - throw new ArgumentNullException("assemblyName"); - } - - var assembly = Assembly.Load(assemblyName.FullName); - - // May throw: - // ArgumentNullException --> Can't get this because we checked it. - // FileNotFoundException - // FileLoadException - // BadImageFormatException - // Now we can load the solver information. - LoadSolverInformationFromAssembly(assembly, typesToExclude); - } - - /// - /// Loads the available objects from the assembly specified by the type. - /// - /// The type in the assembly which should be searched for setup objects. - public static void LoadSolverInformationFromAssembly(Type typeInAssembly) - { - LoadSolverInformationFromAssembly(typeInAssembly, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the type. - /// - /// The type in the assembly which should be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(Type typeInAssembly, params Type[] typesToExclude) - { - if (typeInAssembly == null) - { - throw new ArgumentNullException("typeInAssembly"); - } - - LoadSolverInformationFromAssembly(typeInAssembly.Assembly, typesToExclude); - } - - /// - /// Loads the available objects from the specified assembly. - /// - /// The assembly which will be searched for setup objects. - public static void LoadSolverInformationFromAssembly(Assembly assembly) - { - LoadSolverInformationFromAssembly(assembly, new Type[0]); - } - - /// - /// Loads the available objects from the specified assembly. - /// - /// The assembly which will be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(Assembly assembly, params Type[] typesToExclude) - { - if (assembly == null) - { - throw new ArgumentNullException("assembly"); - } - - if (typesToExclude == null) - { - throw new ArgumentNullException("typesToExclude"); - } - - var excludedTypes = new List(typesToExclude); - - // Load all the types in the assembly - // Find all the types that implement IIterativeSolverSetup - // Create an object of each of these types - // Get the type of the iterative solver that will be instantiated by the setup object - // Check if it's on the excluding list, if so throw the setup object away otherwise keep it. - var interfaceTypes = new List(); - foreach (var type in assembly.GetTypes().Where(type => (!type.IsAbstract && !type.IsEnum && !type.IsInterface && type.IsVisible))) - { - interfaceTypes.Clear(); - interfaceTypes.AddRange(type.GetInterfaces()); - if (!interfaceTypes.Any(match => typeof (IIterativeSolverSetup).IsAssignableFrom(match))) - { - continue; - } - - // See if we actually want this type of iterative solver - IIterativeSolverSetup setup; - try - { - // If something goes wrong we just ignore it and move on with the next type. - // There should probably be a log somewhere indicating that something went wrong? - setup = (IIterativeSolverSetup) Activator.CreateInstance(type); - } - catch (ArgumentException) - { - continue; - } - catch (NotSupportedException) - { - continue; - } - catch (TargetInvocationException) - { - continue; - } -#if !PORTABLE - catch (MethodAccessException) - { - continue; - } - catch (MissingMethodException) - { - continue; - } -#endif - catch (MemberAccessException) - { - continue; - } - catch (TypeLoadException) - { - continue; - } - - if (excludedTypes.Any(match => match.IsAssignableFrom(setup.SolverType) || - match.IsAssignableFrom(setup.PreconditionerType))) - { - continue; - } - - // Ok we want the solver, so store the object - var ratio = setup.SolutionSpeed/setup.Reliability; - if (!SolverSetups.ContainsKey(ratio)) - { - SolverSetups.Add(ratio, new List>()); - } - - var list = SolverSetups[ratio]; - list.Add(setup); - } - } - - #endregion - /// - /// The collection of solvers that will be used to + /// The collection of solvers that will be used /// - readonly List> _solvers = new List>(); + readonly List, IPreConditioner>> _solvers; /// /// The status of the calculation. @@ -350,11 +81,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// IIterativeSolver _currentSolver; - /// - /// Initializes a new instance of the class with the default iterator. - /// - public CompositeSolver() + public CompositeSolver(IEnumerable> solvers) { + _solvers = solvers.Select(setup => new Tuple, IPreConditioner>(setup.CreateSolver(), setup.CreatePreconditioner() ?? new UnitPreconditioner())).ToList(); } /// @@ -380,10 +109,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -394,7 +123,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -419,13 +148,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - // Load the solvers into our own internal data structure - // Once we have solvers we can always reuse them. - if (_solvers.Count == 0) - { - LoadSolvers(); - } - // Create a copy of the solution and result vectors so we can use them // later on var internalInput = input.Clone(); @@ -434,7 +156,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped)) { // Store a reference to the solver so we can stop it. - _currentSolver = solver; + _currentSolver = solver.Item1; try { @@ -442,7 +164,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator.Reset(); // Start the solver - solver.Solve(matrix, internalInput, internalResult, iterator); + solver.Item1.Solve(matrix, internalInput, internalResult, iterator, solver.Item2); _status = iterator.Status; } catch (Exception) @@ -490,26 +212,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers _currentSolver = null; } - /// - /// Load solvers - /// - void LoadSolvers() - { - if (SolverSetups.Count == 0) - { - throw new Exception("IIterativeSolverSetup objects not found"); - } - -#if PORTABLE - foreach (var setup in SolverSetups.OrderBy(solver => solver.Key, new DoubleComparer()).Select(pair => pair.Value).SelectMany(setups => setups)) -#else - foreach (var setup in SolverSetups.Select(pair => pair.Value).SelectMany(setups => setups)) -#endif - { - _solvers.Add(setup.CreateNew()); - } - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -517,10 +219,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -531,7 +233,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -543,9 +245,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs index 6ef79844..04006d84 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs @@ -70,12 +70,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// public sealed class GpBiCg : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to null, in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates the number of BiCGStab steps should be taken /// before switching. @@ -104,27 +98,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public GpBiCg(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets the number of steps taken with the BiCgStab algorithm /// before switching over to the GPBiCG algorithm. @@ -163,15 +136,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } } - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -191,10 +155,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -205,7 +169,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -245,12 +209,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // x_0 is initial guess // Take x_0 = 0 @@ -298,7 +262,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers residuals.Add(temp2, p); // Solve M b_k = p_k - _preconditioner.Approximate(p, temp); + preconditioner.Approximate(p, temp); // s_k = A b_k matrix.Multiply(temp, s); @@ -322,7 +286,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers residuals.Add(temp2, t); // Solve M d_k = t_k - _preconditioner.Approximate(t, temp); + preconditioner.Approximate(t, temp); // c_k = A d_k matrix.Multiply(temp, c); @@ -427,7 +391,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers c.Add(temp2, w); // Get the real value - _preconditioner.Approximate(xtemp, result); + preconditioner.Approximate(xtemp, result); // Now check for convergence if (!ShouldContinue(iterator, iterationNumber, result, input, residuals)) @@ -507,10 +471,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -521,7 +485,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -533,9 +497,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs index 66abbb9a..527e33b3 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs @@ -32,6 +32,7 @@ using System; using System.Collections.Generic; using System.Diagnostics; using System.Linq; +using System.Threading.Tasks; using MathNet.Numerics.Distributions; using MathNet.Numerics.LinearAlgebra.Solvers; using MathNet.Numerics.Properties; @@ -76,12 +77,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// const int DefaultNumberOfStartingVectors = 50; - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// The collection of starting vectors which are used as the basis for the Krylov sub-space. /// @@ -108,27 +103,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public MlkBiCgStab(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets the number of starting vectors. /// @@ -161,15 +135,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers _numberOfStartingVectors = DefaultNumberOfStartingVectors; } - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets a series of orthonormal vectors which will be used as basis for the /// Krylov sub-space. @@ -211,10 +176,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -225,7 +190,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -265,12 +230,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // Choose an initial guess x_0 // Take x_0 = 0 @@ -335,7 +300,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers while (ShouldContinue(iterator, iterationNumber, xtemp, input, residuals)) { // SOLVE M g~_((j-1)k+k) = g_((j-1)k+k) - _preconditioner.Approximate(g[k - 1], gtemp); + preconditioner.Approximate(g[k - 1], gtemp); // w_((j-1)k+k) = A g~_((j-1)k+k) matrix.Multiply(gtemp, w[k - 1]); @@ -355,7 +320,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers residuals.Add(temp, u); // SOLVE M u~_(jk+1) = u_(jk+1) - _preconditioner.Approximate(u, temp1); + preconditioner.Approximate(u, temp1); temp1.CopyTo(utemp); // rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2 @@ -509,7 +474,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers temp2.CopyTo(u); // SOLVE M g~_(jk+i) = g_(jk+i) - _preconditioner.Approximate(g[i], gtemp); + 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); @@ -672,10 +637,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -686,7 +651,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -698,9 +663,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs index c4bd3eb5..2c863fb0 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs @@ -61,12 +61,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// public sealed class TFQMR : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates if the user has stopped the solver. /// @@ -83,36 +77,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public TFQMR(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -131,10 +95,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -145,7 +109,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -185,12 +149,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); var d = new DenseVector(input.Count); var r = DenseVector.OfVector(input); @@ -221,7 +185,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers // Calculate the initial values for v // M temp = yEven - _preconditioner.Approximate(yeven, temp); + preconditioner.Approximate(yeven, temp); // v = A temp matrix.Multiply(temp, v); @@ -253,7 +217,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers yeven.Add(temp1, yodd); // Solve M temp = yOdd - _preconditioner.Approximate(yodd, temp); + preconditioner.Approximate(yodd, temp); // uOdd = A temp matrix.Multiply(temp, uodd); @@ -293,7 +257,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers if (!ShouldContinue(iterator, iterationNumber, result, input, pseudoResiduals)) { // Calculate the real values - _preconditioner.Approximate(x, result); + preconditioner.Approximate(x, result); // Calculate the true residual. Use the temp vector for that // so that we don't pollute the pseudoResidual vector for no @@ -329,7 +293,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers pseudoResiduals.Add(temp1, yeven); // Solve M temp = yOdd - _preconditioner.Approximate(yeven, temp); + preconditioner.Approximate(yeven, temp); // uOdd = A temp matrix.Multiply(temp, ueven); @@ -343,7 +307,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } // Calculate the real values - _preconditioner.Approximate(x, result); + preconditioner.Approximate(x, result); iterationNumber++; } @@ -407,10 +371,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -421,7 +385,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -433,9 +397,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/BiCgStab.cs index 6409319e..009be4cb 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/BiCgStab.cs @@ -29,6 +29,7 @@ // using System; +using MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners; using MathNet.Numerics.LinearAlgebra.Solvers; using MathNet.Numerics.Properties; @@ -65,12 +66,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// public sealed class BiCgStab : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates if the user has stopped the solver. /// @@ -87,36 +82,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public BiCgStab(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -135,10 +100,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient , A. /// The solution , b. /// The result , x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -149,7 +114,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient , A. /// The solution , b. /// The result , x. - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -194,12 +159,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // Compute r_0 = b - Ax_0 for some initial guess x_0 // In this case we take x_0 = vector @@ -263,7 +228,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } // SOLVE Mp~ = p_i // M = preconditioner - _preconditioner.Approximate(vecP, vecPdash); + preconditioner.Approximate(vecP, vecPdash); // nu_i = Ap~ matrix.Multiply(vecPdash, nu); @@ -307,7 +272,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } // SOLVE Ms~ = s - _preconditioner.Approximate(vecS, vecSdash); + preconditioner.Approximate(vecS, vecSdash); // temp = As~ matrix.Multiply(vecSdash, temp); @@ -400,10 +365,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient , A. /// The solution , B. /// The result , X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -414,7 +379,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient , A. /// The solution , B. /// The result , X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -426,9 +391,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/CompositeSolver.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/CompositeSolver.cs index 3552c548..3f1350ec 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/CompositeSolver.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/CompositeSolver.cs @@ -30,9 +30,7 @@ using System; using System.Collections.Generic; -using System.IO; using System.Linq; -using System.Reflection; using MathNet.Numerics.LinearAlgebra.Solvers; using MathNet.Numerics.Properties; @@ -55,276 +53,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// public sealed class CompositeSolver : IIterativeSolver { - #region Internal class - DoubleComparer - - /// - /// An IComparer used to compare double precision floating points. - /// - /// NOTE: The instance of this class is used only in . If C# suppports interface inheritence - /// NOTE: and methods in anonymous types, then this class should be deleted and anonymous type implemented with IComaprer support - /// NOTE: in constructor - public sealed class DoubleComparer : IComparer - { - /// - /// Compares two double values based on the selected comparison method. - /// - /// The first double to compare. - /// The second double to compare. - /// - /// A 32-bit signed integer that indicates the relative order of the objects being compared. The return - /// value has the following meanings: - /// Value Meaning Less than zero This object is less than the other parameter. - /// Zero This object is equal to other. - /// Greater than zero This object is greater than other. - /// - public int Compare(double x, double y) - { - return x.CompareTo(y, 1); - } - } - - #endregion - -#if PORTABLE - private static readonly Dictionary>> SolverSetups = new Dictionary>>(); -#else /// - /// The collection of iterative solver setups. Stored based on the - /// ratio between the relative speed and relative accuracy. + /// The collection of solvers that will be used /// - static readonly SortedList>> SolverSetups = new SortedList>>(new DoubleComparer()); -#endif - - #region Solver information loading methods - - /// - /// Loads all the available objects from the MathNet.Numerics assembly. - /// - public static void LoadSolverInformation() - { - LoadSolverInformation(new Type[0]); - } - - /// - /// Loads the available objects from the MathNet.Numerics assembly. - /// - /// The types that should not be loaded. - public static void LoadSolverInformation(Type[] typesToExclude) - { - LoadSolverInformationFromAssembly(Assembly.GetExecutingAssembly(), typesToExclude); - } - -#if !PORTABLE - /// - /// Loads the available objects from the assembly specified by the file location. - /// - /// The fully qualified path to the assembly. - public static void LoadSolverInformationFromAssembly(string assemblyLocation) - { - LoadSolverInformationFromAssembly(assemblyLocation, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the file location. - /// - /// The fully qualified path to the assembly. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(string assemblyLocation, params Type[] typesToExclude) - { - if (assemblyLocation == null) - { - throw new ArgumentNullException("assemblyLocation"); - } - - if (assemblyLocation.Length == 0) - { - throw new ArgumentException(); - } - - if (!File.Exists(assemblyLocation)) - { - throw new FileNotFoundException(); - } - - // Get the assembly name - var assemblyFileName = Path.GetFileNameWithoutExtension(assemblyLocation); - - // Now load the assembly with an AssemblyName - var assemblyName = new AssemblyName(assemblyFileName); - var assembly = Assembly.Load(assemblyName.FullName); - - // Can't get this because we checked that the file exists. - // FileNotFoundException --> Can't get this because we checked that the file exists. - // FileLoadException - // BadImageFormatException - // Now we can load the solver information. - LoadSolverInformationFromAssembly(assembly, typesToExclude); - } -#endif - - /// - /// Loads the available objects from the assembly specified by the assembly name. - /// - /// The of the assembly that should be searched for setup objects. - public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName) - { - LoadSolverInformationFromAssembly(assemblyName, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the assembly name. - /// - /// The of the assembly that should be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName, params Type[] typesToExclude) - { - if (assemblyName == null) - { - throw new ArgumentNullException("assemblyName"); - } - - var assembly = Assembly.Load(assemblyName.FullName); - - // May throw: - // ArgumentNullException --> Can't get this because we checked it. - // FileNotFoundException - // FileLoadException - // BadImageFormatException - // Now we can load the solver information. - LoadSolverInformationFromAssembly(assembly, typesToExclude); - } - - /// - /// Loads the available objects from the assembly specified by the type. - /// - /// The type in the assembly which should be searched for setup objects. - public static void LoadSolverInformationFromAssembly(Type typeInAssembly) - { - LoadSolverInformationFromAssembly(typeInAssembly, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the type. - /// - /// The type in the assembly which should be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(Type typeInAssembly, params Type[] typesToExclude) - { - if (typeInAssembly == null) - { - throw new ArgumentNullException("typeInAssembly"); - } - - LoadSolverInformationFromAssembly(typeInAssembly.Assembly, typesToExclude); - } - - /// - /// Loads the available objects from the specified assembly. - /// - /// The assembly which will be searched for setup objects. - public static void LoadSolverInformationFromAssembly(Assembly assembly) - { - LoadSolverInformationFromAssembly(assembly, new Type[0]); - } - - /// - /// Loads the available objects from the specified assembly. - /// - /// The assembly which will be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(Assembly assembly, params Type[] typesToExclude) - { - if (assembly == null) - { - throw new ArgumentNullException("assembly"); - } - - if (typesToExclude == null) - { - throw new ArgumentNullException("typesToExclude"); - } - - var excludedTypes = new List(typesToExclude); - - // Load all the types in the assembly - // Find all the types that implement IIterativeSolverSetup - // Create an object of each of these types - // Get the type of the iterative solver that will be instantiated by the setup object - // Check if it's on the excluding list, if so throw the setup object away otherwise keep it. - var interfaceTypes = new List(); - foreach (var type in assembly.GetTypes().Where(type => (!type.IsAbstract && !type.IsEnum && !type.IsInterface && type.IsVisible))) - { - interfaceTypes.Clear(); - interfaceTypes.AddRange(type.GetInterfaces()); - if (!interfaceTypes.Any(match => typeof (IIterativeSolverSetup).IsAssignableFrom(match))) - { - continue; - } - - // See if we actually want this type of iterative solver - IIterativeSolverSetup setup; - try - { - // If something goes wrong we just ignore it and move on with the next type. - // There should probably be a log somewhere indicating that something went wrong? - setup = (IIterativeSolverSetup) Activator.CreateInstance(type); - } - catch (ArgumentException) - { - continue; - } - catch (NotSupportedException) - { - continue; - } - catch (TargetInvocationException) - { - continue; - } -#if !PORTABLE - catch (MethodAccessException) - { - continue; - } - catch (MissingMethodException) - { - continue; - } -#endif - catch (MemberAccessException) - { - continue; - } - catch (TypeLoadException) - { - continue; - } - - if (excludedTypes.Any(match => match.IsAssignableFrom(setup.SolverType) || - match.IsAssignableFrom(setup.PreconditionerType))) - { - continue; - } - - // Ok we want the solver, so store the object - var ratio = setup.SolutionSpeed/setup.Reliability; - if (!SolverSetups.ContainsKey(ratio)) - { - SolverSetups.Add(ratio, new List>()); - } - - var list = SolverSetups[ratio]; - list.Add(setup); - } - } - - #endregion - - /// - /// The collection of solvers that will be used to - /// - readonly List> _solvers = new List>(); + readonly List, IPreConditioner>> _solvers; /// /// The status of the calculation. @@ -342,11 +74,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// IIterativeSolver _currentSolver; - /// - /// Initializes a new instance of the class with the default iterator. - /// - public CompositeSolver() + public CompositeSolver(IEnumerable> solvers) { + _solvers = solvers.Select(setup => new Tuple, IPreConditioner>(setup.CreateSolver(), setup.CreatePreconditioner() ?? new UnitPreconditioner())).ToList(); } /// @@ -372,10 +102,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -386,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -411,13 +141,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - // Load the solvers into our own internal data structure - // Once we have solvers we can always reuse them. - if (_solvers.Count == 0) - { - LoadSolvers(); - } - // Create a copy of the solution and result vectors so we can use them // later on var internalInput = input.Clone(); @@ -426,7 +149,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped)) { // Store a reference to the solver so we can stop it. - _currentSolver = solver; + _currentSolver = solver.Item1; try { @@ -434,7 +157,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator.Reset(); // Start the solver - solver.Solve(matrix, internalInput, internalResult, iterator); + solver.Item1.Solve(matrix, internalInput, internalResult, iterator, solver.Item2); _status = iterator.Status; } catch (Exception) @@ -482,26 +205,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers _currentSolver = null; } - /// - /// Load solvers - /// - void LoadSolvers() - { - if (SolverSetups.Count == 0) - { - throw new Exception("IIterativeSolverSetup objects not found"); - } - -#if PORTABLE - foreach (var setup in SolverSetups.OrderBy(solver => solver.Key, new DoubleComparer()).Select(pair => pair.Value).SelectMany(setups => setups)) -#else - foreach (var setup in SolverSetups.Select(pair => pair.Value).SelectMany(setups => setups)) -#endif - { - _solvers.Add(setup.CreateNew()); - } - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -509,10 +212,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -523,7 +226,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -535,9 +238,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/GpBiCg.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/GpBiCg.cs index 8172911d..127050b3 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/GpBiCg.cs @@ -63,12 +63,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// public sealed class GpBiCg : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to null, in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates the number of BiCGStab steps should be taken /// before switching. @@ -97,27 +91,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public GpBiCg(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets the number of steps taken with the BiCgStab algorithm /// before switching over to the GPBiCG algorithm. @@ -156,15 +129,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } } - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -184,10 +148,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -198,7 +162,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -243,12 +207,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // x_0 is initial guess // Take x_0 = 0 @@ -296,7 +260,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers residuals.Add(temp2, p); // Solve M b_k = p_k - _preconditioner.Approximate(p, temp); + preconditioner.Approximate(p, temp); // s_k = A b_k matrix.Multiply(temp, s); @@ -320,7 +284,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers residuals.Add(temp2, t); // Solve M d_k = t_k - _preconditioner.Approximate(t, temp); + preconditioner.Approximate(t, temp); // c_k = A d_k matrix.Multiply(temp, c); @@ -425,7 +389,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers c.Add(temp2, w); // Get the real value - _preconditioner.Approximate(xtemp, result); + preconditioner.Approximate(xtemp, result); // Now check for convergence if (!ShouldContinue(iterator, iterationNumber, result, input, residuals)) @@ -505,10 +469,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -519,16 +483,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { throw Matrix.DimensionsDontMatch(matrix, input, result); } + if (iterator == null) + { + iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); + } + + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/MlkBiCgStab.cs index 3546ac7a..3054af9c 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/MlkBiCgStab.cs @@ -33,6 +33,7 @@ using System.Collections.Generic; using System.Diagnostics; using System.Linq; using MathNet.Numerics.Distributions; +using MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners; using MathNet.Numerics.LinearAlgebra.Solvers; using MathNet.Numerics.Properties; @@ -68,12 +69,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// const int DefaultNumberOfStartingVectors = 50; - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// The collection of starting vectors which are used as the basis for the Krylov sub-space. /// @@ -100,27 +95,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public MlkBiCgStab(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets the number of starting vectors. /// @@ -153,15 +127,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers _numberOfStartingVectors = DefaultNumberOfStartingVectors; } - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets a series of orthonormal vectors which will be used as basis for the /// Krylov sub-space. @@ -203,10 +168,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -217,7 +182,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -262,12 +227,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // Choose an initial guess x_0 // Take x_0 = 0 @@ -332,7 +297,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers while (ShouldContinue(iterator, iterationNumber, xtemp, input, residuals)) { // SOLVE M g~_((j-1)k+k) = g_((j-1)k+k) - _preconditioner.Approximate(g[k - 1], gtemp); + preconditioner.Approximate(g[k - 1], gtemp); // w_((j-1)k+k) = A g~_((j-1)k+k) matrix.Multiply(gtemp, w[k - 1]); @@ -352,7 +317,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers residuals.Add(temp, u); // SOLVE M u~_(jk+1) = u_(jk+1) - _preconditioner.Approximate(u, temp1); + preconditioner.Approximate(u, temp1); temp1.CopyTo(utemp); // rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2 @@ -506,7 +471,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers temp2.CopyTo(u); // SOLVE M g~_(jk+i) = g_(jk+i) - _preconditioner.Approximate(g[i], gtemp); + 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); @@ -669,10 +634,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -683,7 +648,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -695,9 +660,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/TFQMR.cs index 1ff8e0c5..be9806b8 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/TFQMR.cs @@ -53,12 +53,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// public sealed class TFQMR : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates if the user has stopped the solver. /// @@ -75,36 +69,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public TFQMR(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -123,10 +87,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -137,7 +101,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -182,12 +146,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); var d = new DenseVector(input.Count); var r = DenseVector.OfVector(input); @@ -218,7 +182,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers // Calculate the initial values for v // M temp = yEven - _preconditioner.Approximate(yeven, temp); + preconditioner.Approximate(yeven, temp); // v = A temp matrix.Multiply(temp, v); @@ -250,7 +214,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers yeven.Add(temp1, yodd); // Solve M temp = yOdd - _preconditioner.Approximate(yodd, temp); + preconditioner.Approximate(yodd, temp); // uOdd = A temp matrix.Multiply(temp, uodd); @@ -290,7 +254,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers if (!ShouldContinue(iterator, iterationNumber, result, input, pseudoResiduals)) { // Calculate the real values - _preconditioner.Approximate(x, result); + preconditioner.Approximate(x, result); // Calculate the true residual. Use the temp vector for that // so that we don't pollute the pseudoResidual vector for no @@ -326,7 +290,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers pseudoResiduals.Add(temp1, yeven); // Solve M temp = yOdd - _preconditioner.Approximate(yeven, temp); + preconditioner.Approximate(yeven, temp); // uOdd = A temp matrix.Multiply(temp, ueven); @@ -340,7 +304,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } // Calculate the real values - _preconditioner.Approximate(x, result); + preconditioner.Approximate(x, result); iterationNumber++; } @@ -404,10 +368,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -418,7 +382,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -430,9 +394,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/BiCgStab.cs b/src/Numerics/LinearAlgebra/Double/Solvers/BiCgStab.cs index 303129bf..3e95c79e 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/BiCgStab.cs @@ -65,12 +65,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// public sealed class BiCgStab : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates if the user has stopped the solver. /// @@ -87,36 +81,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public BiCgStab(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -135,10 +99,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient , A. /// The solution , b. /// The result , x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -149,7 +113,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient , A. /// The solution , b. /// The result , x. - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -194,12 +158,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // Compute r_0 = b - Ax_0 for some initial guess x_0 // In this case we take x_0 = vector @@ -263,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } // SOLVE Mp~ = p_i // M = preconditioner - _preconditioner.Approximate(vecP, vecPdash); + preconditioner.Approximate(vecP, vecPdash); // nu_i = Ap~ matrix.Multiply(vecPdash, nu); @@ -307,7 +271,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } // SOLVE Ms~ = s - _preconditioner.Approximate(vecS, vecSdash); + preconditioner.Approximate(vecS, vecSdash); // temp = As~ matrix.Multiply(vecSdash, temp); @@ -400,10 +364,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient , A. /// The solution , B. /// The result , X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -414,7 +378,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient , A. /// The solution , B. /// The result , X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -426,9 +390,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/CompositeSolver.cs b/src/Numerics/LinearAlgebra/Double/Solvers/CompositeSolver.cs index 47712d33..675c1bcb 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/CompositeSolver.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/CompositeSolver.cs @@ -30,9 +30,7 @@ using System; using System.Collections.Generic; -using System.IO; using System.Linq; -using System.Reflection; using MathNet.Numerics.LinearAlgebra.Solvers; using MathNet.Numerics.Properties; @@ -55,273 +53,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// public sealed class CompositeSolver : IIterativeSolver { - #region Internal class - DoubleComparer - - /// - /// An IComparer used to compare double precision floating points. - /// - public sealed class DoubleComparer : IComparer - { - /// - /// Compares two double values based on the selected comparison method. - /// - /// The first double to compare. - /// The second double to compare. - /// - /// A 32-bit signed integer that indicates the relative order of the objects being compared. The return - /// value has the following meanings: - /// Value Meaning Less than zero This object is less than the other parameter. - /// Zero This object is equal to other. - /// Greater than zero This object is greater than other. - /// - public int Compare(double x, double y) - { - return x.CompareTo(y, 1); - } - } - - #endregion - -#if PORTABLE - private static readonly Dictionary>> SolverSetups = new Dictionary>>(); -#else - /// - /// The collection of iterative solver setups. Stored based on the - /// ratio between the relative speed and relative accuracy. - /// - static readonly SortedList>> SolverSetups = new SortedList>>(new DoubleComparer()); -#endif - - #region Solver information loading methods - - /// - /// Loads all the available objects from the MathNet.Numerics assembly. - /// - public static void LoadSolverInformation() - { - LoadSolverInformation(new Type[0]); - } - - /// - /// Loads the available objects from the MathNet.Numerics assembly. - /// - /// The types that should not be loaded. - public static void LoadSolverInformation(Type[] typesToExclude) - { - LoadSolverInformationFromAssembly(Assembly.GetExecutingAssembly(), typesToExclude); - } - -#if !PORTABLE - /// - /// Loads the available objects from the assembly specified by the file location. - /// - /// The fully qualified path to the assembly. - public static void LoadSolverInformationFromAssembly(string assemblyLocation) - { - LoadSolverInformationFromAssembly(assemblyLocation, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the file location. - /// - /// The fully qualified path to the assembly. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(string assemblyLocation, params Type[] typesToExclude) - { - if (assemblyLocation == null) - { - throw new ArgumentNullException("assemblyLocation"); - } - - if (assemblyLocation.Length == 0) - { - throw new ArgumentException(); - } - - if (!File.Exists(assemblyLocation)) - { - throw new FileNotFoundException(); - } - - // Get the assembly name - var assemblyFileName = Path.GetFileNameWithoutExtension(assemblyLocation); - - // Now load the assembly with an AssemblyName - var assemblyName = new AssemblyName(assemblyFileName); - var assembly = Assembly.Load(assemblyName.FullName); - - // Can't get this because we checked that the file exists. - // FileNotFoundException --> Can't get this because we checked that the file exists. - // FileLoadException - // BadImageFormatException - // Now we can load the solver information. - LoadSolverInformationFromAssembly(assembly, typesToExclude); - } -#endif - - /// - /// Loads the available objects from the assembly specified by the assembly name. - /// - /// The of the assembly that should be searched for setup objects. - public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName) - { - LoadSolverInformationFromAssembly(assemblyName, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the assembly name. - /// - /// The of the assembly that should be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName, params Type[] typesToExclude) - { - if (assemblyName == null) - { - throw new ArgumentNullException("assemblyName"); - } - - var assembly = Assembly.Load(assemblyName.FullName); - - // May throw: - // ArgumentNullException --> Can't get this because we checked it. - // FileNotFoundException - // FileLoadException - // BadImageFormatException - // Now we can load the solver information. - LoadSolverInformationFromAssembly(assembly, typesToExclude); - } - - /// - /// Loads the available objects from the assembly specified by the type. - /// - /// The type in the assembly which should be searched for setup objects. - public static void LoadSolverInformationFromAssembly(Type typeInAssembly) - { - LoadSolverInformationFromAssembly(typeInAssembly, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the type. - /// - /// The type in the assembly which should be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(Type typeInAssembly, params Type[] typesToExclude) - { - if (typeInAssembly == null) - { - throw new ArgumentNullException("typeInAssembly"); - } - - LoadSolverInformationFromAssembly(typeInAssembly.Assembly, typesToExclude); - } - - /// - /// Loads the available objects from the specified assembly. - /// - /// The assembly which will be searched for setup objects. - public static void LoadSolverInformationFromAssembly(Assembly assembly) - { - LoadSolverInformationFromAssembly(assembly, new Type[0]); - } - /// - /// Loads the available objects from the specified assembly. + /// The collection of solvers that will be used /// - /// The assembly which will be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(Assembly assembly, params Type[] typesToExclude) - { - if (assembly == null) - { - throw new ArgumentNullException("assembly"); - } - - if (typesToExclude == null) - { - throw new ArgumentNullException("typesToExclude"); - } - - var excludedTypes = new List(typesToExclude); - - // Load all the types in the assembly - // Find all the types that implement IIterativeSolverSetup - // Create an object of each of these types - // Get the type of the iterative solver that will be instantiated by the setup object - // Check if it's on the excluding list, if so throw the setup object away otherwise keep it. - var interfaceTypes = new List(); - foreach (var type in assembly.GetTypes().Where(type => (!type.IsAbstract && !type.IsEnum && !type.IsInterface && type.IsVisible))) - { - interfaceTypes.Clear(); - interfaceTypes.AddRange(type.GetInterfaces()); - if (!interfaceTypes.Any(match => typeof (IIterativeSolverSetup).IsAssignableFrom(match))) - { - continue; - } - - // See if we actually want this type of iterative solver - IIterativeSolverSetup setup; - try - { - // If something goes wrong we just ignore it and move on with the next type. - // There should probably be a log somewhere indicating that something went wrong? - setup = (IIterativeSolverSetup) Activator.CreateInstance(type); - } - catch (ArgumentException) - { - continue; - } - catch (NotSupportedException) - { - continue; - } - catch (TargetInvocationException) - { - continue; - } -#if !PORTABLE - catch (MethodAccessException) - { - continue; - } - catch (MissingMethodException) - { - continue; - } -#endif - catch (MemberAccessException) - { - continue; - } - catch (TypeLoadException) - { - continue; - } - - if (excludedTypes.Any(match => match.IsAssignableFrom(setup.SolverType) || - match.IsAssignableFrom(setup.PreconditionerType))) - { - continue; - } - - // Ok we want the solver, so store the object - var ratio = setup.SolutionSpeed/setup.Reliability; - if (!SolverSetups.ContainsKey(ratio)) - { - SolverSetups.Add(ratio, new List>()); - } - - var list = SolverSetups[ratio]; - list.Add(setup); - } - } - - #endregion - - /// - /// The collection of solvers that will be used to - /// - readonly List> _solvers = new List>(); + readonly List, IPreConditioner>> _solvers; /// /// The status of the calculation. @@ -339,22 +74,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// IIterativeSolver _currentSolver; - /// - /// Initializes a new instance of the class with the default iterator. - /// - public CompositeSolver() - { - } - - /// - /// Gets the status of the iteration once the calculation is finished. - /// - public IterationStatus IterationResult + public CompositeSolver(IEnumerable> solvers) { - get - { - return _status; - } + _solvers = solvers.Select(setup => new Tuple, IPreConditioner>(setup.CreateSolver(), setup.CreatePreconditioner() ?? new UnitPreconditioner())).ToList(); } /// @@ -380,10 +102,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -394,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -419,13 +141,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - // Load the solvers into our own internal data structure - // Once we have solvers we can always reuse them. - if (_solvers.Count == 0) - { - LoadSolvers(); - } - // Create a copy of the solution and result vectors so we can use them // later on var internalInput = input.Clone(); @@ -434,7 +149,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped)) { // Store a reference to the solver so we can stop it. - _currentSolver = solver; + _currentSolver = solver.Item1; try { @@ -442,7 +157,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator.Reset(); // Start the solver - solver.Solve(matrix, internalInput, internalResult, iterator); + solver.Item1.Solve(matrix, internalInput, internalResult, iterator, solver.Item2); _status = iterator.Status; } catch (Exception) @@ -490,26 +205,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers _currentSolver = null; } - /// - /// Load solvers - /// - void LoadSolvers() - { - if (SolverSetups.Count == 0) - { - throw new Exception("IIterativeSolverSetup objects not found"); - } - -#if PORTABLE - foreach (var setup in SolverSetups.OrderBy(solver => solver.Key, new DoubleComparer()).Select(pair => pair.Value).SelectMany(setups => setups)) -#else - foreach (var setup in SolverSetups.Select(pair => pair.Value).SelectMany(setups => setups)) -#endif - { - _solvers.Add(setup.CreateNew()); - } - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -517,10 +212,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -531,7 +226,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -543,9 +238,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/GpBiCg.cs b/src/Numerics/LinearAlgebra/Double/Solvers/GpBiCg.cs index 667d532b..7c5b4104 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/GpBiCg.cs @@ -63,12 +63,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// public sealed class GpBiCg : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to null, in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates the number of BiCGStab steps should be taken /// before switching. @@ -97,27 +91,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public GpBiCg(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets the number of steps taken with the BiCgStab algorithm /// before switching over to the GPBiCG algorithm. @@ -162,15 +135,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } } - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -190,10 +154,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -204,7 +168,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -249,12 +213,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // x_0 is initial guess // Take x_0 = 0 @@ -302,7 +266,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers residuals.Add(temp2, p); // Solve M b_k = p_k - _preconditioner.Approximate(p, temp); + preconditioner.Approximate(p, temp); // s_k = A b_k matrix.Multiply(temp, s); @@ -326,7 +290,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers residuals.Add(temp2, t); // Solve M d_k = t_k - _preconditioner.Approximate(t, temp); + preconditioner.Approximate(t, temp); // c_k = A d_k matrix.Multiply(temp, c); @@ -431,7 +395,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers c.Add(temp2, w); // Get the real value - _preconditioner.Approximate(xtemp, result); + preconditioner.Approximate(xtemp, result); // Now check for convergence if (!ShouldContinue(iterator, iterationNumber, result, input, residuals)) @@ -511,10 +475,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -525,7 +489,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -537,9 +501,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Double/Solvers/MlkBiCgStab.cs index 42b9502b..8dce4e6d 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/MlkBiCgStab.cs @@ -68,12 +68,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// const int DefaultNumberOfStartingVectors = 50; - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// The collection of starting vectors which are used as the basis for the Krylov sub-space. /// @@ -100,27 +94,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public MlkBiCgStab(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets the number of starting vectors. /// @@ -156,15 +129,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers _numberOfStartingVectors = DefaultNumberOfStartingVectors; } - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets a series of orthonormal vectors which will be used as basis for the /// Krylov sub-space. @@ -209,10 +173,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -223,7 +187,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -268,12 +232,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // Choose an initial guess x_0 // Take x_0 = 0 @@ -338,7 +302,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers while (ShouldContinue(iterator, iterationNumber, xtemp, input, residuals)) { // SOLVE M g~_((j-1)k+k) = g_((j-1)k+k) - _preconditioner.Approximate(g[k - 1], gtemp); + preconditioner.Approximate(g[k - 1], gtemp); // w_((j-1)k+k) = A g~_((j-1)k+k) matrix.Multiply(gtemp, w[k - 1]); @@ -358,7 +322,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers residuals.Add(temp, u); // SOLVE M u~_(jk+1) = u_(jk+1) - _preconditioner.Approximate(u, temp1); + preconditioner.Approximate(u, temp1); temp1.CopyTo(utemp); // rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2 @@ -512,7 +476,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers temp2.CopyTo(u); // SOLVE M g~_(jk+i) = g_(jk+i) - _preconditioner.Approximate(g[i], gtemp); + 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); @@ -669,10 +633,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -683,7 +647,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -695,9 +659,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/TFQMR.cs b/src/Numerics/LinearAlgebra/Double/Solvers/TFQMR.cs index a6f2f372..82608e91 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/TFQMR.cs @@ -53,12 +53,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// public sealed class TFQMR : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates if the user has stopped the solver. /// @@ -75,36 +69,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public TFQMR(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -123,10 +87,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -137,7 +101,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -182,12 +146,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); var d = new DenseVector(input.Count); var r = DenseVector.OfVector(input); @@ -218,7 +182,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers // Calculate the initial values for v // M temp = yEven - _preconditioner.Approximate(yeven, temp); + preconditioner.Approximate(yeven, temp); // v = A temp matrix.Multiply(temp, v); @@ -250,7 +214,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers yeven.Add(temp1, yodd); // Solve M temp = yOdd - _preconditioner.Approximate(yodd, temp); + preconditioner.Approximate(yodd, temp); // uOdd = A temp matrix.Multiply(temp, uodd); @@ -290,7 +254,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers if (!ShouldContinue(iterator, iterationNumber, result, input, pseudoResiduals)) { // Calculate the real values - _preconditioner.Approximate(x, result); + preconditioner.Approximate(x, result); // Calculate the true residual. Use the temp vector for that // so that we don't pollute the pseudoResidual vector for no @@ -326,7 +290,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers pseudoResiduals.Add(temp1, yeven); // Solve M temp = yOdd - _preconditioner.Approximate(yeven, temp); + preconditioner.Approximate(yeven, temp); // uOdd = A temp matrix.Multiply(temp, ueven); @@ -340,7 +304,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } // Calculate the real values - _preconditioner.Approximate(x, result); + preconditioner.Approximate(x, result); iterationNumber++; } @@ -404,10 +368,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -418,7 +382,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -430,9 +394,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/BiCgStab.cs b/src/Numerics/LinearAlgebra/Single/Solvers/BiCgStab.cs index 04e5bce3..dad2eecb 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/BiCgStab.cs @@ -65,12 +65,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// public sealed class BiCgStab : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates if the user has stopped the solver. /// @@ -87,36 +81,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public BiCgStab(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -135,10 +99,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient , A. /// The solution , b. /// The result , x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -149,7 +113,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient , A. /// The solution , b. /// The result , x. - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -194,12 +158,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // Compute r_0 = b - Ax_0 for some initial guess x_0 // In this case we take x_0 = vector @@ -263,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } // SOLVE Mp~ = p_i // M = preconditioner - _preconditioner.Approximate(vecP, vecPdash); + preconditioner.Approximate(vecP, vecPdash); // nu_i = Ap~ matrix.Multiply(vecPdash, nu); @@ -307,7 +271,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } // SOLVE Ms~ = s - _preconditioner.Approximate(vecS, vecSdash); + preconditioner.Approximate(vecS, vecSdash); // temp = As~ matrix.Multiply(vecSdash, temp); @@ -400,10 +364,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient , A. /// The solution , B. /// The result , X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -414,7 +378,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient , A. /// The solution , B. /// The result , X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -426,9 +390,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/CompositeSolver.cs b/src/Numerics/LinearAlgebra/Single/Solvers/CompositeSolver.cs index 82823bfd..f06bce73 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/CompositeSolver.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/CompositeSolver.cs @@ -30,9 +30,7 @@ using System; using System.Collections.Generic; -using System.IO; using System.Linq; -using System.Reflection; using MathNet.Numerics.LinearAlgebra.Solvers; using MathNet.Numerics.Properties; @@ -55,276 +53,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// public sealed class CompositeSolver : IIterativeSolver { - #region Internal class - DoubleComparer - - /// - /// An IComparer used to compare double precision floating points. - /// - /// NOTE: The instance of this class is used only in . If C# suppports interface inheritence - /// NOTE: and methods in anonymous types, then this class should be deleted and anonymous type implemented with IComaprer support - /// NOTE: in constructor - public sealed class DoubleComparer : IComparer - { - /// - /// Compares two double values based on the selected comparison method. - /// - /// The first double to compare. - /// The second double to compare. - /// - /// A 32-bit signed integer that indicates the relative order of the objects being compared. The return - /// value has the following meanings: - /// Value Meaning Less than zero This object is less than the other parameter. - /// Zero This object is equal to other. - /// Greater than zero This object is greater than other. - /// - public int Compare(double x, double y) - { - return x.CompareTo(y, 1); - } - } - - #endregion - -#if PORTABLE - private static readonly Dictionary>> SolverSetups = new Dictionary>>(); -#else /// - /// The collection of iterative solver setups. Stored based on the - /// ratio between the relative speed and relative accuracy. + /// The collection of solvers that will be used /// - static readonly SortedList>> SolverSetups = new SortedList>>(new DoubleComparer()); -#endif - - #region Solver information loading methods - - /// - /// Loads all the available objects from the MathNet.Numerics assembly. - /// - public static void LoadSolverInformation() - { - LoadSolverInformation(new Type[0]); - } - - /// - /// Loads the available objects from the MathNet.Numerics assembly. - /// - /// The types that should not be loaded. - public static void LoadSolverInformation(Type[] typesToExclude) - { - LoadSolverInformationFromAssembly(Assembly.GetExecutingAssembly(), typesToExclude); - } - -#if !PORTABLE - /// - /// Loads the available objects from the assembly specified by the file location. - /// - /// The fully qualified path to the assembly. - public static void LoadSolverInformationFromAssembly(string assemblyLocation) - { - LoadSolverInformationFromAssembly(assemblyLocation, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the file location. - /// - /// The fully qualified path to the assembly. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(string assemblyLocation, params Type[] typesToExclude) - { - if (assemblyLocation == null) - { - throw new ArgumentNullException("assemblyLocation"); - } - - if (assemblyLocation.Length == 0) - { - throw new ArgumentException(); - } - - if (!File.Exists(assemblyLocation)) - { - throw new FileNotFoundException(); - } - - // Get the assembly name - var assemblyFileName = Path.GetFileNameWithoutExtension(assemblyLocation); - - // Now load the assembly with an AssemblyName - var assemblyName = new AssemblyName(assemblyFileName); - var assembly = Assembly.Load(assemblyName.FullName); - - // Can't get this because we checked that the file exists. - // FileNotFoundException --> Can't get this because we checked that the file exists. - // FileLoadException - // BadImageFormatException - // Now we can load the solver information. - LoadSolverInformationFromAssembly(assembly, typesToExclude); - } -#endif - - /// - /// Loads the available objects from the assembly specified by the assembly name. - /// - /// The of the assembly that should be searched for setup objects. - public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName) - { - LoadSolverInformationFromAssembly(assemblyName, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the assembly name. - /// - /// The of the assembly that should be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName, params Type[] typesToExclude) - { - if (assemblyName == null) - { - throw new ArgumentNullException("assemblyName"); - } - - var assembly = Assembly.Load(assemblyName.FullName); - - // May throw: - // ArgumentNullException --> Can't get this because we checked it. - // FileNotFoundException - // FileLoadException - // BadImageFormatException - // Now we can load the solver information. - LoadSolverInformationFromAssembly(assembly, typesToExclude); - } - - /// - /// Loads the available objects from the assembly specified by the type. - /// - /// The type in the assembly which should be searched for setup objects. - public static void LoadSolverInformationFromAssembly(Type typeInAssembly) - { - LoadSolverInformationFromAssembly(typeInAssembly, new Type[0]); - } - - /// - /// Loads the available objects from the assembly specified by the type. - /// - /// The type in the assembly which should be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(Type typeInAssembly, params Type[] typesToExclude) - { - if (typeInAssembly == null) - { - throw new ArgumentNullException("typeInAssembly"); - } - - LoadSolverInformationFromAssembly(typeInAssembly.Assembly, typesToExclude); - } - - /// - /// Loads the available objects from the specified assembly. - /// - /// The assembly which will be searched for setup objects. - public static void LoadSolverInformationFromAssembly(Assembly assembly) - { - LoadSolverInformationFromAssembly(assembly, new Type[0]); - } - - /// - /// Loads the available objects from the specified assembly. - /// - /// The assembly which will be searched for setup objects. - /// The types that should not be loaded. - public static void LoadSolverInformationFromAssembly(Assembly assembly, params Type[] typesToExclude) - { - if (assembly == null) - { - throw new ArgumentNullException("assembly"); - } - - if (typesToExclude == null) - { - throw new ArgumentNullException("typesToExclude"); - } - - var excludedTypes = new List(typesToExclude); - - // Load all the types in the assembly - // Find all the types that implement IIterativeSolverSetup - // Create an object of each of these types - // Get the type of the iterative solver that will be instantiated by the setup object - // Check if it's on the excluding list, if so throw the setup object away otherwise keep it. - var interfaceTypes = new List(); - foreach (var type in assembly.GetTypes().Where(type => (!type.IsAbstract && !type.IsEnum && !type.IsInterface && type.IsVisible))) - { - interfaceTypes.Clear(); - interfaceTypes.AddRange(type.GetInterfaces()); - if (!interfaceTypes.Any(match => typeof (IIterativeSolverSetup).IsAssignableFrom(match))) - { - continue; - } - - // See if we actually want this type of iterative solver - IIterativeSolverSetup setup; - try - { - // If something goes wrong we just ignore it and move on with the next type. - // There should probably be a log somewhere indicating that something went wrong? - setup = (IIterativeSolverSetup) Activator.CreateInstance(type); - } - catch (ArgumentException) - { - continue; - } - catch (NotSupportedException) - { - continue; - } - catch (TargetInvocationException) - { - continue; - } -#if !PORTABLE - catch (MethodAccessException) - { - continue; - } - catch (MissingMethodException) - { - continue; - } -#endif - catch (MemberAccessException) - { - continue; - } - catch (TypeLoadException) - { - continue; - } - - if (excludedTypes.Any(match => match.IsAssignableFrom(setup.SolverType) || - match.IsAssignableFrom(setup.PreconditionerType))) - { - continue; - } - - // Ok we want the solver, so store the object - var ratio = setup.SolutionSpeed/setup.Reliability; - if (!SolverSetups.ContainsKey(ratio)) - { - SolverSetups.Add(ratio, new List>()); - } - - var list = SolverSetups[ratio]; - list.Add(setup); - } - } - - #endregion - - /// - /// The collection of solvers that will be used to - /// - readonly List> _solvers = new List>(); + readonly List, IPreConditioner>> _solvers; /// /// The status of the calculation. @@ -342,15 +74,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// IIterativeSolver _currentSolver; - /// - /// Initializes a new instance of the class with the default iterator. - /// - public CompositeSolver() + public CompositeSolver(IEnumerable> solvers) { + _solvers = solvers.Select(setup => new Tuple, IPreConditioner>(setup.CreateSolver(), setup.CreatePreconditioner() ?? new UnitPreconditioner())).ToList(); } /// - /// Stops the solve process. + /// Stops the solve process. /// /// /// Note that it may take an indetermined amount of time for the solver to actually stop the process. @@ -372,10 +102,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -386,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -411,13 +141,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - // Load the solvers into our own internal data structure - // Once we have solvers we can always reuse them. - if (_solvers.Count == 0) - { - LoadSolvers(); - } - // Create a copy of the solution and result vectors so we can use them // later on var internalInput = input.Clone(); @@ -426,7 +149,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped)) { // Store a reference to the solver so we can stop it. - _currentSolver = solver; + _currentSolver = solver.Item1; try { @@ -434,7 +157,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator.Reset(); // Start the solver - solver.Solve(matrix, internalInput, internalResult, iterator); + solver.Item1.Solve(matrix, internalInput, internalResult, iterator, solver.Item2); _status = iterator.Status; } catch (Exception) @@ -482,26 +205,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers _currentSolver = null; } - /// - /// Load solvers - /// - void LoadSolvers() - { - if (SolverSetups.Count == 0) - { - throw new Exception("IIterativeSolverSetup objects not found"); - } - -#if PORTABLE - foreach (var setup in SolverSetups.OrderBy(solver => solver.Key, new DoubleComparer()).Select(pair => pair.Value).SelectMany(setups => setups)) -#else - foreach (var setup in SolverSetups.Select(pair => pair.Value).SelectMany(setups => setups)) -#endif - { - _solvers.Add(setup.CreateNew()); - } - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -509,10 +212,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -523,7 +226,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -535,9 +238,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/GpBiCg.cs b/src/Numerics/LinearAlgebra/Single/Solvers/GpBiCg.cs index 5eae6e5f..5692eab9 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/GpBiCg.cs @@ -63,12 +63,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// public sealed class GpBiCg : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to null, in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates the number of BiCGStab steps should be taken /// before switching. @@ -97,27 +91,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public GpBiCg(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets the number of steps taken with the BiCgStab algorithm /// before switching over to the GPBiCG algorithm. @@ -156,15 +129,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } } - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -184,10 +148,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -198,7 +162,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -243,12 +207,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // x_0 is initial guess // Take x_0 = 0 @@ -296,7 +260,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers residuals.Add(temp2, p); // Solve M b_k = p_k - _preconditioner.Approximate(p, temp); + preconditioner.Approximate(p, temp); // s_k = A b_k matrix.Multiply(temp, s); @@ -320,7 +284,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers residuals.Add(temp2, t); // Solve M d_k = t_k - _preconditioner.Approximate(t, temp); + preconditioner.Approximate(t, temp); // c_k = A d_k matrix.Multiply(temp, c); @@ -425,7 +389,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers c.Add(temp2, w); // Get the real value - _preconditioner.Approximate(xtemp, result); + preconditioner.Approximate(xtemp, result); // Now check for convergence if (!ShouldContinue(iterator, iterationNumber, result, input, residuals)) @@ -505,10 +469,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -519,7 +483,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -531,9 +495,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Single/Solvers/MlkBiCgStab.cs index b0494ff0..21b874d1 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/MlkBiCgStab.cs @@ -67,12 +67,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// const int DefaultNumberOfStartingVectors = 50; - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// The collection of starting vectors which are used as the basis for the Krylov sub-space. /// @@ -99,27 +93,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public MlkBiCgStab(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets the number of starting vectors. /// @@ -155,15 +128,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers _numberOfStartingVectors = DefaultNumberOfStartingVectors; } - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Gets or sets a series of orthonormal vectors which will be used as basis for the /// Krylov sub-space. @@ -208,10 +172,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -222,7 +186,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -267,12 +231,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); // Choose an initial guess x_0 // Take x_0 = 0 @@ -337,7 +301,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers while (ShouldContinue(iterator, iterationNumber, xtemp, input, residuals)) { // SOLVE M g~_((j-1)k+k) = g_((j-1)k+k) - _preconditioner.Approximate(g[k - 1], gtemp); + preconditioner.Approximate(g[k - 1], gtemp); // w_((j-1)k+k) = A g~_((j-1)k+k) matrix.Multiply(gtemp, w[k - 1]); @@ -357,7 +321,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers residuals.Add(temp, u); // SOLVE M u~_(jk+1) = u_(jk+1) - _preconditioner.Approximate(u, temp1); + preconditioner.Approximate(u, temp1); temp1.CopyTo(utemp); // rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2 @@ -511,7 +475,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers temp2.CopyTo(u); // SOLVE M g~_(jk+i) = g_(jk+i) - _preconditioner.Approximate(g[i], gtemp); + 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); @@ -672,10 +636,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -686,7 +650,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -698,9 +662,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/TFQMR.cs b/src/Numerics/LinearAlgebra/Single/Solvers/TFQMR.cs index 6408e04a..dc0db387 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/TFQMR.cs @@ -53,12 +53,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// public sealed class TFQMR : IIterativeSolver { - /// - /// The preconditioner that will be used. Can be set to , in which case the default - /// pre-conditioner will be used. - /// - IPreConditioner _preconditioner; - /// /// Indicates if the user has stopped the solver. /// @@ -75,36 +69,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers { } - /// - /// Initializes a new instance of the class. - /// - /// - /// - /// The main advantages of using a user defined are: - /// - /// It is possible to set the desired convergence limits. - /// - /// It is possible to check the reason for which the solver finished - /// the iterative procedure by calling the property. - /// - /// - /// - /// - /// The that will be used to precondition the matrix equation. - public TFQMR(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - - /// - /// Sets the that will be used to precondition the iterative process. - /// - /// The preconditioner. - public void SetPreconditioner(IPreConditioner preconditioner) - { - _preconditioner = preconditioner; - } - /// /// Stops the solve process. /// @@ -123,10 +87,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null) + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator); + Solve(matrix, vector, result, iterator, preconditioner); return result; } @@ -137,7 +101,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null) + public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure @@ -182,12 +146,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } - if (_preconditioner == null) + if (preconditioner == null) { - _preconditioner = new UnitPreconditioner(); + preconditioner = new UnitPreconditioner(); } - _preconditioner.Initialize(matrix); + preconditioner.Initialize(matrix); var d = new DenseVector(input.Count); var r = DenseVector.OfVector(input); @@ -218,7 +182,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers // Calculate the initial values for v // M temp = yEven - _preconditioner.Approximate(yeven, temp); + preconditioner.Approximate(yeven, temp); // v = A temp matrix.Multiply(temp, v); @@ -250,7 +214,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers yeven.Add(temp1, yodd); // Solve M temp = yOdd - _preconditioner.Approximate(yodd, temp); + preconditioner.Approximate(yodd, temp); // uOdd = A temp matrix.Multiply(temp, uodd); @@ -290,7 +254,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers if (!ShouldContinue(iterator, iterationNumber, result, input, pseudoResiduals)) { // Calculate the real values - _preconditioner.Approximate(x, result); + preconditioner.Approximate(x, result); // Calculate the true residual. Use the temp vector for that // so that we don't pollute the pseudoResidual vector for no @@ -326,7 +290,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers pseudoResiduals.Add(temp1, yeven); // Solve M temp = yOdd - _preconditioner.Approximate(yeven, temp); + preconditioner.Approximate(yeven, temp); // uOdd = A temp matrix.Multiply(temp, ueven); @@ -340,7 +304,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } // Calculate the real values - _preconditioner.Approximate(x, result); + preconditioner.Approximate(x, result); iterationNumber++; } @@ -404,10 +368,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null) + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null) { var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator); + Solve(matrix, input, result, iterator, preconditioner); return result; } @@ -418,7 +382,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null) + public void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null) { if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { @@ -430,9 +394,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers iterator = new Iterator(Iterator.CreateDefaultStopCriteria()); } + if (preconditioner == null) + { + preconditioner = new UnitPreconditioner(); + } + for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, input.Column(column), iterator); + var solution = Solve(matrix, input.Column(column), iterator, preconditioner); foreach (var element in solution.EnumerateNonZeroIndexed()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Solvers/IIterativeSolver.cs b/src/Numerics/LinearAlgebra/Solvers/IIterativeSolver.cs index 5bfcb97a..6f2e95a6 100644 --- a/src/Numerics/LinearAlgebra/Solvers/IIterativeSolver.cs +++ b/src/Numerics/LinearAlgebra/Solvers/IIterativeSolver.cs @@ -53,7 +53,7 @@ namespace MathNet.Numerics.LinearAlgebra.Solvers /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. - Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null); + Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreConditioner preconditioner = null); /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the @@ -62,7 +62,7 @@ namespace MathNet.Numerics.LinearAlgebra.Solvers /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x - void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null); + void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreConditioner preconditioner = null); /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the @@ -71,7 +71,7 @@ namespace MathNet.Numerics.LinearAlgebra.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. - Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null); + Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreConditioner preconditioner = null); /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the @@ -80,6 +80,6 @@ namespace MathNet.Numerics.LinearAlgebra.Solvers /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X - void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null); + void Solve(Matrix matrix, Matrix input, Matrix result, Iterator iterator = null, IPreConditioner preconditioner = null); } } diff --git a/src/Numerics/LinearAlgebra/Solvers/IIterativeSolverSetup.cs b/src/Numerics/LinearAlgebra/Solvers/IIterativeSolverSetup.cs index ad5124e6..d330a57e 100644 --- a/src/Numerics/LinearAlgebra/Solvers/IIterativeSolverSetup.cs +++ b/src/Numerics/LinearAlgebra/Solvers/IIterativeSolverSetup.cs @@ -50,11 +50,14 @@ namespace MathNet.Numerics.LinearAlgebra.Solvers Type PreconditionerType { get; } /// - /// Creates a fully functional iterative solver with the default settings - /// given by this setup. + /// Creates the iterative solver to be used. /// - /// A new . - IIterativeSolver CreateNew(); + IIterativeSolver CreateSolver(); + + /// + /// Creates the preconditioner to be used by default (can be overwritten). + /// + IPreConditioner CreatePreconditioner(); /// /// Gets the relative speed of the solver. diff --git a/src/Numerics/LinearAlgebra/Solvers/SolverSetup.cs b/src/Numerics/LinearAlgebra/Solvers/SolverSetup.cs new file mode 100644 index 00000000..caa79b46 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Solvers/SolverSetup.cs @@ -0,0 +1,138 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Reflection; + +namespace MathNet.Numerics.LinearAlgebra.Solvers +{ + public static class SolverSetup where T : struct, IEquatable, IFormattable + { + /// + /// Loads the available objects from the specified assembly. + /// + /// The assembly which will be searched for setup objects. + /// The types that should not be loaded. + public static IEnumerable> LoadFromAssembly(Assembly assembly, bool ignoreFailed, params Type[] typesToExclude) + { + var excludedTypes = new List(typesToExclude); + var setupInterfaceType = typeof (IIterativeSolverSetup); + var candidates = assembly.GetTypes() + .Where(type => !type.IsAbstract && !type.IsEnum && !type.IsInterface && type.IsVisible) + .Where(type => type.GetInterfaces().Any(setupInterfaceType.IsAssignableFrom)); + + var setups = new List>(); + foreach (var type in candidates) + { + try + { + setups.Add((IIterativeSolverSetup) Activator.CreateInstance(type)); + } + catch + { + if (!ignoreFailed) + { + throw; + } + } + } + + return setups + .Where(s => !excludedTypes.Any(t => t.IsAssignableFrom(s.SolverType) || t.IsAssignableFrom(s.PreconditionerType))) + .OrderBy(s => s.SolutionSpeed/s.Reliability); + } + + /// + /// Loads the available objects from the specified assembly. + /// + /// The assembly which will be searched for setup objects. + public static IEnumerable> LoadFromAssembly(Assembly assembly) + { + return LoadFromAssembly(assembly, true); + } + + /// + /// Loads the available objects from the specified assembly. + /// + /// The type in the assembly which should be searched for setup objects. + /// The types that should not be loaded. + public static IEnumerable> LoadFromAssembly(Type typeInAssembly, bool ignoreFailed, params Type[] typesToExclude) + { + return LoadFromAssembly(typeInAssembly.Assembly, ignoreFailed, typesToExclude); + } + + /// + /// Loads the available objects from the specified assembly. + /// + /// The type in the assembly which should be searched for setup objects. + public static IEnumerable> LoadFromAssembly(Type typeInAssembly) + { + return LoadFromAssembly(typeInAssembly.Assembly, true); + } + + /// + /// Loads the available objects from the specified assembly. + /// + /// The of the assembly that should be searched for setup objects. + /// The types that should not be loaded. + public static IEnumerable> LoadFromAssembly(AssemblyName assemblyName, bool ignoreFailed, params Type[] typesToExclude) + { + return LoadFromAssembly(Assembly.Load(assemblyName.FullName), ignoreFailed, typesToExclude); + } + + /// + /// Loads the available objects from the specified assembly. + /// + /// The of the assembly that should be searched for setup objects. + public static IEnumerable> LoadFromAssembly(AssemblyName assemblyName) + { + return LoadFromAssembly(Assembly.Load(assemblyName.FullName), true); + } + + /// + /// Loads the available objects from the Math.NET Numerics assembly. + /// + /// The types that should not be loaded. + public static IEnumerable> Load(Type[] typesToExclude) + { + return LoadFromAssembly(Assembly.GetExecutingAssembly(), false, typesToExclude); + } + + /// + /// Loads the available objects from the Math.NET Numerics assembly. + /// + public static IEnumerable> Load() + { + return LoadFromAssembly(Assembly.GetExecutingAssembly(), false); + } + } +} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 4088295d..c6999fe8 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -126,6 +126,7 @@ +