C# SparseMatrix.Multiply方法代码示例

        public double[][] ComputeLaplacianBasic(SparseMatrix L, TriMesh mesh)
        {
            if (L == null)
                throw new Exception("Laplacian matrix is null");

            int n = mesh.Vertices.Count; 
            double[] coordinate = new double[n];
            double[][] lap = new double[3][];  
            coordinate = TriMeshUtil.GetX(mesh).ToArray();
            lap[0]=L.Multiply(coordinate); 
            coordinate = TriMeshUtil.GetY(mesh).ToArray();
            lap[1] = L.Multiply(coordinate); 
            coordinate = TriMeshUtil.GetZ(mesh).ToArray();
            lap[2] = L.Multiply(coordinate);
            return lap;
        }

        /// <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 preconditioner, SparseMatrix matrix, Vector vector, Vector result)
        {
            Assert.AreEqual(typeof(Diagonal), preconditioner.GetType(), "#01");

            // Compute M * result = product
            // compare vector and product. Should be equal
            Vector product = new DenseVector(result.Count);
            matrix.Multiply(result, product);
            for (var i = 0; i < product.Count; i++)
            {
                Assert.IsTrue(vector[i].AlmostEqual(product[i], -Epsilon.Magnitude()), "#02-" + i);
            }
        }

        public void SolvePoissonMatrixAndBackMultiply()
        {
            // Create the matrix
            var matrix = new SparseMatrix(100);
            // Assemble the matrix. We assume we're solving the Poisson equation
            // on a rectangular 10 x 10 grid
            const int GridSize = 10;

            // 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
            Vector<Complex> y = new DenseVector(matrix.RowCount, 1);

            // Create an iteration monitor which will keep track of iterative convergence
            var monitor = new Iterator(new IIterationStopCriterium<Complex>[]
                                       {
                                           new IterationCountStopCriterium(MaximumIterations),
                                           new ResidualStopCriterium(ConvergenceBoundary),
                                           new DivergenceStopCriterium(),
                                           new FailureStopCriterium()
                                       });
            var solver = new MlkBiCgStab(monitor);

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

            // 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 is CalculationConverged, "#04");

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

        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
            Vector<float> y = new DenseVector(matrix.RowCount, 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(new IIterationStopCriterium<float>[]
                                           {
                                               new IterationCountStopCriterium(MaximumIterations),
                                               new ResidualStopCriterium(ConvergenceBoundary),
                                               new DivergenceStopCriterium(),
                                               new FailureStopCriterium()
                                           });
                var solver = new MlkBiCgStab(monitor);

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

                if (!(monitor.Status is CalculationConverged))
                {
                    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.IsTrue(Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#04-" + i);
                }

                return;
            }

            Assert.Fail("Solution was not found in 3 tries");
        }

        private void button6_Click(object sender, EventArgs e)
        {
            DateTime dt = DateTime.Now;

            eps = Math.Pow(10, -Convert.ToInt32(tbeps.Text));

            button7.Enabled = false;
            button8.Enabled = false;
            button9.Enabled = false;

            double[,] mat = tomat(tbA.Text);
            lA = new SparseMatrix(mat);
            n=lA.Rows;
            show("Lib A initiated !");

            lB = new SparseMatrix(n, 1);
            for (int i = 0; i < n; i++)
                for (int j = 0; j < n; j++)
                    lB[i, 0] += mat[i, j] * (j + 1);
            tbB.Clear();
            tbB.AppendText(lB.ToString());

            h = new Householder(lA);
            tbR.Clear();
            tbR.AppendText(h.R().ToString());
            tbQ.Clear();
            tbQ.AppendText(h.Q().ToString());

            Matrix lQp = new SparseMatrix(n, n);
            for (int i = 0; i < n; i++)
                for (int j = 0; j < n; j++)
                    lQp[j, i] = h.Q()[i, j];

            Matrix lBp = lQp.Multiply(lB);
            tbBb.Clear();
            tbBb.AppendText(lBp.ToString());

            show("Lib Q,R calculated !");
            button7.Enabled = true;
            button9.Enabled = true;

            tqr = DateTime.Now - dt;
        }