C# DenseVector.ConjugateDotProduct方法代码示例

//.........这里部分代码省略.........

            // Define the temporary vectors
            Vector gtemp = new DenseVector(residuals.Count);

            Vector u = new DenseVector(residuals.Count);
            Vector utemp = new DenseVector(residuals.Count);
            Vector temp = new DenseVector(residuals.Count);
            Vector temp1 = new DenseVector(residuals.Count);
            Vector temp2 = new DenseVector(residuals.Count);

            Vector zd = new DenseVector(residuals.Count);
            Vector zg = new DenseVector(residuals.Count);
            Vector zw = new DenseVector(residuals.Count);

            var d = CreateVectorArray(_startingVectors.Count, residuals.Count);

            // g_0 = r_0
            var g = CreateVectorArray(_startingVectors.Count, residuals.Count);
            residuals.CopyTo(g[k - 1]);

            var w = CreateVectorArray(_startingVectors.Count, residuals.Count);

            // FOR (j = 0, 1, 2 ....)
            var iterationNumber = 0;
            while (ShouldContinue(iterationNumber, xtemp, input, residuals))
            {
                // SOLVE M g~_((j-1)k+k) = g_((j-1)k+k)
                _preconditioner.Approximate(g[k - 1], gtemp);

                // w_((j-1)k+k) = A g~_((j-1)k+k)
                matrix.Multiply(gtemp, w[k - 1]);

                // c_((j-1)k+k) = q^T_1 w_((j-1)k+k)
                c[k - 1] = _startingVectors[0].ConjugateDotProduct(w[k - 1]);
                if (c[k - 1].Real.AlmostEqual(0, 1) && c[k - 1].Imaginary.AlmostEqual(0, 1))
                {
                    throw new Exception("Iterative solver experience a numerical break down");
                }

                // alpha_(jk+1) = q^T_1 r_((j-1)k+k) / c_((j-1)k+k)
                var alpha = _startingVectors[0].ConjugateDotProduct(residuals)/c[k - 1];

                // u_(jk+1) = r_((j-1)k+k) - alpha_(jk+1) w_((j-1)k+k)
                w[k - 1].Multiply(-alpha, temp);
                residuals.Add(temp, u);

                // SOLVE M u~_(jk+1) = u_(jk+1)
                _preconditioner.Approximate(u, temp1);
                temp1.CopyTo(utemp);

                // rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2
                matrix.Multiply(temp1, temp);
                var rho = temp.ConjugateDotProduct(temp);

                // If rho is zero then temp is a zero vector and we're probably
                // about to have zero residuals (i.e. an exact solution).
                // So set rho to 1.0 because in the next step it will turn to zero.
                if (rho.Real.AlmostEqual(0, 1) && rho.Imaginary.AlmostEqual(0, 1))
                {
                    rho = 1.0;
                }

                rho = -u.ConjugateDotProduct(temp)/rho;

                // r_(jk+1) = rho_(j+1) A u~_(jk+1) + u_(jk+1)
                u.CopyTo(residuals);

//.........这里部分代码省略.........

            // w_-1 = 0
            Vector w = new DenseVector(residuals.Count);

            // Define the remaining temporary vectors
            Vector c = new DenseVector(residuals.Count);
            Vector p = new DenseVector(residuals.Count);
            Vector s = new DenseVector(residuals.Count);
            Vector u = new DenseVector(residuals.Count);
            Vector y = new DenseVector(residuals.Count);
            Vector z = new DenseVector(residuals.Count);

            Vector temp = new DenseVector(residuals.Count);
            Vector temp2 = new DenseVector(residuals.Count);
            Vector temp3 = new DenseVector(residuals.Count);

            // for (k = 0, 1, .... )
            var iterationNumber = 0;
            while (ShouldContinue(iterationNumber, xtemp, input, residuals))
            {
                // p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1))
                p.Subtract(u, temp);

                temp.Multiply(beta, temp2);
                residuals.Add(temp2, p);

                // Solve M b_k = p_k
                _preconditioner.Approximate(p, temp);

                // s_k = A b_k
                matrix.Multiply(temp, s);

                // alpha_k = (r*_0 * r_k) / (r*_0 * s_k)
                var alpha = rdash.ConjugateDotProduct(residuals)/rdash.ConjugateDotProduct(s);

                // y_k = t_(k-1) - r_k - alpha_k * w_(k-1) + alpha_k s_k
                s.Subtract(w, temp);
                t.Subtract(residuals, y);

                temp.Multiply(alpha, temp2);
                y.Add(temp2, temp3);
                temp3.CopyTo(y);

                // Store the old value of t in t0
                t.CopyTo(t0);

                // t_k = r_k - alpha_k s_k
                s.Multiply(-alpha, temp2);
                residuals.Add(temp2, t);

                // Solve M d_k = t_k
                _preconditioner.Approximate(t, temp);

                // c_k = A d_k
                matrix.Multiply(temp, c);
                var cdot = c.ConjugateDotProduct(c);

                // cDot can only be zero if c is a zero vector
                // We'll set cDot to 1 if it is zero to prevent NaN's
                // Note that the calculation should continue fine because
                // c.DotProduct(t) will be zero and so will c.DotProduct(y)
                if (cdot.Real.AlmostEqual(0, 1) && cdot.Imaginary.AlmostEqual(0, 1))
                {
                    cdot = 1.0f;
                }

        /// <summary>
        /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
        /// solution vector and x is the unknown vector.
        /// </summary>
        /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
        /// <param name="input">The solution <see cref="Vector"/>, <c>b</c>.</param>
        /// <param name="result">The result <see cref="Vector"/>, <c>x</c>.</param>
        public void Solve(Matrix matrix, Vector input, Vector result)
        {
            // If we were stopped before, we are no longer
            // We're doing this at the start of the method to ensure
            // that we can use these fields immediately.
            _hasBeenStopped = false;

            // Parameters checks
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

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

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

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

            if (result.Count != input.Count)
            {
                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
            }

            if (input.Count != matrix.RowCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
            }

            // Initialize the solver fields
            // Set the convergence monitor
            if (_iterator == null)
            {
                _iterator = Iterator.CreateDefault();
            }

            if (_preconditioner == null)
            {
                _preconditioner = new UnitPreconditioner();
            }

            _preconditioner.Initialize(matrix);

            // Compute r_0 = b - Ax_0 for some initial guess x_0
            // In this case we take x_0 = vector
            // This is basically a SAXPY so it could be made a lot faster
            Vector residuals = new DenseVector(matrix.RowCount);
            CalculateTrueResidual(matrix, residuals, result, input);

            // Choose r~ (for example, r~ = r_0)
            var tempResiduals = residuals.Clone();

            // create seven temporary vectors needed to hold temporary
            // coefficients. All vectors are mangled in each iteration.
            // These are defined here to prevent stressing the garbage collector
            Vector vecP = new DenseVector(residuals.Count);
            Vector vecPdash = new DenseVector(residuals.Count);
            Vector nu = new DenseVector(residuals.Count);
            Vector vecS = new DenseVector(residuals.Count);
            Vector vecSdash = new DenseVector(residuals.Count);
            Vector temp = new DenseVector(residuals.Count);
            Vector temp2 = new DenseVector(residuals.Count);

            // create some temporary float variables that are needed
            // to hold values in between iterations
            Complex32 currentRho = 0;
            Complex32 alpha = 0;
            Complex32 omega = 0;

            var iterationNumber = 0;
            while (ShouldContinue(iterationNumber, result, input, residuals))
            {
                // rho_(i-1) = r~^T r_(i-1) // dotproduct r~ and r_(i-1)
                var oldRho = currentRho;
                currentRho = tempResiduals.ConjugateDotProduct(residuals);

                // if (rho_(i-1) == 0) // METHOD FAILS
                // If rho is only 1 ULP from zero then we fail.
                if (currentRho.Real.AlmostEqual(0, 1) && currentRho.Imaginary.AlmostEqual(0, 1))
                {
                    // Rho-type breakdown
                    throw new Exception("Iterative solver experience a numerical break down");
                }

//.........这里部分代码省略.........