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 <complex>
 
#include <graphblas.hpp>
#include <graphblas/utils/iscomplex.hpp>
 
 
namespace grb {
 
    namespace algorithms {
 
        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 >
        >
        grb::RC conjugate_gradient(
            grb::Vector< IOType > &x,
            const grb::Matrix< NonzeroType > &A,
            const grb::Vector< InputType > &b,
            const size_t max_iterations,
            ResidualType tol,
            size_t &iterations,
            ResidualType &residual,
            grb::Vector< IOType > &r,
            grb::Vector< IOType > &u,
            grb::Vector< IOType > &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."
            );
 
            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 );
 
            // dynamic checks
            {
                const size_t m = grb::nrows( A );
                if( size( x ) != n ) {
                    return MISMATCH;
                }
                if( size( b ) != m ) {
                    return MISMATCH;
                }
                if( size( r ) != n || size( u ) != n || size( temp ) != n ) {
                    std::cerr << "Error: provided workspace vectors are not of the correct "
                        << "length.\n";
                    return 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 ILLEGAL;
                }
 
                // capacities
                if( capacity( x ) != n ) {
                    return ILLEGAL;
                }
                if( capacity( r ) != n || capacity( u ) != n || capacity( temp ) != n ) {
                    return ILLEGAL;
                }
 
                // others
                if( tol <= zero_residual ) {
                    std::cerr << "Error: tolerance input to CG must be strictly positive\n";
                    return ILLEGAL;
                }
            }
 
            // set pure output fields to neutral defaults
            iterations = 0;
            residual = std::numeric_limits< double >::infinity();
 
            // trivial shortcuts
            if( max_iterations == 0 ) {
                return FAILED;
            }
 
            // 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 = SUCCESS;
                if( nnz( x ) != n ) {
                    rc = set< descriptors::invert_mask | descriptors::structural >(
                        x, x, zero
                    );
                }
                if( rc != SUCCESS ) {
                    return rc;
                }
                assert( nnz( x ) == n );
            }
 
            IOType sigma, bnorm, alpha, beta;
 
            // temp = 0
            grb::RC ret = grb::set( temp, 0 );
            assert( ret == SUCCESS );
 
            // temp = A * x
            ret = ret ? ret : grb::mxv< descr_dense >( temp, A, x, ring );
            assert( ret == SUCCESS );
 
            // r = b - temp;
            ret = ret ? ret : grb::set( r, zero );
            ret = ret ? ret : grb::foldl( r, b, ring.getAdditiveMonoid() );
            assert( nnz( r ) == n );
            assert( nnz( temp ) == n );
            ret = ret ? ret : grb::foldl< descr_dense >( r, temp, minus );
            assert( ret == SUCCESS );
            assert( nnz( r ) == n );
 
            // u = r;
            ret = ret ? ret : grb::set( u, r );
            assert( ret == SUCCESS );
 
            // sigma = r' * r;
            sigma = zero;
            if( grb::utils::is_complex< IOType >::value ) {
                ret = ret ? ret : grb::eWiseLambda( [&temp,&r]( const size_t i ) {
                        temp[ i ] = grb::utils::is_complex< IOType >::conjugate( r[ i ] );
                    }, temp
                );
                ret = ret ? ret : grb::dot< descr_dense >( sigma, temp, r, ring );
            } else {
                ret = ret ? ret : grb::dot< descr_dense >( sigma, r, r, ring );
            }
 
            assert( ret == SUCCESS );
 
            // bnorm = b' * b;
            bnorm = zero;
            if( grb::utils::is_complex< IOType >::value ) {
                ret = ret ? ret : grb::eWiseLambda( [&temp,&b]( const size_t i ) {
                        temp[ i ] = grb::utils::is_complex< IOType >::conjugate( b[ i ] );
                    }, temp
                );
                ret = ret ? ret : grb::dot< descr_dense >( bnorm, temp, b, ring );
            } else {
                ret = ret ? ret : grb::dot< descr_dense >( bnorm, b, b, ring );
            }
            assert( ret == SUCCESS );
 
            if( ret == SUCCESS ) {
                tol *= sqrt( grb::utils::is_complex< IOType >::modulus( bnorm ) );
            }
 
            size_t iter = 0;
 
            do {
                assert( iter < max_iterations );
                (void) ++iter;
 
                // temp = 0
                ret = ret ? ret : grb::set( temp, 0 );
                assert( ret == SUCCESS );
 
                // temp = A * u;
                ret = ret ? ret : grb::mxv< descr_dense >( temp, A, u, ring );
                assert( ret == SUCCESS );
 
                // beta = u' * temp
                beta = zero;
                if( grb::utils::is_complex< IOType >::value ) {
                    ret = ret ? ret : grb::eWiseLambda( [&u]( const size_t i ) {
                            u[ i ] = grb::utils::is_complex< IOType >::conjugate( u[ i ] );
                        }, u
                    );
                }
                ret = ret ? ret : grb::dot< descr_dense >( beta, temp, u, ring );
                if( grb::utils::is_complex< IOType >::value ) {
                    ret = ret ? ret : grb::eWiseLambda( [&u]( const size_t i ) {
                            u[ i ] = grb::utils::is_complex< IOType >::conjugate( u[ i ] );
                        }, u
                    );
                }
                assert( ret == SUCCESS );
 
                // alpha = sigma / beta;
                ret = ret ? ret : grb::apply( alpha, sigma, beta, divide );
                assert( ret == SUCCESS );
 
                // x = x + alpha * u;
                ret = ret ? ret : grb::eWiseMul< descr_dense >( x, alpha, u, ring );
                assert( ret == SUCCESS );
 
                // temp = alpha .* temp
                // Warning: operator-based foldr requires temp be dense
                ret = ret ? ret : grb::foldr( alpha, temp, ring.getMultiplicativeMonoid() );
                assert( ret == SUCCESS );
 
                // r = r - temp;
                ret = ret ? ret : grb::foldl< descr_dense >( r, temp, minus );
                assert( ret == SUCCESS );
 
                // beta = r' * r;
                beta = zero;
                if( grb::utils::is_complex< IOType >::value ) {
                    ret = ret ? ret : grb::eWiseLambda( [&temp,&r]( const size_t i ) {
                            temp[ i ] = grb::utils::is_complex< IOType >::conjugate( r[ i ] );
                        }, temp
                    );
                    ret = ret ? ret : grb::dot< descr_dense >( beta, temp, r, ring );
                } else {
                    ret = ret ? ret : grb::dot< descr_dense >( beta, r, r, ring );
                }
                residual = grb::utils::is_complex< IOType >::modulus( beta );
                assert( ret == SUCCESS );
 
                if( ret == SUCCESS ) {
                    if( sqrt( residual ) < tol || iter >= max_iterations ) {
                        break;
                    }
                }
 
                // alpha = beta / sigma;
                ret = ret ? ret : grb::apply( alpha, beta, sigma, divide );
                assert( ret == SUCCESS );
 
                // temp = r + alpha * u;
                ret = ret ? ret : grb::set( temp, r );
                assert( ret == SUCCESS );
                ret = ret ? ret : grb::eWiseMul< descr_dense >( temp, alpha, u, ring );
                assert( ret == SUCCESS );
                assert( nnz( temp ) == size( temp ) );
 
                // u = temp
                std::swap( u, temp );
 
                sigma = beta;
            } while( ret == SUCCESS );
 
            // output that is independent of error code
            iterations = iter;
 
            // return correct error code
            if( ret == SUCCESS ) {
                if( sqrt( residual ) >= tol ) {
                    // did not converge within iterations
                    return FAILED;
                }
            }
            return ret;
        }
 
    } // namespace algorithms
 
} // end namespace grb
 
#endif // end _H_GRB_ALGORITHMS_CONJUGATE_GRADIENT