context: // <copyright file="GpBiCg.cs" company="Math.NET">
  // 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.
  // </copyright>
  
  namespace Nequeo.Science.Math.LinearAlgebra.Double.Solvers.Iterative
  {
      using System;
      using Generic.Solvers.Status;
      using Preconditioners;
      using Properties;
  
      /// <summary>
      /// A Generalized Product Bi-Conjugate Gradient iterative matrix solver.
      /// </summary>
      /// <remarks>
      /// <para>
      /// The Generalized Product Bi-Conjugate Gradient (GPBiCG) solver is an 
      /// alternative version of the Bi-Conjugate Gradient stabilized (CG) solver.
      /// Unlike the CG solver the GPBiCG solver can be used on 
      /// non-symmetric matrices. <br/>
      /// Note that much of the success of the solver depends on the selection of the
      /// proper preconditioner.
      /// </para>
      /// <para>
      /// The GPBiCG algorithm was taken from: <br/>
      /// GPBiCG(m,l): A hybrid of BiCGSTAB and GPBiCG methods with 
      /// efficiency and robustness
      /// <br/>
      /// S. Fujino
      /// <br/>
      /// Applied Numerical Mathematics, Volume 41, 2002, pp 107 - 117
      /// <br/>
      /// </para>
      /// <para>
      /// The example code below provides an indication of the possible use of the
      /// solver.
      /// </para>
      /// </remarks>
      public sealed class GpBiCg : IIterativeSolver
      {
          /// <summary>
          /// The status used if there is no status, i.e. the solver hasn't run yet and there is no
          /// iterator.
          /// </summary>
          private static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined();
  
          /// <summary>
          /// The preconditioner that will be used. Can be set to <c>null</c>, in which case the default
          /// pre-conditioner will be used.
          /// </summary>
          private IPreConditioner _preconditioner;
  
          /// <summary>
          /// The iterative process controller.
          /// </summary>
          private IIterator _iterator;
  
          /// <summary>
          /// Indicates the number of <c>BiCGStab</c> steps should be taken 
          /// before switching.
          /// </summary>
          private int _numberOfBiCgStabSteps = 1;
  
          /// <summary>
          /// Indicates the number of <c>GPBiCG</c> steps should be taken 
          /// before switching.
          /// </summary>
          private int _numberOfGpbiCgSteps = 4;
  
          /// <summary>
          /// Indicates if the user has stopped the solver.
          /// </summary>
          private bool _hasBeenStopped;
  
          /// <summary>
          /// Initializes a new instance of the <see cref="GpBiCg"/> class.
          /// </summary>
          /// <remarks>
          /// When using this constructor the solver will use the <see cref="IIterator"/> with
          /// the standard settings and a default preconditioner.
          /// </remarks>
          public GpBiCg() : this(null, null)
          {
          }
  
          /// <summary>
          /// Initializes a new instance of the <see cref="GpBiCg"/> class.
          /// </summary>
          /// <remarks>
          /// <para>
          /// When using this constructor the solver will use a default preconditioner.
          /// </para>
          /// <para>
          /// The main advantages of using a user defined <see cref="IIterator"/> are:
          /// <list type="number">
          /// <item>It is possible to set the desired convergence limits.</item>
          /// <item>
          /// It is possible to check the reason for which the solver finished 
          /// the iterative procedure by calling the <see cref="IIterator.Status"/> property.
          /// </item>
          /// </list>
          /// </para>
          /// </remarks>
          /// <param name="iterator">The <see cref="IIterator"/> that will be used to monitor the iterative process.</param>
          public GpBiCg(IIterator iterator) : this(null, iterator)
          {
          }
  
          /// <summary>
          /// Initializes a new instance of the <see cref="GpBiCg"/> class.
          /// </summary>
          /// <remarks>
          /// When using this constructor the solver will use the <see cref="IIterator"/> with
          /// the standard settings.
          /// </remarks>
          /// <param name="preconditioner">The <see cref="IPreConditioner"/> that will be used to precondition the matrix equation.</param>
          public GpBiCg(IPreConditioner preconditioner) : this(preconditioner, null)
          {
          }
  
          /// <summary>
          /// Initializes a new instance of the <see cref="GpBiCg"/> class.
          /// </summary>
          /// <remarks>
          /// <para>
          /// The main advantages of using a user defined <see cref="IIterator"/> are:
          /// <list type="number">
          /// <item>It is possible to set the desired convergence limits.</item>
          /// <item>
          /// It is possible to check the reason for which the solver finished 
          /// the iterative procedure by calling the <see cref="IIterator.Status"/> property.
          /// </item>
          /// </list>
          /// </para>
          /// </remarks>
          /// <param name="preconditioner">The <see cref="IPreConditioner"/> that will be used to precondition the matrix equation.</param>
          /// <param name="iterator">The <see cref="IIterator"/> that will be used to monitor the iterative process.</param>
          public GpBiCg(IPreConditioner preconditioner, IIterator iterator)
          {
              _iterator = iterator;
              _preconditioner = preconditioner;
          }
  
          /// <summary>
          /// Gets or sets the number of steps taken with the <c>BiCgStab</c> algorithm
          /// before switching over to the <c>GPBiCG</c> algorithm.
          /// </summary>
          public int NumberOfBiCgStabSteps
          {
              get
              {
                  return _numberOfBiCgStabSteps;
              }
  
              set
              {
                  if (value < 0)
                  {
                      throw new ArgumentOutOfRangeException("value");
                  }
  
                  _numberOfBiCgStabSteps = value;     
              }
          }
  
          /// <summary>
          /// Gets or sets the number of steps taken with the <c>GPBiCG</c> algorithm
          /// before switching over to the <c>BiCgStab</c> algorithm.
          /// </summary>
          public int NumberOfGpBiCgSteps
          {
              get
              {
                  return _numberOfGpbiCgSteps;
              }
  
              set
              {
                  if (value < 0)
                  {
                      throw new ArgumentOutOfRangeException("value");
                  }
  
                  _numberOfGpbiCgSteps = value;
              }
          }
  
          /// <summary>
          /// Sets the <see cref="IPreConditioner"/> that will be used to precondition the iterative process.
          /// </summary>
          /// <param name="preconditioner">The preconditioner.</param>
          public void SetPreconditioner(IPreConditioner preconditioner)
          {
              _preconditioner = preconditioner;
          }
  
          /// <summary>
          /// Sets the <see cref="IIterator"/> that will be used to track the iterative process.
          /// </summary>
          /// <param name="iterator">The iterator.</param>
          public void SetIterator(IIterator iterator)
          {
              _iterator = iterator;
          }
  
          /// <summary>
          /// Gets the status of the iteration once the calculation is finished.
          /// </summary>
          public ICalculationStatus IterationResult
          {
              get 
              {
                  return (_iterator != null) ? _iterator.Status : DefaultStatus;
              }
          }
  
          /// <summary>
          /// Stops the solve process. 
          /// </summary>
          /// <remarks>
          /// Note that it may take an indetermined amount of time for the solver to actually
          /// stop the process.
          /// </remarks>
          public void StopSolve()
          {
              _hasBeenStopped = true;
          }
  
          /// <summary>
          /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
          /// solution vector and x is the unknown vector.
          /// </summary>
          /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
          /// <param name="vector">The solution vector, <c>b</c>.</param>
          /// <returns>The result vector, <c>x</c>.</returns>
          public Vector Solve(Matrix matrix, Vector vector)
          {
              if (vector == null)
              {
                  throw new ArgumentNullException();
              }
  
              Vector result = new DenseVector(matrix.RowCount);
              Solve(matrix, vector, result);
              return result;
          }
  
          /// <summary>
          /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
          /// solution vector and x is the unknown vector.
          /// </summary>
          /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
          /// <param name="input">The solution vector, <c>b</c></param>
          /// <param name="result">The result vector, <c>x</c></param>
          public void Solve(Matrix 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;
  
              // 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);
              }
  
              if (input.Count != matrix.RowCount)
              {
                  throw new ArgumentException(Resources.ArgumentMatrixDimensions);
              }
  
              // Initialize the solver fields
              // Set the convergence monitor
              if (_iterator == null)
              {
                  _iterator = Iterator.CreateDefault();
              }
  
              if (_preconditioner == null)
              {
                  _preconditioner = new UnitPreconditioner();
              }
           
              _preconditioner.Initialize(matrix);
  
              // x_0 is initial guess
              // Take x_0 = 0
              Vector xtemp = new DenseVector(input.Count);
  
              // r_0 = b - Ax_0
              // This is basically a SAXPY so it could be made a lot faster
              Vector residuals = new DenseVector(matrix.RowCount);
              CalculateTrueResidual(matrix, residuals, xtemp, input);
  
              // Define the temporary scalars
              double beta = 0;
              double sigma;
  
              // Define the temporary vectors
              // rDash_0 = r_0
              Vector rdash = new DenseVector(residuals);
  
              // t_-1 = 0
              Vector t = new DenseVector(residuals.Count);
              Vector t0 = new DenseVector(residuals.Count);
  
              // w_-1 = 0
              Vector w = new DenseVector(residuals.Count);
  
              // Define the remaining temporary vectors
              Vector c = new DenseVector(residuals.Count);
  
gt:             Vector p = new DenseVector(residuals.Count);