Quelle r1updt.h Sprache: C

namespace Eigen {

namespace internal {

template <typename Scalar>
void r1updt(
        Matrix< Scalar, Dynamic, Dynamic > &s,
        const Matrix< Scalar, Dynamic, 1> &u,
        std::vector<JacobiRotation<Scalar> > &v_givens,
        std::vector<JacobiRotation<Scalar> > &w_givens,
        Matrix< Scalar, Dynamic, 1> &v,
        Matrix< Scalar, Dynamic, 1> &w,
        bool *sing)
{
    typedef DenseIndex Index;
    const JacobiRotation<Scalar> IdentityRotation = JacobiRotation<Scalar>(1,0);

    /* Local variables */
    const Index m = s.rows();
    const Index n = s.cols();
    Index i, j=1;
    Scalar temp;
    JacobiRotation<Scalar> givens;

    // r1updt had a broader usecase, but we don't use it here. And, more
    // importantly, we can not test it.
    eigen_assert(m==n);
    eigen_assert(u.size()==m);
    eigen_assert(v.size()==n);
    eigen_assert(w.size()==n);

    /* move the nontrivial part of the last column of s into w. */
    w[n-1] = s(n-1,n-1);

    /* rotate the vector v into a multiple of the n-th unit vector */
    /* in such a way that a spike is introduced into w. */
    for (j=n-2; j>=0; --j) {
        w[j] = 0.;
        if (v[j] != 0.) {
            /* determine a givens rotation which eliminates the */
            /* j-th element of v. */
            givens.makeGivens(-v[n-1], v[j]);

            /* apply the transformation to v and store the information */
            /* necessary to recover the givens rotation. */
            v[n-1] = givens.s() * v[j] + givens.c() * v[n-1];
            v_givens[j] = givens;

            /* apply the transformation to s and extend the spike in w. */
            for (i = j; i < m; ++i) {
                temp = givens.c() * s(j,i) - givens.s() * w[i];
                w[i] = givens.s() * s(j,i) + givens.c() * w[i];
                s(j,i) = temp;
            }
        } else
            v_givens[j] = IdentityRotation;
    }

    /* add the spike from the rank 1 update to w. */
    w += v[n-1] * u;

    /* eliminate the spike. */
    *sing = false;
    for (j = 0; j < n-1; ++j) {
        if (w[j] != 0.) {
            /* determine a givens rotation which eliminates the */
            /* j-th element of the spike. */
            givens.makeGivens(-s(j,j), w[j]);

            /* apply the transformation to s and reduce the spike in w. */
            for (i = j; i < m; ++i) {
                temp = givens.c() * s(j,i) + givens.s() * w[i];
                w[i] = -givens.s() * s(j,i) + givens.c() * w[i];
                s(j,i) = temp;
            }

            /* store the information necessary to recover the */
            /* givens rotation. */
            w_givens[j] = givens;
        } else
            v_givens[j] = IdentityRotation;

        /* test for zero diagonal elements in the output s. */
        if (s(j,j) == 0.) {
            *sing = true;
        }
    }
    /* move w back into the last column of the output s. */
    s(n-1,n-1) = w[n-1];

    if (s(j,j) == 0.) {
        *sing = true;
    }
    return;
}

} // end namespace internal

} // end namespace Eigen

quality68%

¤ Dauer der Verarbeitung: 0.1 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.