C# Matrix.ArrayCopy方法代码示例

        /* ------------------------
           Constructor
         * ------------------------ */
        ///<summary>QR Decomposition, computed by Householder reflections.</summary>
        ///<param name="A">Rectangular matrix</param>
        ///<returns>Structure to access R and the Householder vectors and compute Q.</returns>
        public QRDecomposition(Matrix A)
        {
            // Initialize.
              QR = A.ArrayCopy();
              m = A.RowDimension;
              n = A.ColumnDimension;
              Rdiag = new double[n];

              // Main loop.
              for (int k = 0; k < n; k++) {
             // Compute 2-norm of k-th column without under/overflow.
             double nrm = 0;
             for (int i = k; i < m; i++) {
            nrm = Maths.Hypot(nrm,QR[i][k]);
             }

             if (nrm != 0.0) {
            // Form k-th Householder vector.
            if (QR[k][k] < 0) {
               nrm = -nrm;
            }
            for (int i = k; i < m; i++) {
               QR[i][k] /= nrm;
            }
            QR[k][k] += 1.0;

            // Apply transformation to remaining columns.
            for (int j = k+1; j < n; j++) {
               double s = 0.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];
               }
            }
             }
             Rdiag[k] = -nrm;
              }
        }

        /* ------------------------
           Constructor
         * ------------------------ */
        /** LU Decomposition
           @param  A   Rectangular matrix
           @return     Structure to access L, U and piv.
           */
        public LUDecomposition(Matrix A)
        {
            // Use a "left-looking", dot-product, Crout/Doolittle algorithm.

              LU = A.ArrayCopy();
              m = A.RowDimension;
              n = A.ColumnDimension;
              piv = new int[m];
              for (int i = 0; i < m; i++) {
             piv[i] = i;
              }
              pivsign = 1;
              double[] LUrowi;
              double[] LUcolj = new double[m];

              // Outer loop.

              for (int j = 0; j < n; j++) {

             // Make a copy of the j-th column to localize references.

             for (int i = 0; i < m; i++) {
            LUcolj[i] = LU[i][j];
             }

             // Apply previous transformations.

             for (int i = 0; i < m; i++) {
            LUrowi = LU[i];

            // Most of the time is spent in the following dot product.

            int kmax = Math.Min(i,j);
            double s = 0.0;
            for (int k = 0; k < kmax; k++) {
               s += LUrowi[k]*LUcolj[k];
            }

            LUrowi[j] = LUcolj[i] -= s;
             }

             // Find pivot and exchange if necessary.

             int p = j;
             for (int i = j+1; i < m; i++) {
            if (Math.Abs(LUcolj[i]) > Math.Abs(LUcolj[p])) {
               p = i;
            }
             }
             if (p != j) {
            for (int k = 0; k < n; k++) {
               double t = LU[p][k]; LU[p][k] = LU[j][k]; LU[j][k] = t;
            }
            int k2 = piv[p]; piv[p] = piv[j]; piv[j] = k2;
            pivsign = -pivsign;
             }

             // Compute multipliers.

             if (j < m && LU[j][j] != 0.0) {
            for (int i = j+1; i < m; i++) {
               LU[i][j] /= LU[j][j];
            }
             }
              }
        }

        /* ------------------------
           Constructor
         * ------------------------ */
        ///<summary>Construct the singular value decomposition</summary>
        ///<param name="Arg">Rectangular matrix</param>
        ///<returns>Structure to access U, S and V.</returns>
        public SingularValueDecomposition(Matrix Arg)
        {
            // Derived from LINPACK code.
              // Initialize.
              double[][] A = Arg.ArrayCopy();
              m = Arg.RowDimension;
              n = Arg.ColumnDimension;

              /* Apparently the failing cases are only a proper subset of (m<n),
             so let's not throw error.  Correct fix to come later?
              if (m<n) {
              throw new System.ArgumentException("Jama SVD only works for m >= n"); }
              */
              int nu = Math.Min(m,n);
              s = new double [Math.Min(m+1,n)];
              U = new double [m][];
              for (int i = 0; i < m; i++)//Added by KJ
              {
              U[i] = new double[m];
              }
              V = new double [n][];
              for (int i = 0; i < n; i++)//Added by KJ
              {
              V[i] = new double[n];
              }
              double[] e = new double [n];
              double[] work = new double [m];
              bool wantu = true;
              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] = Maths.Hypot(s[k],A[i][k]);
            }
            if (s[k] != 0.0) {
               if (A[k][k] < 0.0) {
                  s[k] = -s[k];
               }
               for (int i = k; i < m; i++) {
                  A[i][k] /= s[k];
               }
               A[k][k] += 1.0;
            }
            s[k] = -s[k];
             }
             for (int j = k+1; j < n; j++) {
            if ((k < nct) & (s[k] != 0.0))  {

            // Apply the transformation.

               double t = 0;
               for (int i = k; i < m; i++) {
                  t += A[i][k]*A[i][j];
               }
               t = -t/A[k][k];
               for (int i = k; i < m; i++) {
                  A[i][j] += t*A[i][k];
               }
            }

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

            e[j] = A[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] = A[i][k];
            }
             }
             if (k < nrt) {

            // Compute the k-th row transformation and place the
            // k-th super-diagonal in e[k].
            // Compute 2-norm without under/overflow.
            e[k] = 0;
            for (int i = k+1; i < n; i++) {
               e[k] = Maths.Hypot(e[k],e[i]);
            }
//.........这里部分代码省略.........

        ///<Sumary>Least squares solution of A*X = B</Sumary>
        ///<param name="B">A Matrix with as many rows as A and any number of columns.</param>
        ///<returns> X that minimizes the two norm of Q*R*X-B.</returns>
        public Matrix Solve(Matrix B)
        {
            if (B.RowDimension != m) {
             throw new ArgumentException("Matrix row dimensions must agree.");
              }
              if (!IsFullRank()) {
             throw new Exception("Matrix is rank deficient.");
              }

              // Copy right hand side
              int nx = B.ColumnDimension;
              double[][] X = B.ArrayCopy();

              // Compute Y = transpose(Q)*B
              for (int k = 0; k < n; k++) {
             for (int j = 0; j < nx; j++) {
            double s = 0.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 < nx; j++) {
            X[k][j] /= Rdiag[k];
             }
             for (int i = 0; i < k; i++) {
            for (int j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*QR[i][k];
            }
             }
              }
              return (new Matrix(X,n,nx).GetMatrix(0,n-1,0,nx-1));
        }

        ///<summary>Solve A*X = B</summary>
        ///<param name="B">A Matrix with as many rows as A and any number of columns.</param>
        ///<returns>X so that L*L'*X = B</returns>
        ///<exception cref="ArgumentException">System.ArgumentException  Matrix row dimensions must agree.</exception>
        ///<exception cref="Exception">System.Exception  Matrix is not symmetric positive definite.</exception>
        public Matrix Solve(Matrix B)
        {
            if (B.RowDimension != n) {
             throw new ArgumentException("Matrix row dimensions must agree.");
              }
              if (!isspd) {
             throw new Exception("Matrix is not symmetric positive definite.");
              }

              // Copy right hand side.
              double[][] X = B.ArrayCopy();
              int nx = B.ColumnDimension;

              // Solve L*Y = B;
              for (int k = 0; k < n; k++) {
            for (int j = 0; j < nx; j++) {
               for (int i = 0; i < k ; i++) {
                   X[k][j] -= X[i][j]*L[k][i];
               }
               X[k][j] /= L[k][k];
            }
              }

              // Solve L'*X = Y;
              for (int k = n-1; k >= 0; k--) {
            for (int j = 0; j < nx; j++) {
               for (int i = k+1; i < n ; i++) {
                   X[k][j] -= X[i][j]*L[i][k];
               }
               X[k][j] /= L[k][k];
            }
              }

              return new Matrix(X,n,nx);
        }