// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // // 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 IStation.Numerics.LinearAlgebra.Solvers; namespace IStation.Numerics.LinearAlgebra.Single.Solvers { /// /// 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 { /// /// The collection of solvers that will be used /// readonly List, IPreconditioner>> _solvers; public CompositeSolver(IEnumerable> solvers) { _solvers = solvers.Select(setup => new Tuple, IPreconditioner>(setup.CreateSolver(), setup.CreatePreconditioner() ?? new UnitPreconditioner())).ToList(); } /// /// 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 /// The iterator to use to control when to stop iterating. /// The preconditioner to use for approximations. public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator, IPreconditioner preconditioner) { if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException("Matrix must be square.", nameof(matrix)); } if (result.Count != input.Count) { throw new ArgumentException("All vectors must have the same dimensionality."); } if (iterator == null) { iterator = new Iterator(); } if (preconditioner == null) { preconditioner = new UnitPreconditioner(); } // Create a copy of the solution and result vectors so we can use them // later on var internalInput = input.Clone(); var internalResult = result.Clone(); foreach (var solver in _solvers) { // Store a reference to the solver so we can stop it. IterationStatus status; try { // Reset the iterator and pass it to the solver iterator.Reset(); // Start the solver solver.Item1.Solve(matrix, internalInput, internalResult, iterator, solver.Item2 ?? preconditioner); status = iterator.Status; } 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); continue; } // There was no fatal breakdown so check the status if (status == IterationStatus.Converged) { // We're done internalResult.CopyTo(result); break; } // We're not done // Either: // - calculation finished without convergence if (status == IterationStatus.StoppedWithoutConvergence) { // 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); } } } } }