public void CanCheckRankOfSquareSingular(int order)
        {
            var matrixA = new DenseMatrix(order, order);
            matrixA[0, 0] = 1;
            matrixA[order - 1, order - 1] = 1;
            for (var i = 1; i < order - 1; i++)
            {
                matrixA[i, i - 1] = 1;
                matrixA[i, i + 1] = 1;
                matrixA[i - 1, i] = 1;
                matrixA[i + 1, i] = 1;
            }

            var factorSvd = matrixA.Svd();

            Assert.AreEqual(factorSvd.Determinant, Complex.Zero);
            Assert.AreEqual(factorSvd.Rank, order - 1);
        }

        public void CanCheckRankOfSquareSingular([Values(10, 50, 100)] int order)
        {
            var A = new DenseMatrix(order, order);
            A[0, 0] = 1;
            A[order - 1, order - 1] = 1;
            for (var i = 1; i < order - 1; i++)
            {
                A[i, i - 1] = 1;
                A[i, i + 1] = 1;
                A[i - 1, i] = 1;
                A[i + 1, i] = 1;
            }

            var factorEvd = A.Evd();

            Assert.AreEqual(factorEvd.Determinant, Complex.Zero);
            Assert.AreEqual(factorEvd.Rank, order - 1);
        }

 public void CalculateYBus()
 {
     Y_Bus = new DenseMatrix(buses.Count, buses.Count);
     foreach (Branch branch in branches)
     {
         int fromBus = Convert.ToInt32(branch.FromBus) - 1;
         int toBus = Convert.ToInt32(branch.ToBus) - 1;
         Complex impedance = new Complex(branch.Resistance,branch.Reactance);
         Complex admittance = 1/impedance;
         if (branch.Ratio != 0)
         {
             double t = 1 / branch.Ratio;
             Y_Bus[fromBus, toBus] = -1 * t * admittance + Y_Bus[fromBus, toBus];
             Y_Bus[toBus, fromBus] = -1 * t * admittance + Y_Bus[toBus, fromBus];
             Y_Bus[fromBus, fromBus] = Y_Bus[fromBus, fromBus] + t * t * admittance;
             Y_Bus[toBus, toBus] = Y_Bus[toBus, toBus] + admittance;
         }
         else
         {
             Complex shuntsusceptance = new Complex(0,0.5*branch.Susceptance);
             Y_Bus[fromBus, fromBus] = Y_Bus[fromBus, fromBus] + admittance + shuntsusceptance;
             Y_Bus[toBus, toBus] = Y_Bus[toBus, toBus] + admittance + shuntsusceptance;
             Y_Bus[fromBus, toBus] = Y_Bus[fromBus, toBus] - admittance;
             Y_Bus[toBus, fromBus] = Y_Bus[toBus, fromBus] - admittance;
         }
     }
     foreach (Bus bus in buses)
     {
         if (bus.ShuntSusceptance != 0)
         {
             Complex busshuntsusceptance = new Complex(0,bus.ShuntSusceptance);
             for (int i = 0; i < Y_Bus.ColumnCount; i++)
             {
                 if (i==(Convert.ToInt32(bus.ID)-1))
                 {
                     Y_Bus[i, i] = Y_Bus[i, i] + busshuntsusceptance;
                 }
             }
         }
     }
     //Console.WriteLine(Y_Bus);
 }

        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 new dense 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 DenseMatrix OfDiagonalVector(Vector<Complex> diagonal)
 {
     var m = new DenseMatrix(diagonal.Count, diagonal.Count);
     m.SetDiagonal(diagonal);
     return m;
 }

 /// <summary>
 /// Create a new dense 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 DenseMatrix OfDiagonalArray(int rows, int columns, Complex[] diagonal)
 {
     var m = new DenseMatrix(rows, columns);
     m.SetDiagonal(diagonal);
     return m;
 }

        public void CanComputeQRFactorWideMatrix()
        {
            var matrix = _matrices["Wide2x3"];
            var r = new Complex[matrix.RowCount*matrix.ColumnCount];
            Array.Copy(matrix.Values, r, r.Length);

            var tau = new Complex[3];
            var q = new Complex[matrix.RowCount*matrix.RowCount];
            Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau);

            var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
            var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
            var a = mq*mr;

            for (var row = 0; row < matrix.RowCount; row++)
            {
                for (var col = 0; col < matrix.ColumnCount; col++)
                {
                    AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
                }
            }
        }

        public void CanSolveForRandomMatrixWhenResultMatrixGivenUsingThinQR(int order)
        {
            var matrixA = Matrix<Complex>.Build.Random(order, order, 1);
            var matrixACopy = matrixA.Clone();
            var factorQR = matrixA.QR(QRMethod.Thin);

            var matrixB = Matrix<Complex>.Build.Random(order, order, 1);
            var matrixBCopy = matrixB.Clone();

            var matrixX = new DenseMatrix(order, order);
            factorQR.Solve(matrixB, matrixX);

            // The solution X row dimension is equal to the column dimension of A
            Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);

            // The solution X has the same number of columns as B
            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);

            var matrixBReconstruct = matrixA * matrixX;

            // Check the reconstruction.
            for (var i = 0; i < matrixB.RowCount; i++)
            {
                for (var j = 0; j < matrixB.ColumnCount; j++)
                {
                    AssertHelpers.AlmostEqual(matrixB[i, j], matrixBReconstruct[i, j], 10);
                }
            }

            // Make sure A didn't change.
            for (var i = 0; i < matrixA.RowCount; i++)
            {
                for (var j = 0; j < matrixA.ColumnCount; j++)
                {
                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
                }
            }

            // Make sure B didn't change.
            for (var i = 0; i < matrixB.RowCount; i++)
            {
                for (var j = 0; j < matrixB.ColumnCount; j++)
                {
                    Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]);
                }
            }
        }

        public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int column)
        {
            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column);
            var matrixACopy = matrixA.Clone();
            var factorSvd = matrixA.Svd();

            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(row, column);
            var matrixBCopy = matrixB.Clone();

            var matrixX = new DenseMatrix(column, column);
            factorSvd.Solve(matrixB, matrixX);

            // The solution X row dimension is equal to the column dimension of A
            Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);

            // The solution X has the same number of columns as B
            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);

            var matrixBReconstruct = matrixA*matrixX;

            // Check the reconstruction.
            for (var i = 0; i < matrixB.RowCount; i++)
            {
                for (var j = 0; j < matrixB.ColumnCount; j++)
                {
                    AssertHelpers.AlmostEqual(matrixB[i, j], matrixBReconstruct[i, j], 9);
                }
            }

            // Make sure A didn't change.
            for (var i = 0; i < matrixA.RowCount; i++)
            {
                for (var j = 0; j < matrixA.ColumnCount; j++)
                {
                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
                }
            }

            // Make sure B didn't change.
            for (var i = 0; i < matrixB.RowCount; i++)
            {
                for (var j = 0; j < matrixB.ColumnCount; j++)
                {
                    Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]);
                }
            }
        }

 public void CholeskyFailsWithNonSquareMatrix()
 {
     var matrix = new DenseMatrix(3, 2);
     Assert.That(() => matrix.Cholesky(), Throws.ArgumentException);
 }

        /// <summary>
        /// Returns the conjugate transpose of this matrix.
        /// </summary>        
        /// <returns>The conjugate transpose of this matrix.</returns>
        public override Matrix<Complex> ConjugateTranspose()
        {
            var ret = new DenseMatrix(_columnCount, _rowCount);
            for (var j = 0; j < _columnCount; j++)
            {
                var index = j * _rowCount;
                for (var i = 0; i < _rowCount; i++)
                {
                    ret._values[(i * _columnCount) + j] = _values[index + i].Conjugate();
                }
            }

            return ret;
        }

        /// <summary>
        /// Subtracts another matrix from this matrix.
        /// </summary>
        /// <param name="other">The matrix to subtract.</param>
        /// <returns>The result of the subtraction.</returns>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
        public override Matrix<Complex> Subtract(Matrix<Complex> other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (other.RowCount != RowCount || other.ColumnCount != ColumnCount)
            {
                throw DimensionsDontMatch<ArgumentOutOfRangeException>(this, other, "other");
            }

            Matrix<Complex> result;
            if (other is DiagonalMatrix)
            {
                result = new DenseMatrix(RowCount, ColumnCount);
            }
            else
            {
                result = new DiagonalMatrix(RowCount, ColumnCount);
            }

            Subtract(other, result);
            return result;
        }

        public void CanSolveUsingSVDSquareMatrix()
        {
            var matrix = _matrices["Square3x3"];
            var a = new Complex[matrix.RowCount*matrix.ColumnCount];
            Array.Copy(matrix.Values, a, a.Length);

            var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
            var x = new Complex[matrix.ColumnCount*2];
            Control.LinearAlgebraProvider.SvdSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x);

            NotModified(3, 3, a, matrix);

            var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
            var mb = matrix*mx;

            AssertHelpers.AlmostEqual(mb[0, 0], b[0], 13);
            AssertHelpers.AlmostEqual(mb[1, 0], b[1], 13);
            AssertHelpers.AlmostEqual(mb[2, 0], b[2], 13);
            AssertHelpers.AlmostEqual(mb[0, 1], b[3], 13);
            AssertHelpers.AlmostEqual(mb[1, 1], b[4], 13);
            AssertHelpers.AlmostEqual(mb[2, 1], b[5], 13);
        }

        public void CanComputeSVDFactorizationOfWideMatrix()
        {
            var matrix = _matrices["Wide2x3"];
            var a = new Complex[matrix.RowCount*matrix.ColumnCount];
            Array.Copy(matrix.Values, a, a.Length);

            var s = new Complex[matrix.RowCount];
            var u = new Complex[matrix.RowCount*matrix.RowCount];
            var vt = new Complex[matrix.ColumnCount*matrix.ColumnCount];

            Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt);

            var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
            for (var index = 0; index < s.Length; index++)
            {
                w[index, index] = s[index];
            }

            var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
            var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
            var result = mU*w*mV;

            AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 14);
            AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 14);
            AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 14);
            AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 14);
            AssertHelpers.AlmostEqualRelative(matrix[0, 2], result[0, 2], 14);
            AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 14);
        }

        public void CanSolveUsingThinQRTallMatrixOnFactoredMatrix()
        {
            var matrix = _matrices["Tall3x2"];
            var a = new Complex[matrix.RowCount*matrix.ColumnCount];
            Array.Copy(matrix.Values, a, a.Length);

            var tau = new Complex[matrix.ColumnCount];
            var r = new Complex[matrix.ColumnCount*matrix.ColumnCount];
            Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau);

            var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
            var x = new Complex[matrix.ColumnCount*2];
            Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, QRMethod.Thin);

            var mb = new DenseMatrix(matrix.RowCount, 2, b);
            var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;

            AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 13);
            AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 13);
            AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 13);
            AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
        }

        public void CanSolveUsingThinQRSquareMatrixOnFactoredMatrix()
        {
            var matrix = _matrices["Square3x3"];
            var a = new Complex[matrix.RowCount*matrix.ColumnCount];
            Array.Copy(matrix.Values, a, a.Length);

            var tau = new Complex[matrix.ColumnCount];
            var r = new Complex[matrix.ColumnCount*matrix.ColumnCount];
            Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau);

            var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
            var x = new Complex[matrix.ColumnCount*2];
            Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, QRMethod.Thin);

            var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
            var mb = matrix*mx;

            AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
            AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
            AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
            AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
            AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
            AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
        }

        public void CanComputeThinQRFactorTallMatrix()
        {
            var matrix = _matrices["Tall3x2"];
            var r = new Complex[matrix.ColumnCount*matrix.ColumnCount];
            var tau = new Complex[3];
            var q = new Complex[matrix.RowCount*matrix.ColumnCount];
            Array.Copy(matrix.Values, q, q.Length);

            Control.LinearAlgebraProvider.ThinQRFactor(q, matrix.RowCount, matrix.ColumnCount, r, tau);

            var mq = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, q);
            var mr = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, r);
            var a = mq*mr;

            for (var row = 0; row < matrix.RowCount; row++)
            {
                for (var col = 0; col < matrix.ColumnCount; col++)
                {
                    AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
                }
            }
        }

        /// <summary>
        /// Initializes a square <see cref="DenseMatrix"/> with all zero's except for ones on the diagonal.
        /// </summary>
        /// <param name="order">the size of the square matrix.</param>
        /// <returns>A dense identity matrix.</returns>
        /// <exception cref="ArgumentException">
        /// If <paramref name="order"/> is less than one.
        /// </exception>
        public static DenseMatrix Identity(int order)
        {
            var m = new DenseMatrix(order);
            for (var i = 0; i < order; i++)
            {
                m._values[(i * order) + i] = 1.0;
            }

            return m;
        }

        /// <summary>
        /// Adds another matrix to this matrix.
        /// </summary>
        /// <param name="other">The matrix to add to this matrix.</param>
        /// <returns>The result of the addition.</returns>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
        public override Matrix<Complex> Add(Matrix<Complex> other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (other.RowCount != RowCount || other.ColumnCount != ColumnCount)
            {
                throw new ArgumentOutOfRangeException("other", Resources.ArgumentMatrixDimensions);
            }

            Matrix<Complex> result;
            if (other is DiagonalMatrix)
            {
                result = new DenseMatrix(RowCount, ColumnCount);
            }
            else
            {
                result = new DiagonalMatrix(RowCount, ColumnCount);
            }

            Add(other, result);
            return result;
        }

        public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int col)
        {
            var matrixA = Matrix<Complex>.Build.RandomPositiveDefinite(row, 1);
            var matrixACopy = matrixA.Clone();
            var chol = matrixA.Cholesky();
            var matrixB = Matrix<Complex>.Build.Random(row, col, 1);
            var matrixBCopy = matrixB.Clone();
            var matrixX = new DenseMatrix(row, col);
            chol.Solve(matrixB, matrixX);

            Assert.AreEqual(matrixB.RowCount, matrixX.RowCount);
            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);

            var matrixBReconstruct = matrixA * matrixX;

            // Check the reconstruction.
            for (var i = 0; i < matrixB.RowCount; i++)
            {
                for (var j = 0; j < matrixB.ColumnCount; j++)
                {
                    AssertHelpers.AlmostEqual(matrixB[i, j], matrixBReconstruct[i, j], 10);
                }
            }

            // Make sure A didn't change.
            for (var i = 0; i < matrixA.RowCount; i++)
            {
                for (var j = 0; j < matrixA.ColumnCount; j++)
                {
                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
                }
            }

            // Make sure B didn't change.
            for (var i = 0; i < matrixB.RowCount; i++)
            {
                for (var j = 0; j < matrixB.ColumnCount; j++)
                {
                    Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]);
                }
            }
        }