fn_polyval.hpp

Go to the documentation of this file.

/*
 * BELFEM -- The Berkeley Lab Finite Element Framework
 * Copyright (c) 2026, The Regents of the University of California,
 * through Lawrence Berkeley National Laboratory (subject to receipt of any required
 * approvals from the U.S. Dept. of Energy).  All rights reserved.
 *
 * Developers: Christian Messe, Gregory Giard
 *
 * See the top-level LICENSE file for the complete license and disclaimer.
 */

 
#ifndef BELFEM_FN_POLYVAL_HPP
#define BELFEM_FN_POLYVAL_HPP
 
#ifdef BELFEM_ARMADILLO
#include "armadillo.hpp"
#endif
 
#include "typedefs.hpp"
#include "assert.h"
#include "cl_Vector.hpp"
 
namespace belfem
{
 
//------------------------------------------------------------------------------

    template < typename T >
    T
    polyval( const Vector< T > & aCoeffs, const T aX )
    {
        const index_t tN = aCoeffs.length();
        T aResult = aCoeffs( 0 );
 
        for( index_t k=1; k<tN; ++k )
        {
            aResult *= aX;
            aResult += aCoeffs( k );
        }
 
        return aResult;
    }
 
//------------------------------------------------------------------------------
 
 
    template < typename T >
    void
    polyval( const Vector< T > & aCoeffs, const Vector< T > & aX, Vector< T > & aY )
    {
#ifdef BELFEM_ARMADILLO
        aY.vector_data() = arma::polyval( aCoeffs.vector_data(), aX.vector_data() );
#elif BELFEM_BLAZE
        index_t tN = aX.length();
        aY.set_size( tN );
 
        if( aCoeffs.length() == 0 )
        {
            aY.vector_data() = 0;
            return;
        }
 
        // Initialize with highest order coefficient
        aY.vector_data() = aCoeffs(0);
 
        // Horner's method with vectorized operations
        for( index_t i = 1; i < aCoeffs.length(); ++i )
        {
            aY.vector_data() = aY.vector_data() * aX.vector_data() + aCoeffs(i);
        }
#else
        index_t tN = aX.vector_data();
        aY.set_size( tN );
        for( uint k=0; k<tN; ++k )
        {
            aY( k ) = polyval( aCoeffs, aX( k ) );
        }
#endif
    }
 
 
 
//------------------------------------------------------------------------------
}
#endif //BELFEM_FN_POLYVAL_HPP