/// <summary>
        /// Initializes a new instance of the <see cref="DenseEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public static DenseEvd Create(DenseMatrix matrix, Symmetricity symmetricity)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matrices for eigenvalues and eigenvectors
            var eigenVectors = DenseMatrix.CreateIdentity(order);
            var blockDiagonal = new DenseMatrix(order);
            var eigenValues = new DenseVector(order);

            bool isSymmetric;
            switch (symmetricity)
            {
                case Symmetricity.Hermitian:
                    isSymmetric = true;
                    break;
                case Symmetricity.Asymmetric:
                    isSymmetric = false;
                    break;
                default:
                    isSymmetric = matrix.IsHermitian();
                    break;
            }

            Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values);

            return new DenseEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric);
        }

        // Computes coefficients of real polynomial (or estimates if values are noised) using Svd
        // Each row of matrix should contain [x, P(x)], at least rank+1 rows
        // Supplied x-es best have magintude in range [1-2], so resulting coefficient matrix is well conditioned
        public static Polynomial EstimatePolynomial(Matrix<float> values, int rank)
        {
            // 1) Create equation Xa = b
            // | x1^n x1^n-1 ... x1 1 | | a0 |   | P(x1) |
            // |                      | | ...| = |       |
            // | xk^n xk^n-1 ... xk 1 | | an |   | P(xk) |
            Matrix<float> X = new DenseMatrix(values.RowCount, rank + 1);
            Vector<float> P = new DenseVector(values.RowCount);

            for(int i = 0; i < values.RowCount; ++i)
            {
                X[i, rank] = 1.0f;
                for(int c = rank - 1; c >= 0; --c)
                {
                    X[i, c] = X[i, c + 1] * values.At(i, 0);
                }
                P[i] = values[i, 1];
            }

            return new Polynomial()
            {
                Coefficents = SvdSolver.Solve(X, P),
                Rank = rank
            };
        }

        public void HeapSortWithDecreasingDoubleArray()
        {
            var sortedIndices = new int[10];
            Vector<Complex> values = new DenseVector(10);
            values[0] = 9;
            values[1] = 8;
            values[2] = 7;
            values[3] = 6;
            values[4] = 5;
            values[5] = 4;
            values[6] = 3;
            values[7] = 2;
            values[8] = 1;
            values[9] = 0;
            for (var i = 0; i < sortedIndices.Length; i++)
            {
                sortedIndices[i] = i;
            }

            IlutpElementSorter.SortDoubleIndicesDecreasing(0, sortedIndices.Length - 1, sortedIndices, values);
            for (var i = 0; i < sortedIndices.Length; i++)
            {
                Assert.AreEqual(i, sortedIndices[i], "#01-" + i);
            }
        }

        /// <summary>
        /// The main method that uses the Ilutp preconditioner.
        /// </summary>
        public void UseSolver()
        {
            // Create a sparse matrix. For now the size will be 10 x 10 elements
            Matrix matrix = CreateMatrix(10);

            // Create the right hand side vector. The size is the same as the matrix
            // and all values will be 2.0.
            Vector rightHandSideVector = new DenseVector(10, 2.0);

            // Create the Ilutp preconditioner
            Ilutp preconditioner = new Ilutp();

            // Set the drop tolerance. All entries with absolute values smaller than this value will be
            // removed from the preconditioner matrices.
            preconditioner.DropTolerance = 1e-5;
            // Set the relative fill level. This indicates how much additional fill we allow. In this case
            // about 200%
            preconditioner.FillLevel = 200;
            // Set the pivot tolerance. This indicates when pivoting is used. In this case we pivot if
            // the largest off-diagonal entry is twice as big as the diagonal entry.
            preconditioner.PivotTolerance = 0.5;

            // Create the actual preconditioner
            preconditioner.Initialize(matrix);

            // Now that all is set we can solve the matrix equation.
            Vector solutionVector = preconditioner.Approximate(rightHandSideVector);

            // Another way to get the values is by using the overloaded solve method
            // In this case the solution vector needs to be of the correct size.
            preconditioner.Approximate(rightHandSideVector, solutionVector);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="DenseEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public DenseEvd(DenseMatrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matrices for eigenvalues and eigenvectors
            MatrixEv = DenseMatrix.Identity(order);
            MatrixD = matrix.CreateMatrix(order, order);
            VectorEv = new DenseVector(order);

            IsSymmetric = true;

            for (var i = 0; IsSymmetric && i < order; i++)
            {
                for (var j = 0; IsSymmetric && j < order; j++)
                {
                    IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
                }
            }

            Control.LinearAlgebraProvider.EigenDecomp(IsSymmetric, order, matrix.Values, ((DenseMatrix) MatrixEv).Values,
                ((DenseVector) VectorEv).Values, ((DenseMatrix) MatrixD).Values);
        }

        public bool ReadFile(string fileName)
        {
            TextReader tr = null;
            try {
                tr = new StreamReader(fileName);
            } catch (Exception e) {
                return false;
            }

            string buffer = tr.ReadLine();
            string[] numbers = buffer.Split(' ');

            N = int.Parse(numbers[0]);
            P = int.Parse(numbers[1]);

            A = new DenseMatrix(P, N);
            for (int i = 0; i < P; i++) {
                buffer = tr.ReadLine();
                numbers = buffer.Split(' ');
                for (int j = 0; j < N; j++) {
                    A[i,j] = int.Parse(numbers[j]);
                }
            }
            B = new DenseVector(P);
            for (int i = 0; i < P; i++) {
                buffer = tr.ReadLine();
                B[i] = int.Parse(buffer);
            }

            return true;
        }

        public void SolveWideMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(2, 3);
            Vector input = new DenseVector(2);

            var solver = new GpBiCg();
            Assert.Throws<ArgumentException>(() => solver.Solve(matrix, input));
        }

        public void SolveLongMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(3, 2);
            Vector input = new DenseVector(3);

            var solver = new BiCgStab();
            Assert.Throws<ArgumentException>(() => solver.Solve(matrix, input));
        }

        public void HeapSortWithDuplicateDoubleEntries()
        {
            var sortedIndices = new int[10];
            Vector<Complex> values = new DenseVector(10);
            values[0] = 1;
            values[1] = 1;
            values[2] = 1;
            values[3] = 1;
            values[4] = 2;
            values[5] = 2;
            values[6] = 2;
            values[7] = 2;
            values[8] = 3;
            values[9] = 4;

            for (var i = 0; i < sortedIndices.Length; i++)
            {
                sortedIndices[i] = i;
            }
            IlutpElementSorter.SortDoubleIndicesDecreasing(0, sortedIndices.Length - 1, sortedIndices, values);
            for (var i = 0; i < sortedIndices.Length; i++)
            {
                if (i == 0)
                {
                    Assert.AreEqual(9, sortedIndices[i], "#01-" + i);
                }
                else
                {
                    if (i == 1)
                    {
                        Assert.AreEqual(8, sortedIndices[i], "#01-" + i);
                    }
                    else
                    {
                        if (i < 6)
                        {
                            if ((sortedIndices[i] != 4) &&
                                (sortedIndices[i] != 5) &&
                                (sortedIndices[i] != 6) &&
                                (sortedIndices[i] != 7))
                            {
                                Assert.Fail("#01-" + i);
                            }
                        }
                        else
                        {
                            if ((sortedIndices[i] != 0) &&
                                (sortedIndices[i] != 1) &&
                                (sortedIndices[i] != 2) &&
                                (sortedIndices[i] != 3))
                            {
                                Assert.Fail("#01-" + i);
                            }
                        }
                    }
                }
            }
        }

 public void CanCallUnaryNegationOperatorOnDenseVector()
 {
     var vector = new DenseVector(this._data);
     var other = -vector;
     for (var i = 0; i < this._data.Length; i++)
     {
         Assert.AreEqual(-this._data[i], other[i]);
     }
 }

 public void CanCallUnaryPlusOperatorOnDenseVector()
 {
     var vector = new DenseVector(Data);
     var other = +vector;
     for (var i = 0; i < Data.Length; i++)
     {
         AssertHelpers.AreEqual(vector[i], other[i]);
     }
 }

        /// <summary>
        /// Create standard vector.
        /// </summary>
        /// <param name="size">Size of the vector.</param>
        /// <returns>New vector.</returns>
        protected Vector CreateStandardBcVector(int size)
        {
            Vector vector = new DenseVector(size);
            for (var i = 0; i < size; i++)
            {
                vector[i] = i + 1;
            }

            return vector;
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="DenseSvd"/> class. This object will compute the
        /// the singular value decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If SVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public static DenseSvd Create(DenseMatrix matrix, bool computeVectors)
        {
            var nm = Math.Min(matrix.RowCount, matrix.ColumnCount);
            var s = new DenseVector(nm);
            var u = new DenseMatrix(matrix.RowCount);
            var vt = new DenseMatrix(matrix.ColumnCount);
            Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix) matrix.Clone()).Values, matrix.RowCount, matrix.ColumnCount, s.Values, u.Values, vt.Values);

            return new DenseSvd(s, u, vt, computeVectors);
        }

        /// <summary>
        /// Creates a new instance of the Vector class.
        /// </summary>
        /// <param name="data">The array to create this vector from.</param>
        /// <returns>The new <c>Vector</c>.</returns>
        protected override Vector<double> CreateVector(IList<double> data)
        {
            var vector = new DenseVector(data.Count);
            for (var index = 0; index < data.Count; index++)
            {
                vector[index] = data[index];
            }

            return vector;
        }

        /// <summary>
        /// The main method that runs the BiCGStab iterative solver.
        /// </summary>
        public void UseSolver()
        {
            // Create a sparse matrix. For now the size will be 10 x 10 elements
            Matrix matrix = CreateMatrix(10);

            // Create the right hand side vector. The size is the same as the matrix
            // and all values will be 2.0.
            Vector rightHandSideVector = new DenseVector(10, 2.0);

            // Create a preconditioner. The possibilities are:
            // 1) No preconditioner - Simply do not provide the solver with a preconditioner.
            // 2) A simple diagonal preconditioner - Create an instance of the Diagonal class.
            // 3) A ILU preconditioner - Create an instance of the IncompleteLu class.
            // 4) A ILU preconditioner with pivoting and drop tolerances - Create an instance of the Ilutp class.

            // Here we'll use the simple diagonal preconditioner.
            // We need a link to the matrix so the pre-conditioner can do it's work.
            IPreConditioner preconditioner = new Diagonal();

            // Create a new iterator. This checks for convergence of the results of the
            // iterative matrix solver.
            // In this case we'll create the default iterator
            IIterator iterator = Iterator.CreateDefault();

            // Create the solver
            BiCgStab solver = new BiCgStab(preconditioner, iterator);

            // Now that all is set we can solve the matrix equation.
            Vector solutionVector = solver.Solve(matrix, rightHandSideVector);

            // Another way to get the values is by using the overloaded solve method
            // In this case the solution vector needs to be of the correct size.
            solver.Solve(matrix, rightHandSideVector, solutionVector);

            // Finally you can check the reason the solver finished the iterative process
            // by calling the SolutionStatus property on the iterator
            ICalculationStatus status = iterator.Status;
            if (status is CalculationCancelled)
                Console.WriteLine("The user cancelled the calculation.");

            if (status is CalculationIndetermined)
                Console.WriteLine("Oh oh, something went wrong. The iterative process was never started.");

            if (status is CalculationConverged)
                Console.WriteLine("Yippee, the iterative process converged.");

            if (status is CalculationDiverged)
                Console.WriteLine("I'm sorry the iterative process diverged.");

            if (status is CalculationFailure)
                Console.WriteLine("Oh dear, the iterative process failed.");

            if (status is CalculationStoppedWithoutConvergence)
                Console.WriteLine("Oh dear, the iterative process did not converge.");
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public static UserEvd Create(Matrix<Complex> matrix, Symmetricity symmetricity)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            var eigenVectors = DenseMatrix.CreateIdentity(order);
            var blockDiagonal = Matrix<Complex>.Build.SameAs(matrix, order, order);
            var eigenValues = new DenseVector(order);

            bool isSymmetric;
            switch (symmetricity)
            {
                case Symmetricity.Hermitian:
                    isSymmetric = true;
                    break;
                case Symmetricity.Asymmetric:
                    isSymmetric = false;
                    break;
                default:
                    isSymmetric = matrix.IsHermitian();
                    break;
            }

            if (isSymmetric)
            {
                var matrixCopy = matrix.ToArray();
                var tau = new Complex[order];
                var d = new double[order];
                var e = new double[order];

                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
                SymmetricDiagonalize(eigenVectors, d, e, order);
                SymmetricUntridiagonalize(eigenVectors, matrixCopy, tau, order);

                for (var i = 0; i < order; i++)
                {
                    eigenValues[i] = new Complex(d[i], e[i]);
                }
            }
            else
            {
                var matrixH = matrix.ToArray();
                NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(eigenVectors, eigenValues, matrixH, order);
            }

            blockDiagonal.SetDiagonal(eigenValues);

            return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="DenseEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public DenseEvd(DenseMatrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            MatrixEv = DenseMatrix.Identity(order);
            MatrixD = matrix.CreateMatrix(order, order);
            VectorEv = new DenseVector(order);

            IsSymmetric = true;

            for (var i = 0; i < order & IsSymmetric; i++)
            {
                for (var j = 0; j < order & IsSymmetric; j++)
                {
                    IsSymmetric &= matrix[i, j] == matrix[j, i].Conjugate();
                }
            }

            if (IsSymmetric)
            {
                var matrixCopy = matrix.ToArray();
                var tau = new Complex[order];
                var d = new double[order];
                var e = new double[order];

                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
                SymmetricDiagonalize(((DenseMatrix)MatrixEv).Data, d, e, order);
                SymmetricUntridiagonalize(((DenseMatrix)MatrixEv).Data, matrixCopy, tau, order);

                for (var i = 0; i < order; i++)
                {
                    VectorEv[i] = new Complex(d[i], e[i]);
                }
            }
            else
            {
                var matrixH = matrix.ToArray();
                NonsymmetricReduceToHessenberg(((DenseMatrix)MatrixEv).Data, matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(((DenseVector)VectorEv).Data, ((DenseMatrix)MatrixEv).Data, matrixH, order);
            }

            MatrixD.SetDiagonal(VectorEv);
        }

 public void CanCreateVectorFromArray()
 {
     var data = new[] { 1.0, 2.0, 3.0, 4.0 };
     var vector = new DenseVector(data);
     Assert.AreSame(data, vector.Data);
     for( var i = 0; i < data.Length; i++)
     {
         Assert.AreEqual(data[i], vector[i]);
     }
     vector[0] = 100.0;
     Assert.AreEqual(100.0, data[0]);
 }

        public void DetermineStatus()
        {
            var criterium = new FailureStopCriterium();
            Assert.IsNotNull(criterium, "There should be a criterium");

            var solution = new DenseVector(new[] { new Complex(3.0, 0), new Complex(2.0, 0), new Complex(1, 0) });
            var source = new DenseVector(new[] { new Complex(1001.0, 0), Complex.Zero, new Complex(2003.0, 0) });
            var residual = new DenseVector(new[] { new Complex(1.0, 0), new Complex(2.0, 0), new Complex(3, 0) });

            criterium.DetermineStatus(5, solution, source, residual);
            Assert.IsInstanceOf(typeof(CalculationRunning), criterium.Status, "Should be running");
        }

        public void ApproximateWithNullVector()
        {
            const int Size = 10;
            var newMatrix = CreateUnitMatrix(Size);
            var vector = CreateStandardBcVector(Size);

            var preconditioner = CreatePreconditioner();
            preconditioner.Initialize(newMatrix);

            Vector<double> result = new DenseVector(vector.Count + 10);
            preconditioner.Approximate(null, result);
        }