首页 > 其他分享 >高性能计算-Intel IPP库ARM移植示例(20)

高性能计算-Intel IPP库ARM移植示例(20)

时间：2024-12-17 17:44:10浏览次数：7  
标签：20 Intel 示例 IPP vx f32 result 2PI vresultIM

1. 简介

Intel® Integrated Performance Primitives，即英特尔集成性能基元（简称IPP），为信号、数据和图像处理特定应用领域，提供simd优化的一组全面的函数库。

本项目将对 exp、cos、sin、tone、Triangle函数用NEON向量化指令实现ARM移植版本，有串行和向量化两个版本。

计算使用泰勒展开，循环加和保证接口要求的精度即可。

移植接口参数可在 IPP 官方文档查看详情，保持一致即可。

项目地址：https://github.com/libo-0379/IPP_ARM

2. 核心代码

#include <stdio.h>
#include <math.h>
#include "../include/ipp.h"

IppStatus ippsExp_32f_A24(const Ipp32f *pSrc, Ipp32f *dst, Ipp32s len)
{
    for(int i=0;i<len/4*4;i+=4)
    {   
        //每个元素泰勒展开项计数，初始化为0
        float32x4_t vn = vdupq_n_f32(0);
        float32x4_t vx = vld1q_f32(pSrc+i);
        //记录当前项的展开值，初始化为第一项的值
        float32x4_t vprior = vdupq_n_f32(1.0);
        //保存每个元素展开的求和
        float32x4_t vsum = vprior;
        while(1)
        {
            //计算当前展开项
            // 1. 计算 prior * x / ++n
            vn = vaddq_f32(vn,vdupq_n_f32(1));
            // 自己实现的 vdivq_f32_neon 对 -20.055 计算误差太大，改用系统提供接口 vdivq_f32
            // float32x4_t vcur = vdivq_f32_neon(vmulq_f32(vprior,vx),vn);
            float32x4_t vcur = vdivq_f32(vmulq_f32(vprior,vx),vn);
            // 2. 将当前项赋值为前项            
            vprior = vcur;
            // printf("0 %f\n",vgetq_lane_f32(vcur,0));
            // 3. 如果所有元素的最大展开项小于精度控制，不再计算
            if(vmaxvq_f32(vabsq_f32(vcur)) - EPSILON < 0)
                break;
            // 4. 将当前展开项加和
            vsum = vaddq_f32(vsum,vcur);
        }
        // printf("n %f\n",vgetq_lane_f32(vn,0));
        vst1q_f32(dst+i,vsum);
    }

    //处理尾部数据
    for(int i=len/4*4;i<len;i++)
    {
        // 初始化第一个值
        int n = 0;
        double prior = 1.0; //前一项数值初始化为第一项的数值
        double sum = prior; //求和保存结果
        Ipp32f x = *(pSrc+i);
        while (1)
        {
            double cur = prior * x / ++n;   //当前加项的数值 = 前项*(src[i]/n)
            prior = cur;
            // printf("%f\n",cur);
            if (fabs(cur)-EPSILON <= 0) //如果小于精度值停止计算
                break;
            sum += cur;
        }
        // printf("n %d\n",n);
        *(dst+i) = sum;
    }
    return ippStsNoErr;
}

IppStatus ippsCos_32f_A24(const Ipp32f *pSrc, Ipp32f *dst, Ipp32s len)
{
    float32x4_t v2PI = vdupq_n_f32(IPP_2PI);
    for(int i=0;i<len/4*4;i+=4)
    {
        //每个元素泰勒展开项计数，初始化为0
        float32x4_t vn = vdupq_n_f32(0);
        float32x4_t vx = vld1q_f32(pSrc+i);
        //记录当前项的展开值，初始化为第一项的值
        float32x4_t vprior = vdupq_n_f32(1.0);
        //保存每个元素展开的求和
        float32x4_t vsum = vprior;
        //将x转换到 [-2π,2π]
        //类型转换向下取整,注意这里是浮点数除法
        float32x4_t vk = vcvtq_f32_s32(vcvtq_s32_f32(vdivq_f32_neon(vx,v2PI)));   
        vx = vfmsq_f32(vx,v2PI,vk);
        while(1)
        {
            //计算当前展开项
            // 1. 计算 prior * (x^2) / (2(++n)*(2n-1))
            vn = vaddq_f32(vn,vdupq_n_f32(1));
            float32x4_t vcur = vmulq_n_f32(vn,2);
            vcur = vdivq_f32_neon(vmulq_f32(vprior,vmulq_f32(vx,vx)),vmulq_f32(vcur,vsubq_f32(vcur,vdupq_n_f32(1.0))));
            //将当前项赋值为前项,不带符号          
            vprior = vcur;
            // 2. 添加正负号
            float32_t n = vgetq_lane_f32(vn,0);
            vcur = vmulq_n_f32(vcur,(int)n%2==0?1:-1);
            // 3. 如果所有元素的展开项绝对值最大项小于精度控制，不再计算
            if(vmaxvq_f32(vabsq_f32(vcur)) - EPSILON <= 0)
                break;
            // 4. 将当前展开项加和
            vsum = vaddq_f32(vsum,vcur);
        }
        printf("n %f\n",vgetq_lane_f32(vn,0));
        vst1q_f32(dst+i,vsum);
    }
    //处理尾部数据
    for(int i=len/4*4;i<len;i++)
    {
        // 初始化第一个值
        int n = 0;
        double prior = 1.0; //前一项数值初始化为第一项的数值
        double sum = prior; //求和保存结果
        Ipp64f x = *(pSrc+i);
        //对x转换到 [-2π,2π]
        Ipp32s k = x/IPP_2PI;
        x -= k*IPP_2PI;
        while (1)
        {
            int32_t temp1 = 2* ++n;
            double cur = prior * (x * x) /(temp1-- * temp1) ;   //当前加项的数值 = 前项*(x^2/(2n*(2n-1)))
            prior = cur;
            if (fabs(cur) - EPSILON <= 0) //如果小于精度值停止计算
                break;
            sum += (n%2==0?1:-1)*cur;   //符号位
        }
        // printf("n %d\n",n);
        *(dst+i) = sum;
    }
    return ippStsNoErr;
}

IppStatus ippsSin_32f_A24(const Ipp32f *pSrc, Ipp32f *dst, Ipp32s len)
{    
    float32x4_t v2PI = vdupq_n_f32(IPP_2PI);
    for(int i=0;i<len/4*4;i+=4)
    {
        //每个元素泰勒展开项计数，初始化为0
        float32x4_t vn = vdupq_n_f32(0);
        float32x4_t vx = vld1q_f32(pSrc+i);
        //将x转换到 [-2π,2π]
        //类型转换向下取整,注意这里是浮点数除法
        float32x4_t vk = vcvtq_f32_s32(vcvtq_s32_f32(vdivq_f32_neon(vx,v2PI)));   
        vx = vfmsq_f32(vx,v2PI,vk);
        // !! 注意这里易错点，应在区间转换后再赋值给 vprior 
        //记录当前项的展开值，初始化为第一项的值
        float32x4_t vprior = vx;
        //保存每个元素展开的求和
        float32x4_t vsum = vprior;
        printf("x0 %f\n",vgetq_lane_f32(vx,0));
        while(1)
        {
            // 1. 计算 prior * (x^2) / (2(++n)*(2n+1))
            vn = vaddq_f32(vn,vdupq_n_f32(1));
            float32x4_t vcur = vmulq_n_f32(vn,2);
            vcur = vdivq_f32_neon(vmulq_f32(vprior,vmulq_f32(vx,vx)),vmulq_f32(vcur,vaddq_f32(vcur,vdupq_n_f32(1))));
            // printf("cur0 %f\n",vgetq_lane_f32(vcur,0));
            //将当前项赋值为前项,不带符号          
            vprior = vcur;
            // 2. 添加正负号
            float32_t n = vgetq_lane_f32(vn,0);
            vcur = vmulq_n_f32(vcur,(int)n%2==0?1:-1);
            // 3.如果所有元素的展开项绝对值最大项小于精度控制，不再计算
            if(vmaxvq_f32(vabsq_f32(vcur)) - EPSILON <= 0)
                break;
            // 4.将当前展开项加和
            vsum = vaddq_f32(vsum,vcur);
        }
        printf("n %f\n",vgetq_lane_f32(vn,0));
        vst1q_f32(dst+i,vsum);
    }
    //处理尾部数据
    for(int i=len/4*4;i<len;i++)
    {
        // 初始化第一个值
        int n = 0;
        Ipp64f x = *(pSrc+i);
        //将x转换到 [-2π,2π]
        Ipp32s k = x/IPP_2PI;
        x -= k*IPP_2PI;
        double prior = x; //前一项数值初始化为第一项的数值
        double sum = prior; //求和保存结果
        while (1)
        {
            int32_t temp1 = 2* ++n;
            double cur = prior * (x * x) /(temp1++ * temp1) ;   //当前加项的数值 = 前项*(x^2/(2n*(2n+1)))
            prior = cur;
            if (fabs(cur) - EPSILON <= 0) //如果小于精度值停止计算
                break;
            sum += (n%2==0?1:-1)*cur;   //符号位
        }
        // printf("n %d\n",n);
        *(dst+i) = sum;
    }
    return ippStsNoErr;
}

IppStatus ippsTone_32f(Ipp32f* pDst, int len, Ipp32f magn, Ipp32f rFreq, Ipp32f*
pPhase, IppHintAlgorithm hint)
{
    if(rFreq<0.0 || (rFreq-0.5)>EPSILON)
        return ippStsToneFreqErr;
    if(*pPhase<0.0 || (*pPhase-IPP_2PI)>=EPSILON)
        return ippStsTonePhaseErr;
    if(magn<=0.0)
        return ippStsToneMagnErr;
    
    float32x4_t vphase = vdupq_n_f32(*pPhase);
    float32x4_t vrFreq= vdupq_n_f32(rFreq);
    for(int i=0;i<len/4*4;i+=4)
    {
        //保存结果
        float32_t result[4] = {0.0};
        float32x4_t vresult = {0.0};
        // IPP_2PI*rFreq*i + phase
        float32_t x[4] = {i,i+1,i+2,i+3};
        float32x4_t vx = vld1q_f32(x);
        vx = vaddq_f32(vphase,vmulq_n_f32(vmulq_f32(vrFreq,vx),IPP_2PI));
        vst1q_f32(x,vx);
        ippsCos_32f_A24(x,result,4);
        vresult = vld1q_f32(result);
        vst1q_f32(pDst+i,vmulq_n_f32(vresult,magn));
    }
    for(int i=len/4*4;i<len;i++)
    {
        float32_t result = 0.0;
        Ipp32f x = IPP_2PI*rFreq*i + (*pPhase);
        ippsCos_32f_A24_s(&x,&result,1);
        *(pDst+i) = result * magn;  //类型强转
    }
    return ippStsNoErr;
}

IppStatus ippsTriangle_32fc(Ipp32fc* pDst, int len, Ipp32f magn, Ipp32f rFreq, Ipp32f
asym, Ipp32f* pPhase)
{
    if(rFreq<0.0 || (rFreq-0.5)>EPSILON)
        return ippStsToneFreqErr;
    if(*pPhase<0.0 || (*pPhase-IPP_2PI)>=EPSILON)
        return ippStsTonePhaseErr;
    if(magn<=0.0)
        return ippStsToneMagnErr;
    
    float32x4_t vH = vdupq_n_f32(IPP_PI+asym);
    printf("h %f\n",IPP_PI+asym);

    float32x4_t vPI = vdupq_n_f32(IPP_PI);
    float32x4_t v2PI = vdupq_n_f32(IPP_2PI);
    float32x4_t v2PI_H = vsubq_f32(v2PI,vH);
    float32x4_t vphase = vdupq_n_f32(*pPhase);
    float32x4_t vrFreq= vdupq_n_f32(rFreq);
    for(int i=0;i<len/4*4;i+=4)
    {
        float32_t x[4] = {i,i+1,i+2,i+3};
        float32x4_t vx = vld1q_f32(x);
        // 计算位置 x = 2*PI*rFreq*n + phase
        vx = vaddq_f32(vphase,vmulq_n_f32(vmulq_f32(vrFreq,vx),IPP_2PI));
        //因为x>0 可以转换到 [0,2π]
        //类型转换向下取整,注意这里是浮点数除法
        float32x4_t vk = vcvtq_f32_s32(vcvtq_s32_f32(vdivq_f32_neon(vx,v2PI)));   
        vx = vfmsq_f32(vx,v2PI,vk);

        printf("x %f\n",vgetq_lane_f32(vx,0));
        
        vst1q_f32(x,vx);
        // for(int i=0;i<4;i++)
        //     printf("x %f\n",x[i]);
        //保存结果包含实部虚部
        float32x4x2_t vresult;
        //实部
        float32x4_t vresultRE;
        float32x4_t vresultRE1;
        float32x4_t vresultRE2;
        //虚部
        float32x4_t vresultIM;
        float32x4_t vresultIM_1;
        float32x4_t vresultIM1;
        float32x4_t vresultIM2;
        float32x4_t vresultIM3;
        //对实部和虚部所有区间都计算最后根据 x 的范围取舍
        //计算实部
        //[0,H] -2*x/H + 1 
        vresultRE1 = vaddq_f32(vdupq_n_f32(1.0),vdivq_f32_neon(vmulq_n_f32(vx,-2),vH));
        //[H,2π] (2(x-PI)-H)/(2PI-H)
        vresultRE2 = vdivq_f32_neon(vsubq_f32(vmulq_n_f32(vsubq_f32(vx,vPI),2),vH),v2PI_H);

        //计算虚部
        //[0,π-H/2] 2*x/(2PI-H)
        vresultIM1 = vdivq_f32_neon(vmulq_n_f32(vx,2),v2PI_H);
        //[π-H/2,π+H/2] -2*(x-PI)/H
        vresultIM2 = vdivq_f32_neon(vmulq_n_f32(vsubq_f32(vx,vPI),-2),vH);
        //[π+H/2,2π] 2(x-2PI)/(2PI-H)
        vresultIM3 = vdivq_f32_neon(vmulq_n_f32(vsubq_f32(vx,v2PI),2),v2PI_H);

        //实部取舍
        uint32x4_t vcompare = vcltq_f32(vx,vH);
    
        printf("c %u\n",vgetq_lane_u32(vcompare,0));

        vresultRE = vbslq_f32(vcompare,vresultRE1,vresultRE2);
        vresultRE = vmulq_n_f32(vresultRE,magn);
        //虚部取舍
        float32x4_t v2_H = vdivq_f32_neon(vH,vdupq_n_f32(2));
        // x< π-H/2 取 vresultIM1,否则取 vresultIM2
        vcompare = vcltq_f32(vx,vsubq_f32(vPI,v2_H));
        vresultIM_1 = vbslq_f32(vcompare,vresultIM1,vresultIM2);
        // x> π+H/2 取 vresultIM3,否则取 vresultIM
        vcompare = vcgtq_f32(vx,vaddq_f32(vPI,v2_H));
        vresultIM = vbslq_f32(vcompare,vresultIM3,vresultIM_1);
        vresultIM = vmulq_n_f32(vresultIM,magn);

        vresult.val[0] = vresultRE;
        vresult.val[1] = vresultIM;
        printf("re %f\n",vgetq_lane_f32(vresultRE,0));
        printf("im %f\n",vgetq_lane_f32(vresultIM,0));
        vst2q_f32((float32_t*)(pDst+i),vresult);
    }
    Ipp32f H = IPP_PI+asym;
    for(int i=len/4*4;i<len;i++)
    {
        Ipp32f x = IPP_2PI*rFreq*i + *pPhase;
        //x 需要转换到 [0,2π]
        Ipp32s k = x/IPP_2PI;
        x -= k*IPP_2PI;
        Ipp32fc result;
        // printf("x %f\n",x);
        //计算实部
        if(x>=0 && (H-x)>=EPSILON)  
             result.re = -2*x/H + 1.0;  // -2*x/H + 1 
        else    
            result.re = (2*(x-IPP_PI)-H)/(IPP_2PI-H);// (2(x-PI)-H)/(2PI-H)
        //计算虚部
        if(x>=0.0 && ((IPP_2PI-H)/2)-x >= EPSILON)  
            result.im = 2*x/(IPP_2PI-H);    // 2*x/(2PI-H)
        else if(x-(IPP_2PI-H)/2 >=EPSILON && (IPP_2PI+H)/2-x>=EPSILON)
            result.im = -2*(x-IPP_PI)/H;    //-2*(x-PI)/H
        else
            result.im = 2*(x-IPP_2PI)/(IPP_2PI-H);  //2(x-2PI)/(2PI-H)
        result.re *= magn;
        result.im *= magn;
        *(pDst+i) = result;
    }

    return ippStsNoErr;
}

3.IPP库测试

官网安装教程：https://www.intel.cn/content/www/cn/zh/developer/tools/oneapi/hpc-toolkit-download.html?packages=hpc-toolkit&hpc-toolkit-os=linux&hpc-toolkit-lin=offline

运行前设置环境变量

source /opt/intel/oneapi/setvars.sh

编译命令

icx test.c -o test -lippvm -lippcore -lipps

4.测试数据

函数 串行单精度 FIELD3 串行双精度 FIELD5 向量化单精度 FIELD7 intel IPP FIELD9
ippsExp_32f_A24 1.976372 0 1.976372 0
0.004269 0.00427 0.004269 0.004271
2527.535645 2527.535645 2527.535645 2527.535645
145801.3281 145801.3438 145801.3281 145801.3438
140399216 140399184 140399184 140399184
ippsCos_32f_A24 0.357278 0.357278 0.357279 0.357278
0.676947 0.67695 0.676952 0.67695
0.01898 0.01898 0.018981 0.01898
0.77985 0.77985 0.779854 0.77985
0.995996 0.995993 0.995993 0.995993
ippsSin_32f_A24 -0.933998 -0.933998 -0.933998 -0.933998
0.73603 0.736029 0.73603 0.736029
0.99982 0.99982 0.99982 0.99982
-0.625968 -0.625967 -0.625965 -0.625966
-0.089438 -0.089436 -0.089436 -0.089436
ippsTone_32f 1.872611 1.872609 1.872612 1.872609
2.394652 2.394654 2.394655 2.394655
-1.740218 -1.74022 -1.74022 -1.74022
-2.490863 -2.490862 -2.490861 -2.490863
1.602509 1.602513 1.602509 1.602513
ippsTriangle_32fc 实部 虚部 实部 虚部 实部 虚部 实部 虚部
-0.218555 -2.920779 -0.218555 -2.920779 -0.218555 -2.920779
2.158905 2.320416 2.158905 2.320415 2.158905 2.320415
0.15116 1.238589 0.15116 1.238589 0.15116 1.238588
-1.856586 -0.769157 -1.856585 -0.769157 -1.856585 -0.769157
-0.615485 -2.776901 -0.615485 -2.776901 -0.615484 -2.776901

5. 结果分析

(1) 用单精度计算单的结果可能误差较大，无法保证单精度十进制6位数的准确，建议使用双精度计算转换成单精度结果；

(2) NEON 提供有 vdiv 除法接口，比使用倒数的版本精度更高。

标签：20,Intel,示例,IPP,vx,f32,result,2PI,vresultIM
From： https://www.cnblogs.com/anluo8/p/18613007

相关文章

  • java面试问题(2024.12.17)
    记录java岗面试问题对java的了解Java是一门面向对象的编程语言，吸收了C++语言中大量的优点，但又抛弃了C++中容易出错的地方，如垃圾回收。Java又是一门平台无关的编程语言，通过java虚拟机（jvm)可以实现一次编译，处处运行。对jvm的了解Java虚拟机，是Java实现跨平台的关键所......
  • 20241217-封装、继承、多态
    1.封装目的在于保护数据安全，隐藏细节。1.1属性的封装//本文属性和方法都定义为静态的，也可不设为静态，由创建对象来调用和访问。publicclassTestCat{ publicstaticvoidmain(String[]args) { Cat.setWeight(2.3f); }}classCat{ privatestaticfloatweigh......
  • 20241213-局部变量和全局变量的思考
    for循环或while循环、方法或方法参数列表里定义的局部变量，在其内的代码块执行完毕后就被销毁了，不能再用了。1.A方法的局部变量a作为B方法的传入参数，在B方法内对该传入参数的运算不会对A方法的局部变量a产生影响。见下代码：publicclassArrayReference{ publicstaticvoid......
  • 2024《小学教材全解精析》1-6年级无偿分享
    前言我相信，从博客园发展到如今，很多人跟我一样，30出头，孩子都开始上小学了。自然而然，都离不开教育，也是最烦恼的！！！教育，是点亮孩子未来的明灯，而优质的学习资料则是这盏明灯中不可或缺的灯芯。在孩子的小学学习旅程中，每一个知识点都如同珍贵的宝石，需要精心雕琢与呵护，才能绽放出璀璨的光......
  • # easygui中所有函数用法示例
    #easygui中所有函数用法示例'''注意事项文件对话框：filesavebox和fileopenbox示例中的filetypes参数可以指定文件类型过滤器，例如["*.txt"]只显示文本文件。运行这个脚本，你将看到easygui提供的各种对话框，并可以测试它们的功能。'''importeasygui#1、msgbox:显......
  • CMake学习2024.12.9问AI的问题的记录
    vim如何将命令行窗口放入一个子页面中在Vim中将命令行窗口放入一个子页面（分割窗口）中，可以使用内置的终端功能。以下是具体步骤：##前提条件确保你使用的Vim版本支持终端功能。Vim8.0及以上版本和Neovim都支持内置终端。如果你使用的是较旧的Vim版本，建议升级或使用Ne......
  • 英国二手市场火爆！英国人今年圣诞将花20.5亿英镑买二手礼物
    随着全球可持续消费理念的普及，英国二手市场迎来了前所未有的增长，尤其在2024年的圣诞季，这一趋势更是表现得淋漓尽致。据数据显示，英国消费者今年圣诞节在二手礼物上的支出预计将达到20.5亿英镑，这一现象背后的原因和市场机遇值得深挖。**独角兽翻译器**一、英国二手市场为何火爆......
  • 2025年网络安全零基础自学全攻略：避开弯路，快速上手！
     自学网络安全是一项充满挑战的任务，但只要遵循合适的学习路径，能够有效避免走弯路，逐步建立知识体系，最终可以在该领域取得成就。本文将为你提供2025年最新的网络安全自学攻略，帮助你高效地规划学习路线，掌握网络安全的核心知识，避免学习中的误区。目录自学网络安全的......
  • webstorm2020破解激活教程
    1、下载WebStorm安装包下载webstorm2021以前的安装包，目前该安装包可能不太好找，笔者在此给大家提供下载方式（包含破解包）：链接:https://pan.baidu.com/s/1W3U94zHqWwAzw5FneF7lsg?pwd=6m4c提取码:6m4c2、安装WebStorm破解版下载好软件之后，直接进行安装，先选择30天免费模式......
  • 个人回忆录---2024软件工程历险记
    这个作业属于哪个课程https://edu.cnblogs.com/campus/fzu/SE2024这个作业要求在哪里https://edu.cnblogs.com/campus/fzu/SE2024/homework/13315这个作业的目标回顾这一学期所完成的软工任务，总结这一学期的收获学号102202108王露洁前言：我对软件工程课程的......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99