/// <summary>
 /// Check the result.
 /// </summary>
 /// <param name="preconditioner">Specific preconditioner.</param>
 /// <param name="matrix">Source matrix.</param>
 /// <param name="vector">Initial vector.</param>
 /// <param name="result">Result vector.</param>
 protected override void CheckResult(IPreConditioner preconditioner, SparseMatrix matrix, Vector vector, Vector result)
 {
     Assert.AreEqual(typeof(UnitPreconditioner), preconditioner.GetType(), "#01");
     // Unit preconditioner is doing nothing. Vector and result should be equal
     for (var i = 0; i < vector.Count; i++)
     {
         Assert.IsTrue(vector[i] == result[i], "#02-" + i);
     }
 }

        /// <summary>
        /// Check the result.
        /// </summary>
        /// <param name="preconditioner">Specific preconditioner.</param>
        /// <param name="matrix">Source matrix.</param>
        /// <param name="vector">Initial vector.</param>
        /// <param name="result">Result vector.</param>
        protected override void CheckResult(IPreConditioner preconditioner, SparseMatrix matrix, Vector vector, Vector result)
        {
            Assert.AreEqual(typeof(Diagonal), preconditioner.GetType(), "#01");

            // Compute M * result = product
            // compare vector and product. Should be equal
            Vector product = new DenseVector(result.Count);
            matrix.Multiply(result, product);
            for (var i = 0; i < product.Count; i++)
            {
                Assert.IsTrue(vector[i].AlmostEqual(product[i], -Epsilon.Magnitude()), "#02-" + i);
            }
        }

 /// <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>
 /// Initializes a new instance of the <see cref="BiCgStab"/> 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 BiCgStab(IPreConditioner preconditioner, IIterator iterator)
 {
     _iterator = iterator;
     _preconditioner = preconditioner;
 }

 /// <summary>
 /// Check the result.
 /// </summary>
 /// <param name="preconditioner">Specific preconditioner.</param>
 /// <param name="matrix">Source matrix.</param>
 /// <param name="vector">Initial vector.</param>
 /// <param name="result">Result vector.</param>
 protected abstract void CheckResult(IPreConditioner preconditioner, SparseMatrix matrix, Vector vector, Vector result);

 /// <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>
 /// Check the result.
 /// </summary>
 /// <param name="preconditioner">Specific preconditioner.</param>
 /// <param name="matrix">Source matrix.</param>
 /// <param name="vector">Initial vector.</param>
 /// <param name="result">Result vector.</param>
 protected abstract void CheckResult(IPreConditioner<float> preconditioner, SparseMatrix matrix, Vector<float> vector, Vector<float> result);

 /// <summary>
 /// Initializes a new instance of the <see cref="TFQMR"/> 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 TFQMR(IPreConditioner preconditioner, IIterator iterator)
 {
     _iterator = iterator;
     _preconditioner = preconditioner;
 }

 protected abstract void CheckResult(IPreConditioner<Complex32> preconditioner, SparseMatrix matrix, Vector<Complex32> vector, Vector<Complex32> 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 <see cref="Matrix"/>, <c>A</c>.</param>
        /// <param name="input">The solution <see cref="Vector"/>, <c>b</c>.</param>
        /// <param name="result">The result <see cref="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;

            // Parameters 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 Matrix.DimensionsDontMatch<ArgumentException>(input, result);
            }

            // Initialize the solver fields
            // Set the convergence monitor
            if (_iterator == null)
            {
                _iterator = Iterator.CreateDefault();
            }

            if (_preconditioner == null)
            {
                _preconditioner = new UnitPreconditioner();
            }
            
            _preconditioner.Initialize(matrix);
            
            // Compute r_0 = b - Ax_0 for some initial guess x_0
            // In this case we take x_0 = vector
            // This is basically a SAXPY so it could be made a lot faster
            Vector residuals = new DenseVector(matrix.RowCount);
            CalculateTrueResidual(matrix, residuals, result, input);

            // Choose r~ (for example, r~ = r_0)
            var tempResiduals = residuals.Clone();

            // create seven temporary vectors needed to hold temporary
            // coefficients. All vectors are mangled in each iteration.
            // These are defined here to prevent stressing the garbage collector
            Vector vecP = new DenseVector(residuals.Count);
            Vector vecPdash = new DenseVector(residuals.Count);
            Vector nu = new DenseVector(residuals.Count);
            Vector vecS = new DenseVector(residuals.Count);
            Vector vecSdash = new DenseVector(residuals.Count);
            Vector temp = new DenseVector(residuals.Count);
            Vector temp2 = new DenseVector(residuals.Count);

            // create some temporary double variables that are needed
            // to hold values in between iterations
            Complex currentRho = 0;
            Complex alpha = 0;
            Complex omega = 0;

            var iterationNumber = 0;
            while (ShouldContinue(iterationNumber, result, input, residuals))
            {
                // rho_(i-1) = r~^T r_(i-1) // dotproduct r~ and r_(i-1)
                var oldRho = currentRho;
                currentRho = tempResiduals.DotProduct(residuals);

                // if (rho_(i-1) == 0) // METHOD FAILS
                // If rho is only 1 ULP from zero then we fail.
                if (currentRho.Real.AlmostEqual(0, 1) && currentRho.Imaginary.AlmostEqual(0, 1))
                {
                    // Rho-type breakdown
                    throw new Exception("Iterative solver experience a numerical break down");
                }

//.........这里部分代码省略.........

        /// <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 Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
            }

            // 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
            float beta = 0;
            float 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);
            Vector p = new DenseVector(residuals.Count);
            Vector s = new DenseVector(residuals.Count);
            Vector u = new DenseVector(residuals.Count);
            Vector y = new DenseVector(residuals.Count);
            Vector z = new DenseVector(residuals.Count);

            Vector temp = new DenseVector(residuals.Count);
            Vector temp2 = new DenseVector(residuals.Count);
            Vector temp3 = new DenseVector(residuals.Count);

            // for (k = 0, 1, .... )
            var iterationNumber = 0;
            while (ShouldContinue(iterationNumber, xtemp, input, residuals))
            {
                // p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1))
//.........这里部分代码省略.........

        /// <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 Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
            }

            // Initialize the solver fields
            // Set the convergence monitor
            if (_iterator == null)
            {
                _iterator = Iterator.CreateDefault();
            }

            if (_preconditioner == null)
            {
                _preconditioner = new UnitPreconditioner();
            }

            _preconditioner.Initialize(matrix);

            var d = new DenseVector(input.Count);
            var r = new DenseVector(input);

            var uodd = new DenseVector(input.Count);
            var ueven = new DenseVector(input.Count);

            var v = new DenseVector(input.Count);
            var pseudoResiduals = new DenseVector(input);

            var x = new DenseVector(input.Count);
            var yodd = new DenseVector(input.Count);
            var yeven = new DenseVector(input);

            // Temp vectors
            var temp = new DenseVector(input.Count);
            var temp1 = new DenseVector(input.Count);
            var temp2 = new DenseVector(input.Count);

            // Initialize
            var startNorm = input.Norm(2);

            // Define the scalars
            double alpha = 0;
            double eta = 0;
            double theta = 0;

            var tau = startNorm;
            var rho = tau * tau;

            // Calculate the initial values for v
            // M temp = yEven
            _preconditioner.Approximate(yeven, temp);

            // v = A temp
            matrix.Multiply(temp, v);

            // Set uOdd
            v.CopyTo(ueven);

            // Start the iteration
            var iterationNumber = 0;
//.........这里部分代码省略.........

        /// <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);

            // Choose an initial guess x_0
            // Take x_0 = 0
            Vector xtemp = new DenseVector(input.Count);

            // Choose k vectors q_1, q_2, ..., q_k
            // Build a new set if:
            // a) the stored set doesn't exist (i.e. == null)
            // b) Is of an incorrect length (i.e. too long)
            // c) The vectors are of an incorrect length (i.e. too long or too short)
            var useOld = false;
            if (_startingVectors != null)
            {
                // We don't accept collections with zero starting vectors so ...
                if (_startingVectors.Count <= NumberOfStartingVectorsToCreate(_numberOfStartingVectors, input.Count))
                {
                    // Only check the first vector for sizing. If that matches we assume the
                    // other vectors match too. If they don't the process will crash
                    if (_startingVectors[0].Count == input.Count)
                    {
                        useOld = true;
                    }
                }
            }

            _startingVectors = useOld ? _startingVectors : CreateStartingVectors(_numberOfStartingVectors, input.Count);

            // Store the number of starting vectors. Not really necessary but easier to type :)
            var k = _startingVectors.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 values
            var c = new Complex[k];

            // Define the temporary vectors
            Vector gtemp = new DenseVector(residuals.Count);

            Vector u = new DenseVector(residuals.Count);
//.........这里部分代码省略.........