C# MatrixF.Clone方法代码示例

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

        //--------------------------------------------------------------
        /// <summary>
        /// Creates the QR decomposition of the given matrix.
        /// </summary>
        /// <param name="matrixA">
        /// The matrix A. (Can be rectangular. NumberOfRows must be ≥ NumberOfColumns.)
        /// </param>
        /// <remarks>
        /// The QR decomposition is computed by Householder reflections.
        /// </remarks>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="matrixA"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// The number of rows must be greater than or equal to the number of columns.
        /// </exception>
        public QRDecompositionF(MatrixF matrixA)
        {
            if (matrixA == null)
            throw new ArgumentNullException("matrixA");
              if (matrixA.NumberOfRows < matrixA.NumberOfColumns)
            throw new ArgumentException("The number of rows must be greater than or equal to the number of columns.", "matrixA");

              // Initialize.
              _qr = matrixA.Clone();
              _m = matrixA.NumberOfRows;
              _n = matrixA.NumberOfColumns;
              _rDiagonal = new float[_n];

              // Main loop.
              for (int k = 0; k < _n; k++)
              {
            // Compute 2-norm of k-th column without under/overflow.
            float norm = 0;
            for (int i = k; i < _m; i++)
              norm = MathHelper.Hypotenuse(norm, _qr[i, k]);

            if (norm != 0)   // TODO: Maybe a comparison with an epsilon tolerance is required here!?
            {
              // Form k-th Householder vector.
              if (_qr[k, k] < 0)
            norm = -norm;
              for (int i = k; i < _m; i++)
            _qr[i, k] /= norm;
              _qr[k, k] += 1;

              // Apply transformation to remaining columns.
              for (int j = k + 1; j < _n; j++)
              {
            float s = 0;
            for (int i = k; i < _m; i++)
              s += _qr[i, k] * _qr[i, j];
            s = -s / _qr[k, k];
            for (int i = k; i < _m; i++)
              _qr[i, j] += s * _qr[i, k];
              }
            }
            _rDiagonal[k] = -norm;
              }
        }

        public void Absolute()
        {
            float[] values = new float[] { -1.0f, -2.0f, -3.0f,
                                     -4.0f, -5.0f, -6.0f,
                                     -7.0f, -8.0f, -9.0f };
              MatrixF m = new MatrixF(3, 3, values, MatrixOrder.RowMajor);

              MatrixF absolute = m.Clone();
              absolute.Absolute();
              for (int i = 0; i < absolute.NumberOfRows; i++)
            for (int j = 0; j < absolute.NumberOfColumns; j++)
              Assert.AreEqual(i * absolute.NumberOfColumns + j + 1, absolute[i, j]);

              absolute = MatrixF.Absolute(m);
              for (int i = 0; i < absolute.NumberOfRows; i++)
            for (int j = 0; j < absolute.NumberOfColumns; j++)
              Assert.AreEqual(i * absolute.NumberOfColumns + j + 1, absolute[i, j]);

              values = new float[] { 1.0f, 2.0f, 3.0f,
                             4.0f, 5.0f, 6.0f,
                             7.0f, 8.0f, 9.0f };
              m = new MatrixF(3, 3, values, MatrixOrder.RowMajor);

              absolute = m.Clone();
              absolute.Absolute();
              for (int i = 0; i < absolute.NumberOfRows; i++)
            for (int j = 0; j < absolute.NumberOfColumns; j++)
              Assert.AreEqual(i * absolute.NumberOfColumns + j + 1, absolute[i, j]);

              absolute = MatrixF.Absolute(m);
              for (int i = 0; i < absolute.NumberOfRows; i++)
            for (int j = 0; j < absolute.NumberOfColumns; j++)
              Assert.AreEqual(i * absolute.NumberOfColumns + j + 1, absolute[i, j]);

              Assert.IsNull(MatrixF.Absolute(null));
        }

        //--------------------------------------------------------------
        /// <summary>
        /// Returns the least squares solution for the equation <c>A * X = B</c>.
        /// </summary>
        /// <param name="matrixB">The matrix B with as many rows as A and any number of columns.</param>
        /// <returns>X with the least squares solution.</returns>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="matrixB"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// The number of rows does not match.
        /// </exception>
        /// <exception cref="MathematicsException">
        /// The matrix A does not have full rank.
        /// </exception>
        public MatrixF SolveLinearEquations(MatrixF matrixB)
        {
            if (matrixB == null)
            throw new ArgumentNullException("matrixB");
              if (matrixB.NumberOfRows != _m)
            throw new ArgumentException("The number of rows does not match.", "matrixB");
              if (HasNumericallyFullRank == false)
            throw new MathematicsException("The matrix does not have full rank.");

              // Copy right hand side
              MatrixF x = matrixB.Clone();

              // Compute Y = transpose(Q)*B
              for (int k = 0; k < _n; k++)
              {
            for (int j = 0; j < matrixB.NumberOfColumns; j++)
            {
              float s = 0;
              for (int i = k; i < _m; i++)
            s += _qr[i, k] * x[i, j];
              s = -s / _qr[k, k];
              for (int i = k; i < _m; i++)
            x[i, j] += s * _qr[i, k];
            }
              }
              // Solve R*X = Y.
              for (int k = _n - 1; k >= 0; k--)
              {
            for (int j = 0; j < matrixB.NumberOfColumns; j++)
              x[k, j] /= _rDiagonal[k];
            for (int i = 0; i < k; i++)
              for (int j = 0; j < matrixB.NumberOfColumns; j++)
            x[i, j] -= x[k, j] * _qr[i, k];
              }
              return x.GetSubmatrix(0, _n - 1, 0, matrixB.NumberOfColumns - 1);
        }

        //--------------------------------------------------------------
        /// <summary>
        /// Solves the equation <c>A * X = B</c>.
        /// </summary>
        /// <param name="matrixB">The matrix B with as many rows as A and any number of columns.</param>
        /// <returns>X, so that <c>A * X = B</c>.</returns>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="matrixB"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// The number of rows does not match.
        /// </exception>
        /// <exception cref="MathematicsException">
        /// The matrix A is not symmetric and positive definite.
        /// </exception>
        public MatrixF SolveLinearEquations(MatrixF matrixB)
        {
            if (matrixB == null)
            throw new ArgumentNullException("matrixB");
              if (matrixB.NumberOfRows != L.NumberOfRows)
            throw new ArgumentException("The number of rows does not match.", "matrixB");
              if (IsSymmetricPositiveDefinite == false)
            throw new MathematicsException("The original matrix A is not symmetric and positive definite.");

              // Initialize x as a copy of B.
              MatrixF x = matrixB.Clone();

              // Solve L*Y = B.
              for (int k = 0; k < L.NumberOfRows; k++)
              {
            for (int j = 0; j < matrixB.NumberOfColumns; j++)
            {
              for (int i = 0; i < k; i++)
            x[k, j] -= x[i, j] * L[k, i];
              x[k, j] /= L[k, k];
            }
              }

              // Solve transpose(L) * X = Y.
              for (int k = L.NumberOfRows - 1; k >= 0; k--)
              {
            for (int j = 0; j < matrixB.NumberOfColumns; j++)
            {
              for (int i = k + 1; i < L.NumberOfRows; i++)
            x[k, j] -= x[i, j] * L[i, k];
              x[k, j] /= L[k, k];
            }
              }
              return x;
        }

        public void Invert()
        {
            Assert.AreEqual(MatrixF.CreateIdentity(3, 3), MatrixF.CreateIdentity(3, 3).Inverse);

              MatrixF m = new MatrixF(new float[,] {{1, 2, 3, 4},
                                            {2, 5, 8, 3},
                                            {7, 6, -1, 1},
                                            {4, 9, 7, 7}});
              MatrixF inverse = m.Clone();
              m.Invert();
              VectorF v = new VectorF(4, 1);
              VectorF w = m * v;
              Assert.IsTrue(VectorF.AreNumericallyEqual(v, inverse * w));
              Assert.IsTrue(MatrixF.AreNumericallyEqual(MatrixF.CreateIdentity(4, 4), m * inverse));

              m = new MatrixF(new float[,] {{1, 2, 3},
                                    {2, 5, 8},
                                    {7, 6, -1},
                                    {4, 9, 7}});
              // To check the pseudo-inverse we use the definition: A*A.Transposed*A = A
              // see http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse
              inverse = m.Clone();
              inverse.Invert();
              Assert.IsTrue(MatrixF.AreNumericallyEqual(m, m * inverse * m));
        }

 public void Clone()
 {
     MatrixF m = new MatrixF(3, 4, rowMajor, MatrixOrder.RowMajor);
       var o = m.Clone();
       Assert.AreEqual(m, o);
 }

        public void TryInvert()
        {
            // Regular, square
              MatrixF m = new MatrixF(new float[,] {{1, 2, 3, 4},
                                            {2, 5, 8, 3},
                                            {7, 6, -1, 1},
                                            {4, 9, 7, 7}});
              MatrixF inverse = m.Clone();
              Assert.AreEqual(true, m.TryInvert());
              Assert.IsTrue(MatrixF.AreNumericallyEqual(MatrixF.CreateIdentity(4, 4), m * inverse));

              // Full column rank, rectangular
              m = new MatrixF(new float[,] {{1, 2, 3},
                                    {2, 5, 8},
                                    {7, 6, -1},
                                    {4, 9, 7}});
              inverse = m.Clone();
              Assert.AreEqual(true, m.TryInvert());
              Assert.IsTrue(MatrixF.AreNumericallyEqual(m, m * inverse * m));

              // singular
              m = new MatrixF(new float[,] {{1, 2, 3},
                                    {2, 5, 8},
                                    {3, 7, 11}});
              inverse = m.Clone();
              Assert.AreEqual(false, m.TryInvert());
        }

        public SingularValueDecompositionF(MatrixF matrixA)
        {
            if (matrixA == null)
            throw new ArgumentNullException("matrixA");

              // Derived from LINPACK code.
              // Initialize.
              _m = matrixA.NumberOfRows;
              _n = matrixA.NumberOfColumns;
              MatrixF matrixAClone = matrixA.Clone();

              if (_m < _n)
            throw new ArgumentException("The number of rows must be greater than or equal to the number of columns.", "matrixA");

              int nu = Math.Min(_m, _n);
              _s = new VectorF(Math.Min(_m + 1, _n));
              _u = new MatrixF(_m, nu);     //Jama getU() returns new Matrix(U,_m,Math.min(_m+1,_n)) ?!
              _v = new MatrixF(_n, _n);
              float[] e = new float[_n];
              float[] work = new float[_m];

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

              // By default, we calculate U and V. To calculate only U or V we can set one of the following
              // two constants to false. (This optimization is not yet tested.)
              const bool wantu = true;
              const bool wantv = true;

              // Reduce A to bidiagonal form, storing the diagonal elements
              // in s and the super-diagonal elements in e.

              int nct = Math.Min(_m - 1, _n);
              int nrt = Math.Max(0, Math.Min(_n - 2, _m));
              for (int k = 0; k < Math.Max(nct, nrt); k++)
              {
            if (k < nct)
            {
              // Compute the transformation for the k-th column and
              // place the k-th diagonal in s[k].
              // Compute 2-norm of k-th column without under/overflow.
              _s[k] = 0;
              for (int i = k; i < _m; i++)
            _s[k] = MathHelper.Hypotenuse(_s[k], matrixAClone[i, k]);

              if (_s[k] != 0)
              {
            if (matrixAClone[k, k] < 0)
              _s[k] = -_s[k];

            for (int i = k; i < _m; i++)
              matrixAClone[i, k] /= _s[k];

            matrixAClone[k, k] += 1;
              }

              _s[k] = -_s[k];
            }
            for (int j = k + 1; j < _n; j++)
            {
              if ((k < nct) && (_s[k] != 0))
              {
            // Apply the transformation.
            float t = 0;
            for (int i = k; i < _m; i++)
              t += matrixAClone[i, k] * matrixAClone[i, j];

            t = -t / matrixAClone[k, k];
            for (int i = k; i < _m; i++)
              matrixAClone[i, j] += t * matrixAClone[i, k];
              }

              // Place the k-th row of A into e for the
              // subsequent calculation of the row transformation.

              e[j] = matrixAClone[k, j];
            }

            if (wantu & (k < nct))
            {
              // Place the transformation in U for subsequent back
              // multiplication.
              for (int i = k; i < _m; i++)
            _u[i, k] = matrixAClone[i, k];
            }

            if (k < nrt)
//.........这里部分代码省略.........

        //--------------------------------------------------------------
        /// <summary>
        /// Creates the LU decomposition of the given matrix.
        /// </summary>
        /// <param name="matrixA">
        /// The matrix A. (Can be rectangular. Number of rows ≥ number of columns.)
        /// </param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="matrixA"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// The number of rows must be greater than or equal to the number of columns.
        /// </exception>
        public LUDecompositionF(MatrixF matrixA)
        {
            if (matrixA == null)
            throw new ArgumentNullException("matrixA");
              if (matrixA.NumberOfColumns > matrixA.NumberOfRows)
            throw new ArgumentException("The number of rows must be greater than or equal to the number of columns.", "matrixA");

              // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
              _lu = matrixA.Clone();
              _m = matrixA.NumberOfRows;
              _n = matrixA.NumberOfColumns;
              _pivotVector = new int[_m];

              // Initialize with the 0 to m-1.
              for (int i = 0; i < _m; i++)
            _pivotVector[i] = i;

              _pivotSign = 1;

              // Outer loop.
              for (int j = 0; j < _n; j++)
              {
            // Make a copy of the j-th column to localize references.
            float[] luColumnJ = new float[_m];
            for (int i = 0; i < _m; i++)
              luColumnJ[i] = _lu[i, j];

            // Apply previous transformations.
            for (int i = 0; i < _m; i++)
            {
              // Most of the time is spent in the following dot product.
              int kmax = Math.Min(i, j);
              float s = 0;
              for (int k = 0; k < kmax; k++)
            s += _lu[i, k] * luColumnJ[k];

              luColumnJ[i] -= s;
              _lu[i, j] = luColumnJ[i];
            }

            // Find pivot and exchange if necessary.
            int p = j;
            for (int i = j + 1; i < _m; i++)
              if (Math.Abs(luColumnJ[i]) > Math.Abs(luColumnJ[p]))
            p = i;

            // Swap lines p and k.
            if (p != j)
            {
              for (int k = 0; k < _n; k++)
              {
            float dummy = _lu[p, k];
            _lu[p, k] = _lu[j, k];
            _lu[j, k] = dummy;
              }

              MathHelper.Swap(ref _pivotVector[p], ref _pivotVector[j]);

              _pivotSign = -_pivotSign;
            }

            // Compute multipliers.
            if (j < _m && _lu[j, j] != 0)
              for (int i = j + 1; i < _m; i++)
            _lu[i, j] /= _lu[j, j];
              }
        }