C# Single.SparseMatrix类代码示例

        public void SolveWideMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(2, 3);
            var input = new DenseVector(2);

            var solver = new TFQMR();
            Assert.Throws<ArgumentException>(() => matrix.SolveIterative(input, solver));
        }

        public void SolveLongMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(3, 2);
            var input = new DenseVector(3);

            var solver = new MlkBiCgStab();
            Assert.Throws<ArgumentException>(() => matrix.SolveIterative(input, solver));
        }

        public void SolveWideMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(2, 3);
            var input = new DenseVector(2);

            var solver = new MlkBiCgStab();
            Assert.That(() => matrix.SolveIterative(input, solver), Throws.ArgumentException);
        }

        public void SolveLongMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(3, 2);
            var input = new DenseVector(3);

            var solver = new GpBiCg();
            Assert.That(() => matrix.SolveIterative(input, solver), Throws.ArgumentException);
        }

        /// <summary>
        /// Create unit matrix.
        /// </summary>
        /// <param name="size">Matrix size.</param>
        /// <returns>New unit matrix.</returns>
        internal SparseMatrix CreateUnitMatrix(int size)
        {
            var matrix = new SparseMatrix(size);
            for (var i = 0; i < size; i++)
            {
                matrix[i, i] = 2;
            }

            return matrix;
        }

        /// <summary>
        /// Check the result.
        /// </summary>
        /// <param name="preconditioner">Specific preconditioner.</param>
        /// <param name="matrix">Source matrix.</param>
        /// <param name="vector">Initial vector.</param>
        /// <param name="result">Result vector.</param>
        protected override void CheckResult(IPreconditioner<float> preconditioner, SparseMatrix matrix, Vector<float> vector, Vector<float> result)
        {
            Assert.AreEqual(typeof (UnitPreconditioner<float>), preconditioner.GetType(), "#01");

            // Unit preconditioner is doing nothing. Vector and result should be equal
            for (var i = 0; i < vector.Count; i++)
            {
                Assert.IsTrue(vector[i] == result[i], "#02-" + i);
            }
        }

        /// <summary>
        /// Check the result.
        /// </summary>
        /// <param name="preconditioner">Specific preconditioner.</param>
        /// <param name="matrix">Source matrix.</param>
        /// <param name="vector">Initial vector.</param>
        /// <param name="result">Result vector.</param>
        protected override void CheckResult(IPreconditioner<float> preconditioner, SparseMatrix matrix, Vector<float> vector, Vector<float> result)
        {
            Assert.AreEqual(typeof (DiagonalPreconditioner), preconditioner.GetType(), "#01");

            // Compute M * result = product
            // compare vector and product. Should be equal
            var product = new DenseVector(result.Count);
            matrix.Multiply(result, product);

            for (var i = 0; i < product.Count; i++)
            {
                Assert.IsTrue(((double) vector[i]).AlmostEqualNumbersBetween(product[i], -Epsilon.Magnitude()), "#02-" + i);
            }
        }

        public void CanAddSparseMatricesBothWays()
        {
            var m1 = new SparseMatrix(1, 3);
            var m2 = SparseMatrix.OfArray(new float[,] { { 0, 1, 1 } });
            var sum1 = m1 + m2;
            var sum2 = m2 + m1;
            Assert.IsTrue(sum1.Equals(m2));
            Assert.IsTrue(sum1.Equals(sum2));

            var sparseResult = new SparseMatrix(1, 3);
            sparseResult.Add(m2, sparseResult);
            Assert.IsTrue(sparseResult.Equals(sum1));

            sparseResult = SparseMatrix.OfArray(new float[,] { { 0, 1, 1 } });
            sparseResult.Add(m1, sparseResult);
            Assert.IsTrue(sparseResult.Equals(sum1));

            sparseResult = SparseMatrix.OfArray(new float[,] { { 0, 1, 1 } });
            m1.Add(sparseResult, sparseResult);
            Assert.IsTrue(sparseResult.Equals(sum1));

            sparseResult = SparseMatrix.OfArray(new float[,] { { 0, 1, 1 } });
            sparseResult.Add(sparseResult, sparseResult);
            Assert.IsTrue(sparseResult.Equals(2*sum1));

            var denseResult = new DenseMatrix(1, 3);
            denseResult.Add(m2, denseResult);
            Assert.IsTrue(denseResult.Equals(sum1));

            denseResult = DenseMatrix.OfArray(new float[,] {{0, 1, 1}});
            denseResult.Add(m1, denseResult);
            Assert.IsTrue(denseResult.Equals(sum1));

            var m3 = DenseMatrix.OfArray(new float[,] {{0, 1, 1}});
            var sum3 = m1 + m3;
            var sum4 = m3 + m1;
            Assert.IsTrue(sum3.Equals(m3));
            Assert.IsTrue(sum3.Equals(sum4));
        }

        /// <summary>
        /// Creates a new  <see cref="SparseMatrix"/> and inserts the given column at the given index.
        /// </summary>
        /// <param name="columnIndex">The index of where to insert the column.</param>
        /// <param name="column">The column to insert.</param>
        /// <returns>A new <see cref="SparseMatrix"/> with the inserted column.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="column "/> is <see langword="null" />. </exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> is &lt; zero or &gt; the number of columns.</exception>
        /// <exception cref="ArgumentException">If the size of <paramref name="column"/> != the number of rows.</exception>
        public override Matrix<float> InsertColumn(int columnIndex, Vector<float> column)
        {
            if (column == null)
            {
                throw new ArgumentNullException("column");
            }

            if (columnIndex < 0 || columnIndex > ColumnCount)
            {
                throw new ArgumentOutOfRangeException("columnIndex");
            }

            if (column.Count != RowCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension, "column");
            }

            var result = new SparseMatrix(RowCount, ColumnCount + 1);

            for (var i = 0; i < columnIndex; i++)
            {
                result.SetColumn(i, Column(i));
            }

            result.SetColumn(columnIndex, column);

            for (var i = columnIndex + 1; i < ColumnCount + 1; i++)
            {
                result.SetColumn(i, Column(i - 1));
            }

            return result;
        }

        /// <summary>
        /// Initializes a square <see cref="SparseMatrix"/> with all zero's except for ones on the diagonal.
        /// </summary>
        /// <param name="order">the size of the square matrix.</param>
        /// <returns>Identity <c>SparseMatrix</c></returns>
        /// <exception cref="ArgumentException">
        /// If <paramref name="order"/> is less than one.
        /// </exception>
        public static SparseMatrix Identity(int order)
        {
            var m = new SparseMatrix(order);
            var mStorage = m.Raw;

            mStorage.ValueCount = order;
            mStorage.Values = new float[order];
            mStorage.ColumnIndices = new int[order];

            for (var i = 0; i < order; i++)
            {
                mStorage.Values[i] = 1f;
                mStorage.ColumnIndices[i] = i;
                mStorage.RowPointers[i] = i;
            }

            return m;
        }

        public void SolvePoissonMatrixAndBackMultiply()
        {
            // Create the matrix
            var matrix = new SparseMatrix(25);

            // Assemble the matrix. We assume we're solving the Poisson equation
            // on a rectangular 5 x 5 grid
            const int GridSize = 5;

            // The pattern is:
            // 0 .... 0 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 0 0 ... 0
            for (var i = 0; i < matrix.RowCount; i++)
            {
                // Insert the first set of -1's
                if (i > (GridSize - 1))
                {
                    matrix[i, i - GridSize] = -1;
                }

                // Insert the second set of -1's
                if (i > 0)
                {
                    matrix[i, i - 1] = -1;
                }

                // Insert the centerline values
                matrix[i, i] = 4;

                // Insert the first trailing set of -1's
                if (i < matrix.RowCount - 1)
                {
                    matrix[i, i + 1] = -1;
                }

                // Insert the second trailing set of -1's
                if (i < matrix.RowCount - GridSize)
                {
                    matrix[i, i + GridSize] = -1;
                }
            }

            // Create the y vector
            var y = DenseVector.Create(matrix.RowCount, i => 1);

            // Due to datatype "float" it can happen that solution will not converge for specific random starting vectors
            // That's why we will do 3 tries
            for (var iteration = 0; iteration <= 3; iteration++)
            {
                // Create an iteration monitor which will keep track of iterative convergence
                var monitor = new Iterator<float>(
                    new IterationCountStopCriterium<float>(MaximumIterations),
                    new ResidualStopCriterium<float>(ConvergenceBoundary),
                    new DivergenceStopCriterium<float>(),
                    new FailureStopCriterium<float>());

                var solver = new MlkBiCgStab();

                // Solve equation Ax = y
                Vector<float> x;
                try
                {
                    x = matrix.SolveIterative(y, solver, monitor);
                }
                catch (Exception)
                {
                    continue;
                }

                if (monitor.Status != IterationStatus.Converged)
                {
                    continue;
                }

                // Now compare the results
                Assert.IsNotNull(x, "#02");
                Assert.AreEqual(y.Count, x.Count, "#03");

                // Back multiply the vector
                var z = matrix.Multiply(x);

                // Now compare the vectors
                for (var i = 0; i < y.Count; i++)
                {
                    Assert.GreaterOrEqual(ConvergenceBoundary, Math.Abs(y[i] - z[i]), "#05-" + i);
                }

                return;
            }
        }

 /// <summary>
 /// Check the result.
 /// </summary>
 /// <param name="preconditioner">Specific preconditioner.</param>
 /// <param name="matrix">Source matrix.</param>
 /// <param name="vector">Initial vector.</param>
 /// <param name="result">Result vector.</param>
 protected abstract void CheckResult(IPreconditioner<float> preconditioner, SparseMatrix matrix, Vector<float> vector, Vector<float> result);

 /// <summary>
 /// Create a new sparse matrix with the diagonal as a copy of the given vector.
 /// This new matrix will be independent from the vector.
 /// A new memory block will be allocated for storing the matrix.
 /// </summary>
 public static SparseMatrix OfDiagonalVector(int rows, int columns, Vector<float> diagonal)
 {
     var m = new SparseMatrix(rows, columns);
     m.SetDiagonal(diagonal);
     return m;
 }

 /// <summary>
 /// Create a new sparse matrix with the diagonal as a copy of the given array.
 /// This new matrix will be independent from the array.
 /// A new memory block will be allocated for storing the matrix.
 /// </summary>
 public static SparseMatrix OfDiagonalArray(int rows, int columns, float[] diagonal)
 {
     var m = new SparseMatrix(rows, columns);
     m.SetDiagonal(diagonal);
     return m;
 }

        public void SolvePoissonMatrixAndBackMultiply()
        {
            // Create the matrix
            var matrix = new SparseMatrix(25);

            // Assemble the matrix. We assume we're solving the Poisson equation
            // on a rectangular 5 x 5 grid
            const int GridSize = 5;

            // The pattern is:
            // 0 .... 0 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 0 0 ... 0
            for (var i = 0; i < matrix.RowCount; i++)
            {
                // Insert the first set of -1's
                if (i > (GridSize - 1))
                {
                    matrix[i, i - GridSize] = -1;
                }

                // Insert the second set of -1's
                if (i > 0)
                {
                    matrix[i, i - 1] = -1;
                }

                // Insert the centerline values
                matrix[i, i] = 4;

                // Insert the first trailing set of -1's
                if (i < matrix.RowCount - 1)
                {
                    matrix[i, i + 1] = -1;
                }

                // Insert the second trailing set of -1's
                if (i < matrix.RowCount - GridSize)
                {
                    matrix[i, i + GridSize] = -1;
                }
            }

            // Create the y vector
            var y = DenseVector.Create(matrix.RowCount, i => 1);

            // Create an iteration monitor which will keep track of iterative convergence
            var monitor = new Iterator<float>(
                new IterationCountStopCriterium<float>(MaximumIterations),
                new ResidualStopCriterium(ConvergenceBoundary),
                new DivergenceStopCriterium(),
                new FailureStopCriterium());

            var solver = new TFQMR();

            // Solve equation Ax = y
            var x = matrix.SolveIterative(y, solver, monitor);

            // Now compare the results
            Assert.IsNotNull(x, "#02");
            Assert.AreEqual(y.Count, x.Count, "#03");

            // Back multiply the vector
            var z = matrix.Multiply(x);

            // Check that the solution converged
            Assert.IsTrue(monitor.Status == IterationStatus.Converged, "#04");

            // Now compare the vectors
            for (var i = 0; i < y.Count; i++)
            {
                Assert.IsTrue(Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#05-" + i);
            }
        }