public void Test1()
        {
            MatrixF a = new MatrixF(new float[,] { { 2, -1, 0},
                                             { -1, 2, -1},
                                             { 0, -1, 2} });

              CholeskyDecompositionF d = new CholeskyDecompositionF(a);

              Assert.AreEqual(true, d.IsSymmetricPositiveDefinite);

              MatrixF l = d.L;

              Assert.AreEqual(0, l[0, 1]);
              Assert.AreEqual(0, l[0, 2]);
              Assert.AreEqual(0, l[1, 2]);

              Assert.IsTrue(MatrixF.AreNumericallyEqual(a, l * l.Transposed));

              Assert.IsTrue(MatrixF.AreNumericallyEqual(a, l * l.Transposed));

              // Check solving of linear equations.
              MatrixF x = new MatrixF(new float[,] { { 1, 2},
                                             { 3, 4},
                                             { 5, 6} });
              MatrixF b = a * x;

              Assert.IsTrue(MatrixF.AreNumericallyEqual(x, d.SolveLinearEquations(b)));
        }

        public void DeterminantException()
        {
            MatrixF a = new MatrixF(new float[,] { { 1, 2 }, { 5, 6 }, { 0, 1 } });

              LUDecompositionF d = new LUDecompositionF(a);
              float det = d.Determinant;
        }

 public void ConstructorException2()
 {
     MatrixF m = new MatrixF(new float[,] {{ 1, 2 },
                                    { 3, 4 },
                                    { 5, 6 }});
       new EigenvalueDecompositionF(m);
 }

        public void TestRandomA()
        {
            RandomHelper.Random = new Random(1);

              for (int i = 0; i < 100; i++)
              {
            // Create A.
            MatrixF a = new MatrixF(3, 3);
            RandomHelper.Random.NextMatrixF(a, 0, 1);

            LUDecompositionF d = new LUDecompositionF(a);

            if (d.IsNumericallySingular == false)
            {
              // Check solving of linear equations.
              MatrixF b = new MatrixF(3, 2);
              RandomHelper.Random.NextMatrixF(b, 0, 1);

              MatrixF x = d.SolveLinearEquations(b);
              MatrixF b2 = a * x;
              Assert.IsTrue(MatrixF.AreNumericallyEqual(b, b2, 0.01f));

              MatrixF aPermuted = d.L * d.U;
              Assert.IsTrue(MatrixF.AreNumericallyEqual(aPermuted, a.GetSubmatrix(d.PivotPermutationVector, 0, 2)));
            }
              }
        }

        //--------------------------------------------------------------
        /// <summary>
        /// Creates the principal component analysis for the given list of points.
        /// </summary>
        /// <param name="points">
        /// The list of data points. All points must have the same 
        /// <see cref="VectorF.NumberOfElements"/>.
        /// </param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="points"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// <paramref name="points"/> is empty.
        /// </exception>
        public PrincipalComponentAnalysisF(IList<VectorF> points)
        {
            if (points == null)
            throw new ArgumentNullException("points");
              if (points.Count == 0)
            throw new ArgumentException("The list of points is empty.");

              // Compute covariance matrix.
              MatrixF covarianceMatrix = StatisticsHelper.ComputeCovarianceMatrix(points);

              // Perform Eigenvalue decomposition.
              EigenvalueDecompositionF evd = new EigenvalueDecompositionF(covarianceMatrix);

              int numberOfElements = evd.RealEigenvalues.NumberOfElements;
              Variances = new VectorF(numberOfElements);
              V = new MatrixF(numberOfElements, numberOfElements);

              // Sort eigenvalues by decreasing value.
              // Since covarianceMatrix is symmetric, we have no imaginary eigenvalues.
              for (int i = 0; i < Variances.NumberOfElements; i++)
              {
            int index = evd.RealEigenvalues.IndexOfLargestElement;

            Variances[i] = evd.RealEigenvalues[index];
            V.SetColumn(i, evd.V.GetColumn(index));

            evd.RealEigenvalues[index] = float.NegativeInfinity;
              }
        }

        public void TestMatricesWithoutFullRank()
        {
            MatrixF a = new MatrixF(3, 3);
              SingularValueDecompositionF svd = new SingularValueDecompositionF(a);
              Assert.AreEqual(0, svd.NumericalRank);
              Assert.AreEqual(svd.SingularValues[0], svd.Norm2);
              float condNumber = svd.ConditionNumber;

              a = new MatrixF(new float[,] { { 1, 2, 3 }, { 4, 5, 6 }, { 4, 5, 6 } });
              svd = new SingularValueDecompositionF(a);
              Assert.AreEqual(2, svd.NumericalRank);
              Assert.AreEqual(svd.SingularValues[0], svd.Norm2);
              Assert.IsTrue(MatrixF.AreNumericallyEqual(a, svd.U * svd.S * svd.V.Transposed));
              svd = new SingularValueDecompositionF(a.Transposed);
              Assert.AreEqual(2, svd.NumericalRank);
              Assert.IsTrue(MatrixF.AreNumericallyEqual(a.Transposed, svd.U * svd.S * svd.V.Transposed));
              Assert.IsTrue(MatrixF.AreNumericallyEqual(a.Transposed, svd.U * svd.S * svd.V.Transposed)); // Repeat to test with cached values.
              Assert.AreEqual(svd.SingularValues[0], svd.Norm2);
              condNumber = svd.ConditionNumber;

              a = new MatrixF(new float[,] { { 1, 2 }, { 1, 2 }, { 1, 2 } });
              svd = new SingularValueDecompositionF(a);
              Assert.AreEqual(1, svd.NumericalRank);
              Assert.IsTrue(MatrixF.AreNumericallyEqual(a, svd.U * svd.S * svd.V.Transposed));
              Assert.AreEqual(svd.SingularValues[0], svd.Norm2);
              condNumber = svd.ConditionNumber;
        }

 public void AddMatrix()
 {
     MatrixF m1 = new MatrixF(3, 4, rowMajor, MatrixOrder.RowMajor);
       MatrixF m2 = new MatrixF(3, 4, rowMajor, MatrixOrder.RowMajor) * (-3);
       MatrixF result = MatrixF.Add(m1, m2);
       for (int i = 0; i < 12; i++)
     Assert.AreEqual(-rowMajor[i] * 2, result[i]);
 }

        public void SolveLinearEquationsException1()
        {
            // Create A.
              MatrixF a = new MatrixF(new float[,] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, -9 } });

              CholeskyDecompositionF decomp = new CholeskyDecompositionF(a);
              decomp.SolveLinearEquations(null);
        }

        public void Test1()
        {
            MatrixF a = new MatrixF(new float[,] {{ 1, -1, 4 },
                                           { 3, 2, -1 },
                                           { 2, 1, -1}});
              EigenvalueDecompositionF d = new EigenvalueDecompositionF(a);

              Assert.IsTrue(MatrixF.AreNumericallyEqual(a * d.V, d.V * d.D));
        }

        public void Test2()
        {
            MatrixF a = new MatrixF(new float[,] {{ 0, 1, 2 },
                                           { 1, 4, 3 },
                                           { 2, 3, 5}});
              EigenvalueDecompositionF d = new EigenvalueDecompositionF(a);

              Assert.IsTrue(MatrixF.AreNumericallyEqual(a, d.V * d.D * d.V.Transposed));
        }

        public void SolveLinearEquationsException3()
        {
            // Create A.
              MatrixF a = new MatrixF(new float[,] { { 1, 2, 3 }, { 4, 5, 6 }, { 1, 2, 3 } });
              MatrixF b = new MatrixF(new float[,] { { 1, 2, 3 }, { 4, 5, 6 }, {7, 8, -9 } });

              QRDecompositionF decomp = new QRDecompositionF(a);
              decomp.SolveLinearEquations(b);
        }

        //--------------------------------------------------------------
        /// <summary>
        /// Creates the eigenvalue decomposition of the given matrix.
        /// </summary>
        /// <param name="matrixA">The square matrix A.</param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="matrixA"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// <paramref name="matrixA"/> is non-square (rectangular).
        /// </exception>
        public EigenvalueDecompositionF(MatrixF matrixA)
        {
            if (matrixA == null)
            throw new ArgumentNullException("matrixA");
              if (matrixA.IsSquare == false)
            throw new ArgumentException("The matrix A must be square.", "matrixA");

              _n = matrixA.NumberOfColumns;
              _d = new VectorF(_n);
              _e = new VectorF(_n);

              _isSymmetric = matrixA.IsSymmetric;

              if (_isSymmetric)
              {
            _v = matrixA.Clone();

            // Tridiagonalize.
            ReduceToTridiagonal();

            // Diagonalize.
            TridiagonalToQL();
              }
              else
              {
            _v = new MatrixF(_n, _n);

            // Abort if A contains NaN values.
            // If we continue with NaN values, we run into an infinite loop.
            for (int i = 0; i < _n; i++)
            {
              for (int j = 0; j < _n; j++)
              {
            if (Numeric.IsNaN(matrixA[i, j]))
            {
              _e.Set(float.NaN);
              _v.Set(float.NaN);
              _d.Set(float.NaN);
              return;
            }
              }
            }

            // Storage of nonsymmetric Hessenberg form.
            MatrixF matrixH = matrixA.Clone();
            // Working storage for nonsymmetric algorithm.
            float[] ort = new float[_n];

            // Reduce to Hessenberg form.
            ReduceToHessenberg(matrixH, ort);

            // Reduce Hessenberg to real Schur form.
            HessenbergToRealSchur(matrixH);
              }
        }

        public void Determinant()
        {
            MatrixF a = new MatrixF(new float[,] { { 1, 2, 3, 4 }, { 5, 6, 7, 8 }, { 0, 1, 2, 0 }, { 1, 0, 1, 0 } });

              LUDecompositionF d = new LUDecompositionF(a);
              Assert.AreEqual(false, d.IsNumericallySingular);
              Assert.IsTrue(Numeric.AreEqual(-24, d.Determinant));

              MatrixF aPermuted = d.L * d.U;
              Assert.IsTrue(MatrixF.AreNumericallyEqual(aPermuted, a.GetSubmatrix(d.PivotPermutationVector, 0, 3)));
        }

        public void Test1()
        {
            MatrixF A = new MatrixF(new float[,] { { 4 } });
              VectorF b = new VectorF(new float[] { 20 });

              SorMethodF solver = new SorMethodF();
              VectorF x = solver.Solve(A, null, b);

              Assert.IsTrue(VectorF.AreNumericallyEqual(new VectorF(1, 5), x));
              Assert.AreEqual(2, solver.NumberOfIterations);
        }

        public void SolveWithDefaultInitialGuess()
        {
            MatrixF A = new MatrixF(new float[,] { { 4 } });
              VectorF b = new VectorF(new float[] { 20 });

              JacobiMethodF solver = new JacobiMethodF();
              VectorF x = solver.Solve(A, b);

              Assert.IsTrue(VectorF.AreNumericallyEqual(new VectorF(1, 5), x));
              Assert.AreEqual(2, solver.NumberOfIterations);
        }

        public void Test3()
        {
            MatrixF A = new MatrixF(new float[,] { { 2, 0 },
                                             { 0, 2 }});
              VectorF b = new VectorF(new float[] { 20, 28 });

              GaussSeidelMethodF solver = new GaussSeidelMethodF();
              VectorF x = solver.Solve(A, null, b);

              Assert.IsTrue(VectorF.AreNumericallyEqual(b / 2, x));
              Assert.AreEqual(2, solver.NumberOfIterations);
        }

        public void Test2()
        {
            MatrixF A = new MatrixF(new float[,] { { 1, 0 },
                                             { 0, 1 }});
              VectorF b = new VectorF(new float[] { 20, 28 });

              JacobiMethodF solver = new JacobiMethodF();
              VectorF x = solver.Solve(A, null, b);

              Assert.IsTrue(VectorF.AreNumericallyEqual(b, x));
              Assert.AreEqual(2, solver.NumberOfIterations);
        }

        public void Test4()
        {
            MatrixF A = new MatrixF(new float[,] { { -12, 2 },
                                             { 2, 3 }});
              VectorF b = new VectorF(new float[] { 20, 28 });

              SorMethodF solver = new SorMethodF();
              VectorF x = solver.Solve(A, null, b);

              VectorF solution = MatrixF.SolveLinearEquations(A, b);
              Assert.IsTrue(VectorF.AreNumericallyEqual(solution, x));
        }

        public void Test5()
        {
            MatrixF A = new MatrixF(new float[,] { { -21, 2, -4, 0 },
                                             { 2, 3, 0.1f, -1 },
                                             { 2, 10, 111.1f, -11 },
                                             { 23, 112, 111.1f, -143 }});
              VectorF b = new VectorF(new float[] { 20, 28, -12, 0.1f });

              GaussSeidelMethodF solver = new GaussSeidelMethodF();
              VectorF x = solver.Solve(A, null, b);

              VectorF solution = MatrixF.SolveLinearEquations(A, b);
              Assert.IsTrue(VectorF.AreNumericallyEqual(solution, x));
        }

        /// <summary>
        /// Solves the specified linear system of equations <i>Ax=b</i>.
        /// </summary>
        /// <param name="matrixA">The matrix A.</param>
        /// <param name="initialX">
        /// The initial guess for x. If this value is <see langword="null"/>, a zero vector will be used
        /// as initial guess.
        /// </param>
        /// <param name="vectorB">The vector b.</param>
        /// <returns>The solution vector x.</returns>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="matrixA"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="vectorB"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// <paramref name="matrixA"/> is not a square matrix.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// The number of elements of <paramref name="initialX"/> does not match.
        /// </exception>
        public override VectorF Solve(MatrixF matrixA, VectorF initialX, VectorF vectorB)
        {
            // TODO: We can possible improve the method by reordering after each step.
              // This can be done randomly or we sort by the "convergence" of the elements.
              // See book Physics-Based Animation.

              NumberOfIterations = 0;

              if (matrixA == null)
            throw new ArgumentNullException("matrixA");
              if (vectorB == null)
            throw new ArgumentNullException("vectorB");
              if (matrixA.IsSquare == false)
            throw new ArgumentException("Matrix A must be a square matrix.", "matrixA");
              if (matrixA.NumberOfRows != vectorB.NumberOfElements)
            throw new ArgumentException("The number of rows of A and b do not match.");
              if (initialX != null && initialX.NumberOfElements != vectorB.NumberOfElements)
            throw new ArgumentException("The number of elements of the initial guess for x and b do not match.");

              VectorF xOld = initialX ?? new VectorF(vectorB.NumberOfElements);
              VectorF xNew = new VectorF(vectorB.NumberOfElements);
              bool isConverged = false;

              // Make iterations until max iteration count or the result has converged.
              for (int i = 0; i < MaxNumberOfIterations && !isConverged; i++)
              {
            for (int j = 0; j < vectorB.NumberOfElements; j++)
            {
              float delta = 0;
              for (int k = 0; k < j; k++)
            delta += matrixA[j, k] * xOld[k];

              for (int k = j + 1; k < vectorB.NumberOfElements; k++)
            delta += matrixA[j, k] * xOld[k];

              xNew[j] = (vectorB[j] - delta) / matrixA[j, j];
            }

            // Test convergence
            isConverged = VectorF.AreNumericallyEqual(xOld, xNew, Epsilon);

            xOld = xNew.Clone();
            NumberOfIterations = i + 1;
              }

              return xNew;
        }