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

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

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

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

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

            var solver = new TFQMR();
            Assert.That(() => matrix.SolveIterative(input, solver), Throws.ArgumentException);
        }

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

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

        /// <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<Complex> preconditioner, SparseMatrix matrix, Vector<Complex> vector, Vector<Complex> result)
        {
            Assert.AreEqual(typeof(UnitPreconditioner<Complex>), 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>
        /// Create unit matrix.
        /// </summary>
        /// <param name="size">Matrix size.</param>
        /// <returns>New unit matrix.</returns>
        internal SparseMatrix CreateUnitMatrix(int size)
        {
            var matrix = new SparseMatrix(size);
            for (var i = 0; i < size; i++)
            {
                matrix[i, i] = 2;
            }

            return matrix;
        }

        /// <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<Complex> preconditioner, SparseMatrix matrix, Vector<Complex> vector, Vector<Complex> result)
        {
            Assert.AreEqual(typeof(Diagonal), preconditioner.GetType(), "#01");

            // Compute M * result = product
            // compare vector and product. Should be equal
            var product = new DenseVector(result.Count);
            matrix.Multiply(result, product);

            for (var i = 0; i < product.Count; i++)
            {
                Assert.IsTrue(vector[i].Real.AlmostEqual(product[i].Real, -Epsilon.Magnitude()), "#02-" + i);
                Assert.IsTrue(vector[i].Imaginary.AlmostEqual(product[i].Imaginary, -Epsilon.Magnitude()), "#03-" + i);
            }
        }

        public void CanAddSparseMatricesBothWays()
        {
            var m1 = new SparseMatrix(1, 3);
            var m2 = SparseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            var sum1 = m1 + m2;
            var sum2 = m2 + m1;
            Assert.IsTrue(sum1.Equals(m2));
            Assert.IsTrue(sum1.Equals(sum2));

            var sparseResult = new SparseMatrix(1, 3);
            sparseResult.Add(m2, sparseResult);
            Assert.IsTrue(sparseResult.Equals(sum1));

            sparseResult = SparseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            sparseResult.Add(m1, sparseResult);
            Assert.IsTrue(sparseResult.Equals(sum1));

            sparseResult = SparseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            m1.Add(sparseResult, sparseResult);
            Assert.IsTrue(sparseResult.Equals(sum1));

            sparseResult = SparseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            sparseResult.Add(sparseResult, sparseResult);
            Assert.IsTrue(sparseResult.Equals(2*sum1));

            var denseResult = new DenseMatrix(1, 3);
            denseResult.Add(m2, denseResult);
            Assert.IsTrue(denseResult.Equals(sum1));

            denseResult = DenseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            denseResult.Add(m1, denseResult);
            Assert.IsTrue(denseResult.Equals(sum1));

            var m3 = DenseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            var sum3 = m1 + m3;
            var sum4 = m3 + m1;
            Assert.IsTrue(sum3.Equals(m3));
            Assert.IsTrue(sum3.Equals(sum4));
        }

        /// <summary>
        /// Create a matrix based on this vector in row form (one single row).
        /// </summary>
        /// <returns>This vector as a row matrix.</returns>
        public override Matrix<Complex> ToRowMatrix()
        {
            var matrix = new SparseMatrix(1, Count);
            for (var i = 0; i < NonZerosCount; i++)
            {
                matrix.At(0, _nonZeroIndices[i], _nonZeroValues[i]);
            }

            return matrix;
        }

        /// <summary>
        /// Create a matrix based on this vector in column form (one single column).
        /// </summary>
        /// <returns>This vector as a column matrix.</returns>
        public override Matrix<Complex> ToColumnMatrix()
        {
            var matrix = new SparseMatrix(Count, 1);
            for (var i = 0; i < NonZerosCount; i++)
            {
                matrix.At(_nonZeroIndices[i], 0, _nonZeroValues[i]);
            }

            return matrix;
        }

 /// <summary>
 /// Create a new sparse matrix with the diagonal as a copy of the given array.
 /// This new matrix will be independent from the array.
 /// A new memory block will be allocated for storing the matrix.
 /// </summary>
 public static SparseMatrix OfDiagonalArray(int rows, int columns, Complex[] diagonal)
 {
     var m = new SparseMatrix(rows, columns);
     m.SetDiagonal(diagonal);
     return m;
 }

 /// <summary>
 /// Create a new sparse matrix with the diagonal as a copy of the given vector.
 /// This new matrix will be independent from the vector.
 /// A new memory block will be allocated for storing the matrix.
 /// </summary>
 public static SparseMatrix OfDiagonalVector(int rows, int columns, Vector<Complex> diagonal)
 {
     var m = new SparseMatrix(rows, columns);
     m.SetDiagonal(diagonal);
     return m;
 }

        /// <summary>
        /// Can multiply a matrix with matrix.
        /// </summary>
        /// <param name="nameA">Matrix A name.</param>
        /// <param name="nameB">Matrix B name.</param>
        public override void CanMultiplyMatrixWithMatrixIntoResult(string nameA, string nameB)
        {
            var matrixA = TestMatrices[nameA];
            var matrixB = TestMatrices[nameB];
            var matrixC = new SparseMatrix(matrixA.RowCount, matrixB.ColumnCount);
            matrixA.Multiply(matrixB, matrixC);

            Assert.AreEqual(matrixC.RowCount, matrixA.RowCount);
            Assert.AreEqual(matrixC.ColumnCount, matrixB.ColumnCount);

            for (var i = 0; i < matrixC.RowCount; i++)
            {
                for (var j = 0; j < matrixC.ColumnCount; j++)
                {
                    AssertHelpers.AlmostEqualRelative(matrixA.Row(i)*matrixB.Column(j), matrixC[i, j], 15);
                }
            }
        }

        public void CanSubtractSparseMatricesBothWays()
        {
            var m1 = new SparseMatrix(1, 3);
            var m2 = SparseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            var diff1 = m1 - m2;
            var diff2 = m2 - m1;
            Assert.IsTrue(diff1.Equals(m2.Negate()));
            Assert.IsTrue(diff1.Equals(diff2.Negate()));

            var sparseResult = new SparseMatrix(1, 3);
            sparseResult.Subtract(m2, sparseResult);
            Assert.IsTrue(sparseResult.Equals(diff1));

            sparseResult = SparseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            sparseResult.Subtract(m1, sparseResult);
            Assert.IsTrue(sparseResult.Equals(diff2));

            sparseResult = SparseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            m1.Subtract(sparseResult, sparseResult);
            Assert.IsTrue(sparseResult.Equals(diff1));

            sparseResult = SparseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            sparseResult.Subtract(sparseResult, sparseResult);
            Assert.IsTrue(sparseResult.Equals(0*diff1));

            var denseResult = new DenseMatrix(1, 3);
            denseResult.Subtract(m2, denseResult);
            Assert.IsTrue(denseResult.Equals(diff1));

            denseResult = DenseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            denseResult.Subtract(m1, denseResult);
            Assert.IsTrue(denseResult.Equals(diff2));

            var m3 = DenseMatrix.OfArray(new Complex[,] {{0, 1, 1}});
            var diff3 = m1 - m3;
            var diff4 = m3 - m1;
            Assert.IsTrue(diff3.Equals(m3.Negate()));
            Assert.IsTrue(diff3.Equals(diff4.Negate()));
        }

 public void CanCreateLargeMatrix()
 {
     const int Order = 1000000;
     var matrix = new SparseMatrix(Order);
     Assert.AreEqual(Order, matrix.RowCount);
     Assert.AreEqual(Order, matrix.ColumnCount);
     Assert.DoesNotThrow(() => matrix[0, 0] = 1);
 }

        public void CanCreateLargeSparseMatrix()
        {
            var matrix = new SparseMatrix(500, 1000);
            var nonzero = 0;
            var rnd = new System.Random(0);

            for (var i = 0; i < matrix.RowCount; i++)
            {
                for (var j = 0; j < matrix.ColumnCount; j++)
                {
                    var value = rnd.Next(10)*rnd.Next(10)*rnd.Next(10)*rnd.Next(10)*rnd.Next(10);
                    if (value != 0)
                    {
                        nonzero++;
                    }

                    matrix[i, j] = value;
                }
            }

            Assert.AreEqual(matrix.NonZerosCount, nonzero);
        }

        /// <summary>
        /// Create a matrix based on this vector in row form (one single row).
        /// </summary>
        /// <returns>This vector as a row matrix.</returns>
        public override Matrix<Complex> ToRowMatrix()
        {
            var matrix = new SparseMatrix(1, Count);
            for (var i = 0; i < _storage.ValueCount; i++)
            {
                matrix.At(0, _storage.Indices[i], _storage.Values[i]);
            }

            return matrix;
        }

        /// <summary>
        /// Create a matrix based on this vector in column form (one single column).
        /// </summary>
        /// <returns>This vector as a column matrix.</returns>
        public override Matrix<Complex> ToColumnMatrix()
        {
            var matrix = new SparseMatrix(Count, 1);
            for (var i = 0; i < _storage.ValueCount; i++)
            {
                matrix.At(_storage.Indices[i], 0, _storage.Values[i]);
            }

            return matrix;
        }

        /// <summary>
        /// Outer product of two vectors
        /// </summary>
        /// <param name="u">First vector</param>
        /// <param name="v">Second vector</param>
        /// <returns>Matrix M[i,j] = u[i]*v[j] </returns>
        /// <exception cref="ArgumentNullException">If the u vector is <see langword="null" />.</exception> 
        /// <exception cref="ArgumentNullException">If the v vector is <see langword="null" />.</exception> 
        public static Matrix<Complex> OuterProduct(SparseVector u, SparseVector v)
        {
            if (u == null)
            {
                throw new ArgumentNullException("u");
            }

            if (v == null)
            {
                throw new ArgumentNullException("v");
            }

            var matrix = new SparseMatrix(u.Count, v.Count);
            for (var i = 0; i < u.NonZerosCount; i++)
            {
                for (var j = 0; j < v.NonZerosCount; j++)
                {
                    if (u._nonZeroIndices[i] == v._nonZeroIndices[j])
                    {
                        matrix.At(i, j, u._nonZeroValues[i] * v._nonZeroValues[j]);
                    }
                }
            }

            return matrix;
        }

        /// <summary>
        /// Returns the transpose of this matrix.
        /// </summary>        
        /// <returns>The transpose of this matrix.</returns>
        public override Matrix<Complex> Transpose()
        {
            var ret = new SparseMatrix(ColumnCount, RowCount)
            {
                _columnIndices = new int[NonZerosCount],
                _nonZeroValues = new Complex[NonZerosCount]
            };

            // Do an 'inverse' CopyTo iterate over the rows
            for (var i = 0; i < _rowIndex.Length; i++)
            {
                // Get the begin / end index for the current row
                var startIndex = _rowIndex[i];
                var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;

                // Get the values for the current row
                if (startIndex == endIndex)
                {
                    // Begin and end are equal. There are no values in the row, Move to the next row
                    continue;
                }

                for (var j = startIndex; j < endIndex; j++)
                {
                    ret.SetValueAt(_columnIndices[j], i, _nonZeroValues[j]);
                }
            }

            return ret;
        }