conjugate_gradient.hpp

Go to the documentation of this file.

/*
 *   Copyright 2021 Huawei Technologies Co., Ltd.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef _H_GRB_ALGORITHMS_CONJUGATE_GRADIENT
#define _H_GRB_ALGORITHMS_CONJUGATE_GRADIENT

#include <cstdio>
#include <cmath>

#include <graphblas.hpp>
#include <graphblas/utils/iscomplex.hpp>


namespace grb {

    namespace algorithms {

        template<
            Descriptor descr = descriptors::no_operation,
            bool preconditioned = true,
            typename IOType,
            typename ResidualType,
            typename NonzeroType,
            typename InputType,
            class Ring = Semiring<
                grb::operators::add< IOType >, grb::operators::mul< IOType >,
                grb::identities::zero, grb::identities::one
            >,
            class Minus = operators::subtract< IOType >,
            class Divide = operators::divide< IOType >,
            typename RSI, typename NZI, Backend backend
        >
        grb::RC preconditioned_conjugate_gradient(
            grb::Vector< IOType, backend > &x,
            const grb::Matrix< NonzeroType, backend, RSI, RSI, NZI > &A,
            const grb::Vector< InputType, backend > &b,
            const std::function<
                    grb::RC(
                        grb::Vector< IOType, backend > &,
                        const grb::Vector< IOType, backend > &
                        )
                    > &Minv,
            const size_t max_iterations,
            ResidualType tol,
            size_t &iterations,
            ResidualType &residual,
            grb::Vector< IOType, backend > &r,
            grb::Vector< IOType, backend > &u,
            grb::Vector< IOType, backend > &temp,
            grb::Vector< IOType, backend > &temp_precond,
            const Ring &ring = Ring(),
            const Minus &minus = Minus(),
            const Divide &divide = Divide()
        ) {
            // static checks
            static_assert( std::is_floating_point< ResidualType >::value,
                "Can only use the CG algorithm with floating-point residual "
                "types." ); // unless some different norm were used: issue #89
            static_assert( !( descr & descriptors::no_casting ) || (
                    std::is_same< IOType, ResidualType >::value &&
                    std::is_same< IOType, NonzeroType >::value &&
                    std::is_same< IOType, InputType >::value
                ), "One or more of the provided containers have differing element types "
                "while the no-casting descriptor has been supplied"
            );
            static_assert( !( descr & descriptors::no_casting ) || (
                    std::is_same< NonzeroType, typename Ring::D1 >::value &&
                    std::is_same< IOType, typename Ring::D2 >::value &&
                    std::is_same< InputType, typename Ring::D3 >::value &&
                    std::is_same< InputType, typename Ring::D4 >::value
                ), "no_casting descriptor was set, but semiring has incompatible domains "
                "with the given containers."
            );
            static_assert( !( descr & descriptors::no_casting ) || (
                    std::is_same< InputType, typename Minus::D1 >::value &&
                    std::is_same< InputType, typename Minus::D2 >::value &&
                    std::is_same< InputType, typename Minus::D3 >::value
                ), "no_casting descriptor was set, but given minus operator has "
                "incompatible domains with the given containers."
            );
            static_assert( !( descr & descriptors::no_casting ) || (
                    std::is_same< ResidualType, typename Divide::D1 >::value &&
                    std::is_same< ResidualType, typename Divide::D2 >::value &&
                    std::is_same< ResidualType, typename Divide::D3 >::value
                ), "no_casting descriptor was set, but given divide operator has "
                "incompatible domains with the given tolerance type."
            );
            static_assert( std::is_floating_point< ResidualType >::value,
                "Require floating-point residual type."
            );

            constexpr const Descriptor descr_dense = descr | descriptors::dense;
            const ResidualType zero_residual = ring.template getZero< ResidualType >();
            const IOType zero = ring.template getZero< IOType >();
            const size_t n = grb::ncols( A );

            // retrieve conditional buffers
            grb::Vector< IOType, backend > &z = preconditioned ? temp_precond : r;

            // dynamic checks
            {
                const size_t m = grb::nrows( A );
                if( grb::size( x ) != n ) {
                    std::cerr << "Error: initial solution guess and output vector length ("
                        << size( x ) << ") does not match matrix size (" << m << ").\n";
                    return grb::MISMATCH;
                }
                if( grb::size( b ) != m ) {
                    std::cerr << "Error: right-hand side size (" << grb::size( b ) << ") does "
                        << "not match matrix size (" << m << ").\n";
                    return grb::MISMATCH;
                }
                if( grb::size( r ) != n || grb::size( u ) != n || grb::size( temp ) != n ) {
                    std::cerr << "Error: provided workspace vectors are not of the correct "
                        << "length.\n";
                    return grb::MISMATCH;
                }
                if( preconditioned && grb::size( temp_precond ) != n ) {
                    std::cerr << "Error: (left) preconditioner workspace vector does not have "
                        << "the correct length.\n";
                    return grb::MISMATCH;
                }
                if( m != n ) {
                    std::cerr << "Warning: grb::algorithms::conjugate_gradient requires "
                        << "square input matrices, but a non-square input matrix was "
                        << "given instead.\n";
                    return grb::ILLEGAL;
                }

                // capacities
                if( grb::capacity( x ) != n ) {
                    return grb::ILLEGAL;
                }
                if( grb::capacity( r ) != n || grb::capacity( u ) != n ||
                    grb::capacity( temp ) != n
                ) {
                    return grb::ILLEGAL;
                }
                if( preconditioned && grb::capacity( temp_precond ) != n ) {
                    return grb::ILLEGAL;
                }

                // others
                if( tol <= zero_residual ) {
                    std::cerr << "Error: tolerance input to CG must be strictly positive\n";
                    return grb::ILLEGAL;
                }
                if( max_iterations == 0 ) {
                    std::cerr << "Error: at least one CG iteration must be requested\n";
                    return grb::ILLEGAL;
                }
            }

            // set pure output fields to neutral defaults
            iterations = 0;
            residual = std::numeric_limits< double >::infinity();

            // make x and b structurally dense (if not already) so that the remainder
            // algorithm can safely use the dense descriptor for faster operations
            {
                RC rc = grb::SUCCESS;
                if( nnz( x ) != n ) {
                    rc = grb::set< descriptors::invert_mask | descriptors::structural >(
                        x, x, zero
                    );
                }
                if( rc != grb::SUCCESS ) { return rc; }
                assert( nnz( x ) == n );
            }

            IOType sigma, bnorm, alpha, beta;

            // r = b - temp;
            grb::RC ret = grb::set( temp, 0 ); assert( ret == grb::SUCCESS );
            ret = ret ? ret : grb::mxv< descr_dense >( temp, A, x, ring );
            assert( ret == grb::SUCCESS );
            ret = ret ? ret : grb::set( r, zero ); assert( ret == grb::SUCCESS );
            // note: no dense descriptor since we actually allow sparse b
            ret = ret ? ret : grb::foldl( r, b, ring.getAdditiveMonoid() );
            // from here onwards, r, temp, x are dense and will remain so
            assert( nnz( r ) == n );
            assert( nnz( temp ) == n );
            ret = ret ? ret : grb::foldl< descr_dense >( r, temp, minus );
            assert( ret == grb::SUCCESS );

            // z = M^-1r
            if( preconditioned ) {
                ret = ret ? ret : grb::set( z, 0 ); // also ensures z is dense, henceforth
                ret = ret ? ret : Minv( z, r );
            } // else, z equals r (by reference)

            // u = z;
            ret = ret ? ret : grb::set( u, z );
            assert( ret == grb::SUCCESS );
            // from here onwards, u is dense; i.e., all vectors are dense from now on,
            // and we can freely use the dense descriptor in the subsequent

            // sigma = r' * z;
            sigma = zero;
            ret = ret ? ret : grb::dot< descr_dense >(
                    sigma,
                    r, z,
                    ring.getAdditiveMonoid(),
                    grb::operators::conjugate_right_mul< IOType >()
                );

            assert( ret == grb::SUCCESS );

            // bnorm = b' * b;
            bnorm = zero;
            ret = ret ? ret : grb::dot< descr_dense >(
                bnorm,
                b, b,
                ring.getAdditiveMonoid(),
                grb::operators::conjugate_left_mul< IOType >() );
            assert( ret == grb::SUCCESS );

            // get effective tolerance and exit on any error during prelude
            if( ret == grb::SUCCESS ) {
                tol *= std::sqrt( grb::utils::is_complex< IOType >::modulus( bnorm ) );
            } else {
                std::cerr << "Warning: preconditioned CG caught error during prelude ("
                    << grb::toString( ret ) << ")\n";
                return ret;
            }

            // all OK; perform main iterations
            size_t iter = 0;
            do {
                assert( iter < max_iterations );
                (void) ++iter;

                // temp = A * u;
                ret = ret ? ret : grb::set< descr_dense >( temp, 0 );
                assert( ret == grb::SUCCESS );
                ret = ret ? ret : grb::mxv< descr_dense >( temp, A, u, ring );
                assert( ret == grb::SUCCESS );

                // beta = (A * u)' * u;
                beta = zero;
                ret = ret ? ret : grb::dot< descr_dense >(
                        beta,
                        temp, u,
                        ring.getAdditiveMonoid(),
                        grb::operators::conjugate_right_mul< IOType >()
                    );
                assert( ret == grb::SUCCESS );

                // alpha = sigma / beta;
                ret = ret ? ret : grb::apply( alpha, sigma, beta, divide );
                assert( ret == grb::SUCCESS );

                // x = x + alpha * u;
                ret = ret ? ret : grb::eWiseMul< descr_dense >( x, alpha, u, ring );
                assert( ret == grb::SUCCESS );

                // r = r - alpha .* temp = r - alpha .* (A * u);
                ret = ret ? ret : grb::foldr< descr_dense >( alpha, temp,
                    ring.getMultiplicativeMonoid() );
                assert( ret == grb::SUCCESS );
                ret = ret ? ret : grb::foldl< descr_dense >( r, temp, minus );
                assert( ret == grb::SUCCESS );

                // get residual. In the preconditioned case, the resulting scalar is *not*
                // used for subsequent operations. Therefore, we first compute the residual
                // using alpha as a temporary scalar
                alpha = zero;
                ret = ret ? ret : grb::dot< descr_dense >(
                        alpha,
                        r, r,
                        ring.getAdditiveMonoid(),
                        grb::operators::conjugate_left_mul< IOType >()
                    );
                assert( ret == grb::SUCCESS );
                residual = grb::utils::is_complex< IOType >::modulus( alpha );

                // check residual
                if( ret == grb::SUCCESS ) {
                    if( sqrt( residual ) < tol || iter >= max_iterations ) { break; }
                }

                // apply preconditioner action (if required), and compute beta for the
                // preconditioned case
                // z = M^-1r
                // beta = r' * z
                if( preconditioned ) {
                    beta = zero;
                    ret = ret ? ret : Minv( z, r ); assert( ret == grb::SUCCESS );
                    ret = ret ? ret : grb::dot< descr_dense >(
                            beta,
                            r, z,
                            ring.getAdditiveMonoid(),
                            grb::operators::conjugate_right_mul< IOType >()
                        );
                    assert( ret == grb::SUCCESS );
                } else {
                    beta = alpha;
                }

                // alpha = beta / sigma;
                ret = ret ? ret : grb::apply( alpha, beta, sigma, divide );
                assert( ret == grb::SUCCESS );

                // u_next = z + beta * u_previous;
                ret = ret ? ret : grb::foldr< descr_dense >( alpha, u,
                    ring.getMultiplicativeMonoid() );
                assert( ret == grb::SUCCESS );
                ret = ret ? ret : grb::foldr< descr_dense >( z, u,
                    ring.getAdditiveMonoid() );
                assert( ret == grb::SUCCESS );

                sigma = beta;
            } while( ret == grb::SUCCESS );

            // output that is independent of error code
            iterations = iter;

            // return correct error code
            if( ret == grb::SUCCESS ) {
                if( std::sqrt( residual ) >= tol ) {
                    // did not converge within iterations
                    return grb::FAILED;
                }
            }
            return ret;
        }

        template<
            Descriptor descr = descriptors::no_operation,
            typename IOType,
            typename ResidualType,
            typename NonzeroType,
            typename InputType,
            class Ring = Semiring<
                grb::operators::add< IOType >, grb::operators::mul< IOType >,
                grb::identities::zero, grb::identities::one
            >,
            class Minus = operators::subtract< IOType >,
            class Divide = operators::divide< IOType >,
            typename RSI, typename NZI, Backend backend
        >
        grb::RC conjugate_gradient(
            grb::Vector< IOType, backend > &x,
            const grb::Matrix< NonzeroType, backend, RSI, RSI, NZI > &A,
            const grb::Vector< InputType, backend > &b,
            const size_t max_iterations,
            ResidualType tol,
            size_t &iterations,
            ResidualType &residual,
            grb::Vector< IOType, backend > &r,
            grb::Vector< IOType, backend > &u,
            grb::Vector< IOType, backend > &temp,
            const Ring &ring = Ring(),
            const Minus &minus = Minus(),
            const Divide &divide = Divide()
        ) {
            // static checks
            static_assert( std::is_floating_point< ResidualType >::value,
                "Can only use the CG algorithm with floating-point residual "
                "types." ); // unless some different norm were used: issue #89
            static_assert( !( descr & descriptors::no_casting ) || (
                    std::is_same< IOType, ResidualType >::value &&
                    std::is_same< IOType, NonzeroType >::value &&
                    std::is_same< IOType, InputType >::value
                ), "One or more of the provided containers have differing element types "
                "while the no-casting descriptor has been supplied"
            );
            static_assert( !( descr & descriptors::no_casting ) || (
                    std::is_same< NonzeroType, typename Ring::D1 >::value &&
                    std::is_same< IOType, typename Ring::D2 >::value &&
                    std::is_same< InputType, typename Ring::D3 >::value &&
                    std::is_same< InputType, typename Ring::D4 >::value
                ), "no_casting descriptor was set, but semiring has incompatible domains "
                "with the given containers."
            );
            static_assert( !( descr & descriptors::no_casting ) || (
                    std::is_same< InputType, typename Minus::D1 >::value &&
                    std::is_same< InputType, typename Minus::D2 >::value &&
                    std::is_same< InputType, typename Minus::D3 >::value
                ), "no_casting descriptor was set, but given minus operator has "
                "incompatible domains with the given containers."
            );
            static_assert( !( descr & descriptors::no_casting ) || (
                    std::is_same< ResidualType, typename Divide::D1 >::value &&
                    std::is_same< ResidualType, typename Divide::D2 >::value &&
                    std::is_same< ResidualType, typename Divide::D3 >::value
                ), "no_casting descriptor was set, but given divide operator has "
                "incompatible domains with the given tolerance type."
            );
            static_assert( std::is_floating_point< ResidualType >::value,
                "Require floating-point residual type."
            );

            // create a dummy preconditioner and buffer that will never be used
            std::function<
                grb::RC(
                    grb::Vector< IOType, backend >&,
                    const grb::Vector< IOType, backend >&
                )
            > dummy_preconditioner =
                [](
                    grb::Vector< IOType, backend >&,
                    const grb::Vector< IOType, backend >&
                ) -> grb::RC {
                    return grb::FAILED;
                };
            grb::Vector< IOType, backend > dummy_buffer( 0 );

            // call PCG with preconditioning disabled
            return preconditioned_conjugate_gradient< descr, false >(
                    x, A, b,
                    dummy_preconditioner,
                    max_iterations, tol,
                    iterations, residual,
                    r, u, temp, dummy_buffer,
                    ring, minus, divide
                );
        }

    } // namespace algorithms

} // end namespace grb

#endif // end _H_GRB_ALGORITHMS_CONJUGATE_GRADIENT