trustregiongsstep.cc

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:

template<class MatrixType, class VectorType>
inline
void TrustRegionGSStep<MatrixType, VectorType>::iterate()
{
    const MatrixType& mat = *this->mat_;
    VectorType& x = *this->x_;
    const std::vector<BoxConstraint<Field,BlockSize> >& obstacles = *this->obstacles_;

    for (size_t i=0; i<x.size(); i++) {

        // Dirichlet nodes in this block
        std::bitset<BlockSize> ignoreNodes(0);
        if (this->hasIgnore())
            ignoreNodes = this->ignore()[i];

        if (ignoreNodes.all())
            continue;

        VectorBlock blockResidual;
        this->residual(i, blockResidual);

        mat[i][i].umv(x[i], blockResidual);

        // scalar gauss-seidel
        for (size_t j=0; j<BlockSize; j++) {

            if (ignoreNodes[j])
                continue;

            // Compute residual
            Field r = blockResidual[j] - mat[i][i][j] * x[i];
            r += mat[i][i][j][j] * x[i][j];

            const Field& diag = mat[i][i][j][j];

            // Find local minimum
            if (diag > 0) {
                // 1d problem is convex
                x[i][j] = r / diag;

                x[i][j] = obstacles[i][j].projectIn(x[i][j]);

            } else  {

                // If the corresponding dof was truncated, then compute no correction
                if (diag>-1e-13 && std::fabs(r)<1e-13)
                    continue;

                // 1d problem is concave or linear.
                // Minimum is attained at one of the boundaries
                Field lBound = obstacles[i].lower(j);
                Field uBound = obstacles[i].upper(j);

                Field lValue = 0.5*diag*lBound*lBound - r*lBound;
                Field uValue = 0.5*diag*uBound*uBound - r*uBound;

                x[i][j] = (lValue < uValue) ? lBound : uBound;

            }

        }

    }

}