C# DenseVector.Add方法代码示例

//.........这里部分代码省略.........
            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 double variables that are needed
            // to hold values in between iterations
            Complex currentRho = 0;
            Complex alpha = 0;
            Complex 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.DotProduct(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");
                }

                if (iterationNumber != 0)
                {
                    // beta_(i-1) = (rho_(i-1)/rho_(i-2))(alpha_(i-1)/omega(i-1))
                    var beta = (currentRho / oldRho) * (alpha / omega);

                    // p_i = r_(i-1) + beta_(i-1)(p_(i-1) - omega_(i-1) * nu_(i-1))
                    nu.Multiply(-omega, temp);
                    vecP.Add(temp, temp2);
                    temp2.CopyTo(vecP);

                    vecP.Multiply(beta, vecP);
                    vecP.Add(residuals, temp2);
                    temp2.CopyTo(vecP);
                }
                else
                {
                    // p_i = r_(i-1)
                    residuals.CopyTo(vecP);
                }

                // SOLVE Mp~ = p_i // M = preconditioner
                _preconditioner.Approximate(vecP, vecPdash);
                
                // nu_i = Ap~
                matrix.Multiply(vecPdash, nu);

                // alpha_i = rho_(i-1)/ (r~^T nu_i) = rho / dotproduct(r~ and nu_i)
                alpha = currentRho * 1 / tempResiduals.DotProduct(nu);

                // s = r_(i-1) - alpha_i nu_i
                nu.Multiply(-alpha, temp);
                residuals.Add(temp, vecS);

                // Check if we're converged. If so then stop. Otherwise continue;
                // Calculate the temporary result. 
                // Be careful not to change any of the temp vectors, except for
                // temp. Others will be used in the calculation later on.
                // x_i = x_(i-1) + alpha_i * p^_i + s^_i
                vecPdash.Multiply(alpha, temp);
                temp.Add(vecSdash, temp2);

//.........这里部分代码省略.........
            var rho = tau * tau;

            // Calculate the initial values for v
            // M temp = yEven
            _preconditioner.Approximate(yeven, temp);

            // v = A temp
            matrix.Multiply(temp, v);

            // Set uOdd
            v.CopyTo(ueven);

            // Start the iteration
            var iterationNumber = 0;
            while (ShouldContinue(iterationNumber, result, input, pseudoResiduals))
            {
                // First part of the step, the even bit
                if (IsEven(iterationNumber))
                {
                    // sigma = (v, r)
                    var sigma = v.DotProduct(r);
                    if (sigma.AlmostEqual(0, 1))
                    {
                        // FAIL HERE
                        _iterator.IterationCancelled();
                        break;
                    }

                    // alpha = rho / sigma
                    alpha = rho / sigma;

                    // yOdd = yEven - alpha * v
                    v.Multiply(-alpha, temp1);
                    yeven.Add(temp1, yodd);

                    // Solve M temp = yOdd
                    _preconditioner.Approximate(yodd, temp);

                    // uOdd = A temp
                    matrix.Multiply(temp, uodd);
                }

                // The intermediate step which is equal for both even and
                // odd iteration steps.
                // Select the correct vector
                var uinternal = IsEven(iterationNumber) ? ueven : uodd;
                var yinternal = IsEven(iterationNumber) ? yeven : yodd;

                // pseudoResiduals = pseudoResiduals - alpha * uOdd
                uinternal.Multiply(-alpha, temp1);
                pseudoResiduals.Add(temp1, temp2);
                temp2.CopyTo(pseudoResiduals);

                // d = yOdd + theta * theta * eta / alpha * d
                d.Multiply(theta * theta * eta / alpha, temp);
                yinternal.Add(temp, d);

                // theta = ||pseudoResiduals||_2 / tau
                theta = pseudoResiduals.Norm(2) / tau;
                var c = 1 / Math.Sqrt(1 + (theta * theta));

                // tau = tau * theta * c
                tau *= theta * c;

                // eta = c^2 * alpha
                eta = c * c * alpha;

//.........这里部分代码省略.........
            float sigma;

            // Define the temporary vectors
            // rDash_0 = r_0
            Vector rdash = new DenseVector(residuals);

            // t_-1 = 0
            Vector t = new DenseVector(residuals.Count);
            Vector t0 = new DenseVector(residuals.Count);

            // 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.DotProduct(residuals) / rdash.DotProduct(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.DotProduct(c);

//.........这里部分代码省略.........
            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].DotProduct(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].DotProduct(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.DotProduct(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.DotProduct(temp) / rho;

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

                // Reuse temp
                temp.Multiply(rho, temp);
                residuals.Add(temp, temp2);
                temp2.CopyTo(residuals);

                // x_(jk+1) = x_((j-1)k_k) - rho_(j+1) u~_(jk+1) + alpha_(jk+1) g~_((j-1)k+k)
                utemp.Multiply(-rho, temp);
                xtemp.Add(temp, temp2);
                temp2.CopyTo(xtemp);