C# Matrix.ToArray方法代码示例

        public static void InverseDeterminant(Matrix matrix, out Matrix inverse, out double determinant)
        {
            int n = matrix.Rows;

            if (matrix.Columns != n)
            {
                throw new ArgumentException("The matrix isn't a square matrix.");
            }

            double[,] a = matrix.ToArray();

            if (!trfac.spdmatrixcholesky(ref a, n, false))
            {
                throw new ArithmeticException();
            }

            determinant = matdet.spdmatrixcholeskydet(ref a, n);

            int info = 0;
            matinv.matinvreport rep = new matinv.matinvreport();
            matinv.spdmatrixcholeskyinverse(ref a, n, false, ref info, ref rep);

            for (int i = 0; i < n; i++)
            {
                for (int j = i + 1; j < n; j++)
                {
                    a[i, j] = a[j, i];
                }
            }

            inverse = new Matrix(a);
        }

        private LUDecomposition(Matrix matrix)
        {
            m = matrix.Rows;
            n = matrix.Columns;
            q = Math.Min(m, n);

            a = matrix.ToArray();
            pivots = new int[0];

            if (n == 0 && m == 0)
            {
                // Ignore degenerate case.
                return;
            }

            trfac.rmatrixlu(ref a, m, n, ref pivots);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public static UserEvd Create(Matrix<float> matrix)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            var eigenVectors = Matrix<float>.Build.SameAs(matrix, order, order);
            var blockDiagonal = Matrix<float>.Build.SameAs(matrix, order, order);
            var eigenValues = new LinearAlgebra.Complex.DenseVector(order);

            var isSymmetric = true;

            for (var i = 0; isSymmetric && i < order; i++)
            {
                for (var j = 0; isSymmetric && j < order; j++)
                {
                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i);
                }
            }

            var d = new float[order];
            var e = new float[order];

            if (isSymmetric)
            {
                matrix.CopyTo(eigenVectors);
                d = eigenVectors.Row(order - 1).ToArray();

                SymmetricTridiagonalize(eigenVectors, d, e, order);
                SymmetricDiagonalize(eigenVectors, d, e, order);
            }
            else
            {
                var matrixH = matrix.ToArray();

                NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(eigenVectors, matrixH, d, e, order);
            }

            for (var i = 0; i < order; i++)
            {
                blockDiagonal.At(i, i, d[i]);

                if (e[i] > 0)
                {
                    blockDiagonal.At(i, i + 1, e[i]);
                }
                else if (e[i] < 0)
                {
                    blockDiagonal.At(i, i - 1, e[i]);
                }
            }

            for (var i = 0; i < order; i++)
            {
                eigenValues[i] = new Complex(d[i], e[i]);
            }

            return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public UserEvd(Matrix<Complex32> matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            MatrixEv = DenseMatrix.Identity(order);
            MatrixD = matrix.CreateMatrix(order, order);
            VectorEv = new LinearAlgebra.Complex.DenseVector(order);
           
            IsSymmetric = true;

            for (var i = 0; IsSymmetric && i < order; i++)
            {
                for (var j = 0; IsSymmetric && j < order; j++)
                {
                    IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
                }
            }

            if (IsSymmetric)
            {
                var matrixCopy = matrix.ToArray();
                var tau = new Complex32[order];
                var d = new float[order];
                var e = new float[order];

                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
                SymmetricDiagonalize(d, e, order);
                SymmetricUntridiagonalize(matrixCopy, tau, order);

                for (var i = 0; i < order; i++)
                {
                    VectorEv[i] = new Complex(d[i], e[i]);
                }
            }
            else
            {
                var matrixH = matrix.ToArray();
                NonsymmetricReduceToHessenberg(matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(matrixH, order);
            }

            for (var i = 0; i < VectorEv.Count; i++)
            {
                MatrixD.At(i, i, (Complex32)VectorEv[i]);
            }
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public static UserEvd Create(Matrix<Complex> matrix)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            var eigenVectors = DenseMatrix.Identity(order);
            var blockDiagonal = matrix.CreateMatrix(order, order);
            var eigenValues = new DenseVector(order);

            var isSymmetric = true;

            for (var i = 0; isSymmetric && i < order; i++)
            {
                for (var j = 0; isSymmetric && j < order; j++)
                {
                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
                }
            }

            if (isSymmetric)
            {
                var matrixCopy = matrix.ToArray();
                var tau = new Complex[order];
                var d = new double[order];
                var e = new double[order];

                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
                SymmetricDiagonalize(eigenVectors, d, e, order);
                SymmetricUntridiagonalize(eigenVectors, matrixCopy, tau, order);

                for (var i = 0; i < order; i++)
                {
                    eigenValues[i] = new Complex(d[i], e[i]);
                }
            }
            else
            {
                var matrixH = matrix.ToArray();
                NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(eigenVectors, eigenValues, matrixH, order);
            }

            blockDiagonal.SetDiagonal(eigenValues);

            return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public static UserEvd Create(Matrix<Complex32> matrix, Symmetricity symmetricity)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            var eigenVectors = DenseMatrix.CreateIdentity(order);
            var blockDiagonal = Matrix<Complex32>.Build.SameAs(matrix, order, order);
            var eigenValues = new LinearAlgebra.Complex.DenseVector(order);

            bool isSymmetric;
            switch (symmetricity)
            {
                case Symmetricity.Hermitian:
                    isSymmetric = true;
                    break;
                case Symmetricity.Asymmetric:
                    isSymmetric = false;
                    break;
                default:
                    isSymmetric = matrix.IsHermitian();
                    break;
            }

            if (isSymmetric)
            {
                var matrixCopy = matrix.ToArray();
                var tau = new Complex32[order];
                var d = new float[order];
                var e = new float[order];

                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
                SymmetricDiagonalize(eigenVectors, d, e, order);
                SymmetricUntridiagonalize(eigenVectors, matrixCopy, tau, order);

                for (var i = 0; i < order; i++)
                {
                    eigenValues[i] = new Complex(d[i], e[i]);
                }
            }
            else
            {
                var matrixH = matrix.ToArray();
                NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(eigenVectors, eigenValues, matrixH, order);
            }

            for (var i = 0; i < eigenValues.Count; i++)
            {
                blockDiagonal.At(i, i, (Complex32) eigenValues[i]);
            }

            return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public UserEvd(Matrix<Complex> matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            MatrixEv = DenseMatrix.Identity(order);
            MatrixD = matrix.CreateMatrix(order, order);
            VectorEv = new DenseVector(order);

            IsSymmetric = true;

            for (var i = 0; i < order & IsSymmetric; i++)
            {
                for (var j = 0; j < order & IsSymmetric; j++)
                {
                    IsSymmetric &= matrix[i, j] == matrix[j, i].Conjugate();
                }
            }

            if (IsSymmetric)
            {
                var matrixCopy = matrix.ToArray();
                var tau = new Complex[order];
                var d = new double[order];
                var e = new double[order];

                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
                SymmetricDiagonalize(d, e, order);
                SymmetricUntridiagonalize(matrixCopy, tau, order);

                for (var i = 0; i < order; i++)
                {
                    VectorEv[i] = new Complex(d[i], e[i]);
                }
            }
            else
            {
                var matrixH = matrix.ToArray();
                NonsymmetricReduceToHessenberg(matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(matrixH, order);
            }

            MatrixD.SetDiagonal(VectorEv);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public UserEvd(Matrix<float> matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            MatrixEv = matrix.CreateMatrix(order, order);
            MatrixD = matrix.CreateMatrix(order, order);
            VectorEv = new LinearAlgebra.Complex.DenseVector(order);
           
            IsSymmetric = true;

            for (var i = 0; IsSymmetric && i < order; i++)
            {
                for (var j = 0; IsSymmetric && j < order; j++)
                {
                    IsSymmetric &= matrix.At(i, j) == matrix.At(j, i);
                }
            }

            var d = new float[order];
            var e = new float[order];

            if (IsSymmetric)
            {
                matrix.CopyTo(MatrixEv);
                d = MatrixEv.Row(order - 1).ToArray();

                SymmetricTridiagonalize(d, e, order);
                SymmetricDiagonalize(d, e, order);
            }
            else
            {
                var matrixH = matrix.ToArray();

                NonsymmetricReduceToHessenberg(matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(matrixH, d, e, order);
            }

            for (var i = 0; i < order; i++)
            {
                MatrixD.At(i, i, d[i]);

                if (e[i] > 0)
                {
                    MatrixD.At(i, i + 1, e[i]);
                }
                else if (e[i] < 0)
                {
                    MatrixD.At(i, i - 1, e[i]);
                }
            }

            for (var i = 0; i < order; i++)
            {
                VectorEv[i] = new Complex(d[i], e[i]);
            }
        }

        public void PrintMatrix(Matrix matrix, float[,] matrixAsArray, string message)
        {
            float[,] matrixArr;
            if (matrixAsArray != null)
            {
                matrixArr = (float[,])matrixAsArray.Clone();
            }
            else
            {
                matrixArr = matrix.ToArray();
            }

            string path = MyConstants.IMAGE_PROCESSING_TEXT_FILES + "KlTransofmData" + Now.Millisecond.ToString() + ".txt";
            using (StreamWriter sw = System.IO.File.Exists(Server.MapPath(path)) ? System.IO.File.AppendText(Server.MapPath(path)) : System.IO.File.CreateText(Server.MapPath(path)))
            {
                sw.WriteLine("-----" + message + "-------");
                for (int i = 0; i < matrixArr.GetLength(0); i++)
                {
                    for (int j = 0; j < matrixArr.GetLength(1); j++)
                    {
                        // DeleteFileIfExist(path);
                        sw.Write(" " + matrixArr[i, j]);
                    }
                    sw.WriteLine();
                }
                sw.WriteLine("-----" + message + " end ");
            }



        }

 public void          Parameter                     (string name, Matrix      matrix)
 {
     int handle = Gl.glGetHandleARB(Gl.GL_PROGRAM_OBJECT_ARB);
     int var    = Gl.glGetUniformLocation(this.programHandle, name);
     Gl.glUniformMatrix4fv(var, 1, 1, matrix.ToArray());
 }

        /// <summary>
        /// Updates the underlying interpolation mode, if needed.
        /// </summary>
        private void UpdateModel()
        {
            switch (this.interpolationType)
            {
                // Right now only the Least Squares requires this operation to be done.
                case EInterpolationType.LEAST_SQUARES:
                    this.quadraticModel = new Fairmat.Optimization.QuadraticModel();

                    // Unroll matrix and coordinate vectors in order to make it suitable
                    // for the Quadratic model implementation.
                    int n = this.values.R * this.values.C;
                    Matrix xy = new Matrix(n, 2);
                    Vector z = new Vector(n);
                    int count = 0;
                    for (int x = 0; x < this.coordinatesX.Length; x++)
                    {
                        for (int y = 0; y < this.coordinatesY.Length; y++)
                        {
                            xy[count, Range.All] = ((Matrix)new Vector() { this[x, -1], this[-1, y] }).T;
                            z[count] = this[x, y];
                            count++;
                        }
                    }

                    this.quadraticModel.Estimate(xy.ToArray() as double[,], z.ToArray() as double[]);
                    break;
                default:
                    break;
            }
        }

        /// <summary>Check for symmetry, then construct the eigenvalue decomposition</summary>
        public EigenvalueDecomposition(Matrix matrix)
        {
            Contract.Assert(matrix != null);

            if (!matrix.IsSquare)
            {
                throw new InvalidOperationException();
            }

            double[][] A = matrix.ToArray();
            n = matrix.Width;
            V = new double[n][];
            for (int i = 0; i < n; i++)
            {
                V[i] = new double[n];
            }
            d = new double[n];
            e = new double[n];

            issymmetric = true;
            for (int j = 0; (j < n) & issymmetric; j++)
            {
                for (int i = 0; (i < n) & issymmetric; i++)
                {
                    issymmetric = (A[i][j] == A[j][i]);
                }
            }

            if (issymmetric)
            {
                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        V[i][j] = A[i][j];
                    }
                }

                // Tridiagonalize.
                tred2();

                // Diagonalize.
                tql2();
            }
            else
            {
                H = new double[n][];
                for (int i2 = 0; i2 < n; i2++)
                {
                    H[i2] = new double[n];
                }
                ort = new double[n];

                for (int j = 0; j < n; j++)
                {
                    for (int i = 0; i < n; i++)
                    {
                        H[i][j] = A[i][j];
                    }
                }

                // Reduce to Hessenberg form.
                orthes();

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

        public static void TestInverse()
        {
            //----------------------
            //| 0.18 | 0.41 | 0.14 |
            //| 0.60 | 0.24 | 0.30 |
            //| 0.57 | 0.99 | 0.97 |
            //----------------------

            Matrix matrix = new Matrix(3, 3);
            Matrix matrix2 = new Matrix(4, 4);
            matrix.SetValue(0, 0, 0.18);
            matrix.SetValue(0, 1, 0.41);
            matrix.SetValue(0, 2, 0.14);
            matrix.SetValue(1, 0, 0.60);
            matrix.SetValue(1, 1, 0.24);
            matrix.SetValue(1, 2, 0.30);
            matrix.SetValue(2, 0, 0.57);
            matrix.SetValue(2, 1, 0.99);
            matrix.SetValue(2, 2, 0.97);
            double[,] test = matrix.ToArray();
            for (uint i = 0; i < matrix.Columns; i++)
            {
                for (uint j = 0; j < matrix.Rows; j++)
                {
                    matrix2.SetValue(i + 1, j + 1, matrix.GetValue(i, j));
                }
            }

            //LU分解による方法
            Matrix inv = new Matrix(3, 3);
            int sig;
            Permutation perm = new Permutation(3);
            perm.Initialize();
            LinearAlgebra.LUDecomposition(ref matrix, ref perm, out sig);
            LinearAlgebra.LUInvert(matrix, perm, ref inv);
            for (uint i = 0; i < inv.Columns; i++)
            {
                for (uint j = 0; j < inv.Rows; j++)
                {
                    Console.Write(inv.GetValue(i, j).ToString("F4").PadLeft(8) + " | ");
                }
                Console.WriteLine();
            }
            Console.WriteLine();

            //部分行列のテスト
            perm.Initialize();
            Matrix inv2 = new Matrix(4, 4);
            MatrixView mView = new MatrixView(matrix2, 1, 1, 3, 3);
            MatrixView mViewINV = new MatrixView(inv2, 0, 1, 3, 3);
            LinearAlgebra.LUDecomposition(ref mView, ref perm, out sig);
            LinearAlgebra.LUInvert(mView, perm, ref mViewINV);
            for (uint i = 0; i < mViewINV.ColumnSize; i++)
            {
                for (uint j = 0; j < mViewINV.RowSize; j++)
                {
                    Console.Write(mViewINV.GetValue(i, j).ToString("F4").PadLeft(8) + " | ");
                }
                Console.WriteLine();
            }
            Console.WriteLine();

            for (uint i = 0; i < inv2.Columns; i++)
            {
                for (uint j = 0; j < inv2.Rows; j++)
                {
                    Console.Write(inv2.GetValue(i, j).ToString("F4").PadLeft(8) + " | ");
                }
                Console.WriteLine();
            }

            Console.Read();
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public static UserEvd Create(Matrix<double> matrix, Symmetricity symmetricity)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            var eigenVectors = Matrix<double>.Build.SameAs(matrix, order, order);
            var blockDiagonal = Matrix<double>.Build.SameAs(matrix, order, order);
            var eigenValues = new LinearAlgebra.Complex.DenseVector(order);

            bool isSymmetric;
            switch (symmetricity)
            {
                case Symmetricity.Symmetric:
                case Symmetricity.Hermitian:
                    isSymmetric = true;
                    break;
                case Symmetricity.Asymmetric:
                    isSymmetric = false;
                    break;
                default:
                    isSymmetric = matrix.IsSymmetric();
                    break;
            }

            var d = new double[order];
            var e = new double[order];

            if (isSymmetric)
            {
                matrix.CopyTo(eigenVectors);
                d = eigenVectors.Row(order - 1).ToArray();

                SymmetricTridiagonalize(eigenVectors, d, e, order);
                SymmetricDiagonalize(eigenVectors, d, e, order);
            }
            else
            {
                var matrixH = matrix.ToArray();

                NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(eigenVectors, matrixH, d, e, order);
            }

            for (var i = 0; i < order; i++)
            {
                blockDiagonal.At(i, i, d[i]);

                if (e[i] > 0)
                {
                    blockDiagonal.At(i, i + 1, e[i]);
                }
                else if (e[i] < 0)
                {
                    blockDiagonal.At(i, i - 1, e[i]);
                }
            }

            for (var i = 0; i < order; i++)
            {
                eigenValues[i] = new Complex(d[i], e[i]);
            }

            return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric);
        }

 public void ToArrayTest()
 {
     double[,] data = null; // TODO: инициализация подходящего значения
     Matrix target = new Matrix(data); // TODO: инициализация подходящего значения
     double[,] expected = null; // TODO: инициализация подходящего значения
     double[,] actual;
     actual = target.ToArray();
     Assert.AreEqual(expected, actual);
     Assert.Inconclusive("Проверьте правильность этого метода теста.");
 }