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
  //
  // 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.
  // </copyright>
  
  using System;
  using Nequeo.Science.Math.LinearAlgebra.Solvers;
  using Nequeo.Science.Math.Properties;
  
  namespace Nequeo.Science.Math.LinearAlgebra.Complex.Solvers
  {
  
  #if NOSYSNUMERICS
      using Complex = Numerics.Complex;
  #else
      using Complex = System.Numerics.Complex;
  #endif
  
      /// <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<Complex>
      {
          /// <summary>
          /// Indicates the number of <c>BiCGStab</c> steps should be taken 
          /// before switching.
          /// </summary>
          int _numberOfBiCgStabSteps = 1;
  
          /// <summary>
          /// Indicates the number of <c>GPBiCG</c> steps should be taken 
          /// before switching.
          /// </summary>
          int _numberOfGpbiCgSteps = 4;
  
          /// <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>
          /// Calculates the <c>true</c> residual of the matrix equation Ax = b according to: residual = b - Ax
          /// </summary>
          /// <param name="matrix">Instance of the <see cref="Matrix"/> A.</param>
          /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
          /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
          /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
          static void CalculateTrueResidual(Matrix<Complex> matrix, Vector<Complex> residual, Vector<Complex> x, Vector<Complex> b)
          {
              // -Ax = residual
              matrix.Multiply(x, residual);
              residual.Multiply(-1, residual);
  
              // residual + b
              residual.Add(b, residual);
          }
  
          /// <summary>
          /// Decide if to do steps with BiCgStab
          /// </summary>
          /// <param name="iterationNumber">Number of iteration</param>
          /// <returns><c>true</c> if yes, otherwise <c>false</c></returns>
          bool ShouldRunBiCgStabSteps(int iterationNumber)
          {
              // Run the first steps as BiCGStab
              // The number of steps past a whole iteration set
              var difference = iterationNumber % (_numberOfBiCgStabSteps + _numberOfGpbiCgSteps);
  
              // Do steps with BiCGStab if:
              // - The difference is zero or more (i.e. we have done zero or more complete cycles)
              // - The difference is less than the number of BiCGStab steps that should be taken
              return (difference >= 0) && (difference < _numberOfBiCgStabSteps);
          }
  
          /// <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>
          /// <param name="iterator">The iterator to use to control when to stop iterating.</param>
          /// <param name="preconditioner">The preconditioner to use for approximations.</param>
          public void Solve(Matrix<Complex> matrix, Vector<Complex> input, Vector<Complex> result, Iterator<Complex> iterator, IPreconditioner<Complex> preconditioner)
          {
              if (matrix.RowCount != matrix.ColumnCount)
              {
                  throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix");
              }
  
              if (input.Count != matrix.RowCount || result.Count != input.Count)
              {
                  throw Matrix.DimensionsDontMatch<ArgumentException>(matrix, input, result);
              }
  
              if (iterator == null)
              {
                  iterator = new Iterator<Complex>();
              }
  
              if (preconditioner == null)
              {
                  preconditioner = new UnitPreconditioner<Complex>();
              }
  
              preconditioner.Initialize(matrix);
  
              // x_0 is initial guess
              // Take x_0 = 0
              var xtemp = new DenseVector(input.Count);
  
              // r_0 = b - Ax_0
              // This is basically a SAXPY so it could be made a lot faster
              var residuals = new DenseVector(matrix.RowCount);
              CalculateTrueResidual(matrix, residuals, xtemp, input);
  
              // Define the temporary scalars
              Complex beta = 0;
  
              // Define the temporary vectors
              // rDash_0 = r_0
              var rdash = DenseVector.OfVector(residuals);
  
              // t_-1 = 0
              var t = new DenseVector(residuals.Count);
              var t0 = new DenseVector(residuals.Count);
  
              // w_-1 = 0
              var w = new DenseVector(residuals.Count);
  
              // Define the remaining temporary vectors
              var c = new DenseVector(residuals.Count);
              var p = new DenseVector(residuals.Count);
              var s = new DenseVector(residuals.Count);
              var u = new DenseVector(residuals.Count);
              var y = new DenseVector(residuals.Count);
              var z = new DenseVector(residuals.Count);
  
              var temp = new DenseVector(residuals.Count);
              var temp2 = new DenseVector(residuals.Count);
              var temp3 = new DenseVector(residuals.Count);
  
              // for (k = 0, 1, .... )
              var iterationNumber = 0;
              while (iterator.DetermineStatus(iterationNumber, xtemp, input, residuals) == IterationStatus.Continue)
              {
                  // p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1))
                  p.Subtract(u, temp);
  
                  temp.Multiply(beta, temp2);
                  residuals.Add(temp2, p);
  
                  // Solve M b_k = p_k
                  preconditioner.Approximate(p, temp);
  
                  // s_k = A b_k
                  matrix.Multiply(temp, s);
  
                  // alpha_k = (r*_0 * r_k) / (r*_0 * s_k)
                  var alpha = rdash.ConjugateDotProduct(residuals)/rdash.ConjugateDotProduct(s);
  
                  // y_k = t_(k-1) - r_k - alpha_k * w_(k-1) + alpha_k s_k
                  s.Subtract(w, temp);
                  t.Subtract(residuals, y);
  
                  temp.Multiply(alpha, temp2);
                  y.Add(temp2, temp3);
                  temp3.CopyTo(y);
  
                  // Store the old value of t in t0
                  t.CopyTo(t0);
  
                  // t_k = r_k - alpha_k s_k
                  s.Multiply(-alpha, temp2);
                  residuals.Add(temp2, t);
  
                  // Solve M d_k = t_k
                  preconditioner.Approximate(t, temp);
  
                  // c_k = A d_k
                  matrix.Multiply(temp, c);
  
gt:                 var cdot = c.ConjugateDotProduct(c);