C# Matrix.CreateMatrix方法代码示例

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

            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
            }

            MatrixR = matrix.Clone();
            MatrixQ = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount);

            for (var i = 0; i < matrix.RowCount; i++)
            {
                MatrixQ.At(i, i, 1.0f);
            }

            var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount);
            var u = new Complex32[minmn][];
            for (var i = 0; i < minmn; i++)
            {
                u[i] = GenerateColumn(MatrixR, i, i);
                ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
            }

            for (var i = minmn - 1; i >= 0; i--)
            {
                ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
            }
        }

        public void TestCreateMatrix()
        {
            var matrix = new Matrix<int>(2, 2);

            var other_matrix = matrix.CreateMatrix(3, 3);
            Assert.IsInstanceOfType(matrix.GetType(), other_matrix);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserQR"/> class. This object will compute the
        /// QR factorization when the constructor is called and cache it's factorization.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        public UserQR(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
            }

            MatrixR = matrix.Clone();
            MatrixQ = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount);

            for (var i = 0; i < matrix.RowCount; i++)
            {
                MatrixQ.At(i, i, 1.0);
            }

            var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount);
            var u = new double[minmn][];
            for (var i = 0; i < minmn; i++)
            {
                u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
                ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
            }

            for (var i = minmn - 1; i >= 0; i--)
            {
                ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
            }
        }

        /// <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<Complex32> matrix)
        {
            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.CreateMatrix(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).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(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="UserGramSchmidt"/> class. This object creates an unitary matrix 
        /// using the modified Gram-Schmidt method.
        /// </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 <paramref name="matrix"/> row count is less then column count</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is rank deficient</exception>
        public UserGramSchmidt(Matrix<Complex32> matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
            }

            MatrixQ = matrix.Clone();
            MatrixR = matrix.CreateMatrix(matrix.ColumnCount, matrix.ColumnCount);

            for (var k = 0; k < MatrixQ.ColumnCount; k++)
            {
                var norm = MatrixQ.Column(k).Norm(2).Real;
                if (norm == 0.0f)
                {
                    throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient);
                }

                MatrixR.At(k, k, norm);
                for (var i = 0; i < MatrixQ.RowCount; i++)
                {
                    MatrixQ.At(i, k, MatrixQ.At(i, k) / norm);
                }

                for (var j = k + 1; j < MatrixQ.ColumnCount; j++)
                {
                    var dot = Complex32.Zero;
                    for (int i = 0; i < MatrixQ.RowCount; i++)
                    {
                        dot += MatrixQ.Column(k)[i].Conjugate() * MatrixQ.Column(j)[i];
                    }
                    
                    MatrixR.At(k, j, dot);
                    for (var i = 0; i < MatrixQ.RowCount; i++)
                    {
                        var value = MatrixQ.At(i, j) - (MatrixQ.At(i, k) * dot);
                        MatrixQ.At(i, j, value);
                    }
                }
            }
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserGramSchmidt"/> class. This object creates an unitary matrix 
        /// using the modified Gram-Schmidt method.
        /// </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 <paramref name="matrix"/> row count is less then column count</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is rank deficient</exception>
        public static UserGramSchmidt Create(Matrix<Complex32> matrix)
        {
            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            var q = matrix.Clone();
            var r = matrix.CreateMatrix(matrix.ColumnCount, matrix.ColumnCount);

            for (var k = 0; k < q.ColumnCount; k++)
            {
                var norm = q.Column(k).L2Norm().Real;
                if (norm == 0.0f)
                {
                    throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient);
                }

                r.At(k, k, norm);
                for (var i = 0; i < q.RowCount; i++)
                {
                    q.At(i, k, q.At(i, k) / norm);
                }

                for (var j = k + 1; j < q.ColumnCount; j++)
                {
                    var dot = Complex32.Zero;
                    for (int i = 0; i < q.RowCount; i++)
                    {
                        dot += q.Column(k)[i].Conjugate() * q.Column(j)[i];
                    }

                    r.At(k, j, dot);
                    for (var i = 0; i < q.RowCount; i++)
                    {
                        var value = q.At(i, j) - (q.At(i, k) * dot);
                        q.At(i, j, value);
                    }
                }
            }

            return new UserGramSchmidt(q, r);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserGramSchmidt"/> class. This object creates an orthogonal matrix 
        /// using the modified Gram-Schmidt method.
        /// </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 <paramref name="matrix"/> row count is less then column count</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is rank deficient</exception>
        public UserGramSchmidt(Matrix<double> matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
            }

            MatrixQ = matrix.Clone();
            MatrixR = matrix.CreateMatrix(matrix.ColumnCount, matrix.ColumnCount);

            for (var k = 0; k < MatrixQ.ColumnCount; k++)
            {
                var norm = MatrixQ.Column(k).Norm(2);
                if (norm == 0.0)
                {
                    throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient);
                }

                MatrixR.At(k, k, norm);
                for (var i = 0; i < MatrixQ.RowCount; i++)
                {
                    MatrixQ.At(i, k, MatrixQ.At(i, k) / norm);
                }

                for (var j = k + 1; j < MatrixQ.ColumnCount; j++)
                {
                    var dot = MatrixQ.Column(k).DotProduct(MatrixQ.Column(j));
                    MatrixR.At(k, j, dot);
                    for (var i = 0; i < MatrixQ.RowCount; i++)
                    {
                        var value = MatrixQ.At(i, j) - (MatrixQ.At(i, k) * dot);
                        MatrixQ.At(i, j, value);
                    }
                }
            }
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserGramSchmidt"/> class. This object creates an orthogonal matrix 
        /// using the modified Gram-Schmidt method.
        /// </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 <paramref name="matrix"/> row count is less then column count</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is rank deficient</exception>
        public static UserGramSchmidt Create(Matrix<float> matrix)
        {
            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            var q = matrix.Clone();
            var r = matrix.CreateMatrix(matrix.ColumnCount, matrix.ColumnCount);

            for (var k = 0; k < q.ColumnCount; k++)
            {
                var norm = q.Column(k).L2Norm();
                if (norm == 0.0)
                {
                    throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient);
                }

                r.At(k, k, norm);
                for (var i = 0; i < q.RowCount; i++)
                {
                    q.At(i, k, q.At(i, k) / norm);
                }

                for (var j = k + 1; j < q.ColumnCount; j++)
                {
                    var dot = q.Column(k).DotProduct(q.Column(j));
                    r.At(k, j, dot);
                    for (var i = 0; i < q.RowCount; i++)
                    {
                        var value = q.At(i, j) - (q.At(i, k) * dot);
                        q.At(i, j, value);
                    }
                }
            }

            return new UserGramSchmidt(q, r);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserQR"/> class. This object will compute the
        /// QR factorization when the constructor is called and cache it's factorization.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="method">The QR factorization method to use.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        public static UserQR Create(Matrix<Complex> matrix, QRMethod method = QRMethod.Full)
        {
            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            Matrix<Complex> q;
            Matrix<Complex> r;

            var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount);
            var u = new Complex[minmn][];

            if (method == QRMethod.Full)
            {
                r = matrix.Clone();
                q = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount);

                for (var i = 0; i < matrix.RowCount; i++)
                {
                    q.At(i, i, 1.0f);
                }

                for (var i = 0; i < minmn; i++)
                {
                    u[i] = GenerateColumn(r, i, i);
                    ComputeQR(u[i], r, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
                }

                for (var i = minmn - 1; i >= 0; i--)
                {
                    ComputeQR(u[i], q, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
                }
            }
            else
            {
                q = matrix.Clone();

                for (var i = 0; i < minmn; i++)
                {
                    u[i] = GenerateColumn(q, i, i);
                    ComputeQR(u[i], q, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
                }

                r = q.SubMatrix(0, matrix.ColumnCount, 0, matrix.ColumnCount);
                q.Clear();

                for (var i = 0; i < matrix.ColumnCount; i++)
                {
                    q.At(i, i, 1.0f);
                }

                for (var i = minmn - 1; i >= 0; i--)
                {
                    ComputeQR(u[i], q, i, matrix.RowCount, i, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
                }
            }

            return new UserQR(q, r, method);
        }

        /// <summary>
        /// Puts the strictly lower triangle of this matrix into the result matrix.
        /// </summary>
        /// <param name="result">Where to store the lower triangle.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <see langword="null" />.</exception>
        /// <exception cref="ArgumentException">If the result matrix's dimensions are not the same as this matrix.</exception>
        public override void StrictlyLowerTriangle(Matrix<Complex> result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
            {
                throw DimensionsDontMatch<ArgumentException>(this, result);
            }

            if (ReferenceEquals(this, result))
            {
                var tmp = result.CreateMatrix(result.RowCount, result.ColumnCount);
                StrictlyLowerTriangle(tmp);
                tmp.CopyTo(result);
            }
            else
            {
                result.Clear();
                StrictlyLowerTriangleImpl(result);
            }
        }

        /// <summary>
        /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
        /// solution matrix and X is the unknown matrix.
        /// </summary>
        /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
        /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
        /// <returns>The result <see cref="Matrix"/>, <c>X</c>.</returns>
        public Matrix Solve(Matrix matrix, Matrix input)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (input == null)
            {
                throw new ArgumentNullException("input");
            }

            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
            Solve(matrix, input, result);
            return result;
        }

        /// <summary>
        /// Multiplies this matrix with transpose of another matrix and returns the result.
        /// </summary>
        /// <param name="other">The matrix to multiply with.</param>
        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>        
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
        /// <returns>The result of multiplication.</returns>
        public override Matrix<float> TransposeAndMultiply(Matrix<float> other)
        {
            var otherDiagonal = other as DiagonalMatrix;
            if (otherDiagonal == null)
            {
                return base.TransposeAndMultiply(other);
            }

            if (ColumnCount != otherDiagonal.ColumnCount)
            {
                throw DimensionsDontMatch<ArgumentException>(this, otherDiagonal);
            }

            var result = other.CreateMatrix(RowCount, other.RowCount);
            TransposeAndMultiply(other, result);
            return result;
        }

        /// <summary>
        /// Multiplies this matrix with another matrix and returns the result.
        /// </summary>
        /// <param name="other">The matrix to multiply with.</param>
        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>        
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
        /// <returns>The result of multiplication.</returns>
        public override Matrix<float> Multiply(Matrix<float> other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (ColumnCount != other.RowCount)
            {
                throw DimensionsDontMatch<ArgumentException>(this, other);
            }

            var result = other.CreateMatrix(RowCount, other.ColumnCount);
            Multiply(other, result);
            return result;
        }