C# MatrixD.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 EigenvalueDecompositionD(MatrixD 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 VectorD(_n);
              _e = new VectorD(_n);

              _isSymmetric = matrixA.IsSymmetric;

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

            // Tridiagonalize.
            ReduceToTridiagonal();

            // Diagonalize.
            TridiagonalToQL();
              }
              else
              {
            _v = new MatrixD(_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(double.NaN);
              _v.Set(double.NaN);
              _d.Set(double.NaN);
              return;
            }
              }
            }

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

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

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

		/// <summary>
		/// 部分ピボット選択付きガウス消去法により連立一次方程式の解を求める
		/// </summary>
		/// <param name="a"></param>
		/// <param name="b"></param>
		/// <returns></returns>
		private static VectorD SolveByGEPP(MatrixD a, VectorD b)
		{
			int r = 0;
			int c = 0;
			var a2 = a.Clone();
			var b2 = b.Clone();

			while (r < a2.Rows - 1)
			{
				var maxA = a2.GetColumn(c).Where((x, i) => i >= r).Max();
				var idxMaxA = a2.GetRowIndex((x, i) => (i >= r) && (x[c] == maxA));
				if (idxMaxA != r)
				{
					a2.SwapRow(r, idxMaxA);
					b2.Swap(r, idxMaxA);
				}

				for (int i = r + 1; i < a2.Rows; i++)
				{
					var x = a2[i, c] * -1.0;
					a2[i, c] = 0;
					for (int j = c + 1; j < a2.Columns; j++)
					{
						a2[i, j] += a2[r, j] * (x / a2[r, c]);
					}

					b2[i] += b2[r] * (x / a2[r, c]);
				}

				r++;
				c++;
			}

			var n = b2.Length;
			var result = new VectorD(new double[n]);
			for (int k = n - 1; k >= 0; k--)
			{
				var sum = 0.0;
				for (int j = k + 1; j < n; j++)
				{
					sum += a2[k, j] * result[j];
				}
				result[k] = (b2[k] - sum) / a2[k, k];
			}

			return result;
		}

        //--------------------------------------------------------------
        /// <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 QRDecompositionD(MatrixD 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 double[_n];

              // Main loop.
              for (int k = 0; k < _n; k++)
              {
            // Compute 2-norm of k-th column without under/overflow.
            double 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++)
              {
            double 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;
              }
        }

		/// <summary>
		/// ガウス消去法により連立一次方程式の解を求める
		/// </summary>
		/// <param name="a"></param>
		/// <param name="b"></param>
		/// <returns></returns>
		private static VectorD SolveByGE(MatrixD a, VectorD b)
		{
			int r = 0;
			int c = 0;
			var a2 = a.Clone();
			var b2 = b.Clone();

			while (r < a2.Rows - 1)
			{
				for (int i = r + 1; i < a2.Rows; i++)
				{
					var x = a2[i, c] * -1.0;
					a2[i, c] = 0;
					for (int j = c + 1; j < a2.Columns; j++)
					{
						a2[i, j] += a2[r, j] * (x / a2[r, c]);
					}

					b2[i] += b2[r] * (x / a2[r, c]);
				}

				r++;
				c++;
			}

			var n = b2.Length;
			var result = new VectorD(new double[n]);
			for (int k = n - 1; k >= 0; k--)
			{
				var sum = 0.0;
				for (int j = k + 1; j < n; j++)
				{
					sum += a2[k, j] * result[j];
				}
				result[k] = (b2[k] - sum) / a2[k, k];
			}

			return result;
		}

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

              MatrixD 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 = MatrixD.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 double[] { 1.0, 2.0, 3.0,
                             4.0, 5.0, 6.0,
                             7.0, 8.0, 9.0 };
              m = new MatrixD(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 = MatrixD.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(MatrixD.Absolute(null));
        }

        //--------------------------------------------------------------
        /// <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 LUDecompositionD(MatrixD 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.
            double[] luColumnJ = new double[_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);
              double 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++)
              {
            double 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];
              }
        }

        //--------------------------------------------------------------
        /// <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 MatrixD SolveLinearEquations(MatrixD 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
              MatrixD x = matrixB.Clone();

              // Compute Y = transpose(Q)*B
              for (int k = 0; k < _n; k++)
              {
            for (int j = 0; j < matrixB.NumberOfColumns; j++)
            {
              double 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);
        }

        public SingularValueDecompositionD(MatrixD matrixA)
        {
            if (matrixA == null)
            throw new ArgumentNullException("matrixA");

              // Derived from LINPACK code.
              // Initialize.
              _m = matrixA.NumberOfRows;
              _n = matrixA.NumberOfColumns;
              MatrixD 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 VectorD(Math.Min(_m + 1, _n));
              _u = new MatrixD(_m, nu);     //Jama getU() returns new Matrix(U,_m,Math.min(_m+1,_n)) ?!
              _v = new MatrixD(_n, _n);
              double[] e = new double[_n];
              double[] work = new double[_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(double.NaN);
            _v.Set(double.NaN);
            _s.Set(double.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.
            double 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)
//.........这里部分代码省略.........

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

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

              m = new MatrixD(new double[,] {{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(MatrixD.AreNumericallyEqual(m, m * inverse * m));
        }

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

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

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

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

        //--------------------------------------------------------------
        /// <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 MatrixD SolveLinearEquations(MatrixD 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.
              MatrixD 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;
        }