// // 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-2010 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. // namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative { using System; using System.Collections.Generic; using System.IO; using System.Linq; using System.Reflection; using Generic.Solvers.Status; using Properties; /// /// A composite matrix solver. The actual solver is made by a sequence of /// matrix solvers. /// /// /// /// Solver based on:
/// Faster PDE-based simulations using robust composite linear solvers
/// S. Bhowmicka, P. Raghavan a,*, L. McInnes b, B. Norris
/// Future Generation Computer Systems, Vol 20, 2004, pp 373–387
///
/// /// Note that if an iterator is passed to this solver it will be used for all the sub-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 /// /// The default status used if the solver is not running. /// private static readonly ICalculationStatus NonRunningStatus = new CalculationIndetermined(); /// /// The default status used if the solver is running. /// private static readonly ICalculationStatus RunningStatus = new CalculationRunning(); #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. /// private 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 /// private readonly List _solvers = new List(); /// /// The status of the calculation. /// private ICalculationStatus _status = NonRunningStatus; /// /// The iterator that is used to control the iteration process. /// private IIterator _iterator; /// /// A flag indicating if the solver has been stopped or not. /// private bool _hasBeenStopped; /// /// The solver that is currently running. Reference is used to be able to stop the /// solver if the user cancels the solve process. /// private IIterativeSolver _currentSolver; /// /// Initializes a new instance of the class with the default iterator. /// public CompositeSolver() : this(null) { } /// /// Initializes a new instance of the class with the specified iterator. /// /// The iterator that will be used to control the iteration process. public CompositeSolver(IIterator iterator) { _iterator = iterator; } /// /// Sets the IIterator that will be used to track the iterative process. /// /// The iterator. public void SetIterator(IIterator iterator) { _iterator = iterator; } /// /// Gets the status of the iteration once the calculation is finished. /// public ICalculationStatus IterationResult { get { return _status; } } /// /// Stops the solve process. /// /// /// Note that it may take an indetermined amount of time for the solver to actually stop the process. /// public void StopSolve() { _hasBeenStopped = true; if (_currentSolver != null) { _currentSolver.StopSolve(); } } /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. /// /// The coefficient matrix, A. /// The solution vector, b. /// The result vector, x. public Vector Solve(Matrix matrix, Vector vector) { if (vector == null) { throw new ArgumentNullException(); } Vector result = new DenseVector(matrix.RowCount); Solve(matrix, vector, result); return result; } /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. /// /// The coefficient matrix, A. /// The solution vector, b /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result) { // If we were stopped before, we are no longer // We're doing this at the start of the method to ensure // that we can use these fields immediately. _hasBeenStopped = false; _currentSolver = null; // Error checks if (matrix == null) { throw new ArgumentNullException("matrix"); } if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } if (input == null) { throw new ArgumentNullException("input"); } if (result == null) { throw new ArgumentNullException("result"); } if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } // Initialize the solver fields // Set the convergence monitor if (_iterator == null) { _iterator = Iterator.CreateDefault(); } // 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 = (Vector)input.Clone(); var internalResult = (Vector)result.Clone(); foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped)) { // Store a reference to the solver so we can stop it. _currentSolver = solver; try { // Reset the iterator and pass it to the solver _iterator.ResetToPrecalculationState(); solver.SetIterator(_iterator); // Start the solver solver.Solve(matrix, internalInput, internalResult); } catch (Exception) { // The solver broke down. // Log a message about this // Switch to the next preconditioner. // Reset the solution vector to the previous solution input.CopyTo(internalInput); _status = RunningStatus; continue; } // There was no fatal breakdown so check the status if (_iterator.Status is CalculationConverged) { // We're done internalResult.CopyTo(result); break; } // We're not done // Either: // - calculation finished without convergence if (_iterator.Status is CalculationStoppedWithoutConvergence) { // Copy the internal result to the result vector and // continue with the calculation. internalResult.CopyTo(result); } else { // - calculation failed --> restart with the original vector // - calculation diverged --> restart with the original vector // - Some unknown status occurred --> To be safe restart. input.CopyTo(internalInput); } } // Inside the loop we already copied the final results (if there are any) // So no need to do that again. // Clean up // No longer need the current solver _currentSolver = null; // Set the final status _status = _iterator.Status; } /// /// Load solvers /// private 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. /// /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X. public Matrix Solve(Matrix matrix, Matrix input) { if (matrix == null) { throw new ArgumentNullException("matrix"); } if (input == null) { throw new ArgumentNullException("input"); } var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount); Solve(matrix, input, result); return result; } /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. /// /// The coefficient matrix, A. /// The solution matrix, B. /// The result matrix, X public void Solve(Matrix matrix, Matrix input, Matrix result) { if (matrix == null) { throw new ArgumentNullException("matrix"); } if (input == null) { throw new ArgumentNullException("input"); } if (result == null) { throw new ArgumentNullException("result"); } if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount) { throw Matrix.DimensionsDontMatch(matrix, input, result); } for (var column = 0; column < input.ColumnCount; column++) { var solution = Solve(matrix, (Vector)input.Column(column)); foreach (var element in solution.GetIndexedEnumerator()) { result.At(element.Item1, column, element.Item2); } } } } }