Ubuntu SSE指令集编程实例---复数乘法与共轭乘法

最新推荐文章于 2023-08-30 23:48:05 发布

光光学长

最新推荐文章于 2023-08-30 23:48:05 发布

阅读量4.4k

点赞数 8

分类专栏： Linux编程文章标签： SSE指令集并行计算 SIMD 编译

本文链接：https://blog.csdn.net/jxwxg/article/details/53091376

版权

Linux编程专栏收录该内容

2 篇文章 0 订阅

订阅专栏

Ubuntu SSE指令集编程实例—复数乘法与共轭乘法

如有相关问题，欢迎随时讨论交流 jxwxg@foxmail.com

1. SSE指令集介绍

　　1966 年，Micheal Flynn 提出根据指令和数据流的概念对计算机系统结构进行分类。他将计算机划分为四种基本类型：SISD，SIMD，MISD，MIMD。

SISD，即单指令单数据流，传统的冯诺依曼体系结构的计算机都属于这种类型，其体系结构如下图所示。
SISD

MISD，即多指令流单数据流，这种体系结构对单一数据进行多条不同指令的处理，比较适用于分类问题，其体系结构如下图所示。

MISD

SIMD，即单指令流多数据流，这种体系结构对不同的多组数据采用相同指令来处理，比较适合并行算法的实现，在语音、图像等信号的处理中应用广泛，其体系结构如下图所示。

SISD

MIMD，即多指令流多数据流，这种体系结构拥有若干数据流和若干指令，可以异步地同时执行多种操作，功能最为强大，也最难设计和实现，其系统结构如下图所示。

MIMD

　　SSE（Streaming SIMD Extensions）指令是Intel 公司发布的SIMD的扩展指令，它可以方便地实现并行计算，其优势包括：
　　　1. 更高分辨率的图像浏览和处理功能;
　　　2. 语音识别占用更少CPU 资源；
　　　3. 具备更高精度和更快响应速度。

　　SSE 是MMX（MultiMedia eXtensions）的超集，MMX 支持8 个64bit 的MMX 寄存器，支持在打包的字，字节，双字整数上的SIMD 操作。SSE 指令集的四个版本概要如下表所示。

SSE Version

2. SSE指令集原理

　　计算机CPU 的结构如下图所示，CPU 首先从指令计数器寄存器中查看指令地址，然后从相应的指令地址中取出相应的指令存入指令寄存器中。接着CPU 根据操作数地址从内存中取出相应的操作数放入寄存器组中，最后运算器根据操作指令对操作数进行运算。SIMD 就是从内存中取一条指令放在指令寄存器中，再从内存中取一组数据放在SIMD 寄存器中，然后用这条指令并行处理这一组数据。
CPU 结构

　　SSE2 以后就能一次处理128 比特数据，这128 比特数据可能的存储方式是16*8bit，8*16bit，4*32bit，2*64bit，1*128bit。如下图所示。
SSE2

　　多数处理器都提供了SIMD指令集,并增加了专门用于执行SIMD指令的SIMD寄存器。SIMD 寄存器的位宽多为32 位、64 位或者128 位,这个位宽被称为“字”。这些SIMD寄存器可以划分为若干个较小的相互独立的单元,称为“亚字”。SIMD 指令对SIMD 寄存器中的所有亚字同时执行相同的操作。下图给出了一个典型的SIMD 指令执行模型: 两组数据各包括四个相互独立的单元,每对单元执行相同的操作,可同时得到四个运算结果。
SSE Operation

　　例如将4 组int 型（32 比特）数据分别执行加运算，CPU 首先从内存中取一条指令（加指令）放入指令寄存器。然后从内存中取出4 组int 型数据放入SIMD 寄存器中，再通过运算器并行的把相对应寄存器中的数据相加，计算结果存入寄存器中。这实现了一条指令并行处理4 组int 型数据。

3. 常用SSE指令集

　　SSE 指令集包含整型算术指令，整型移位指令，整型逻辑指令，整型比较指令，整型混杂指令，整型读取存储指令，浮点型算术指令，浮点型逻辑指令，浮点型比较指令，浮点型读取存储指令等。此外，计算机处理整数速度比处理浮点数快。

　　查看SSE 指令时，可以参考以下规律：
　　　1. SSE 指令名字包含三段。以_mm 开头，中间是所做运算的简写，最后指明操作数的单位。
例如_mm_add_epi16(_m128i a,_m128i b)表示将a 和b 中对应位置的16bit 有符号整数分别相加，其中 addadd 表示加操作， epi16epi16 表示操作数的单位是 16bit 有符号的整数。

　　　2. 操作数的单位一般有 8bit ，16bit ，32bit ，64bit 64bit 四种，对应四种相应的的 SSE 指令，如_mm_add_epi8 ，_mm_add_epi16，_mm_add_epi32 ，_mm_add_epi64 。

　　　3. epi16epi16 表示 16bit 有符号整数， epu16epu16 表示 16bit 无符号整数。 lo 表示高， hi 表示高，例如高，例如_mm_unpacklo_epi16(_m128i a,_m128i b) 表示将寄存器 a和寄存器和寄存器 b的低 64bit 数以 16bit 为单位交织在一块。

　　一些常用的SSE指令如下表所示。

指令类型	指令	功能
整型算术指令	_mm_add_epi16(_m128i a,_m128i b)	返回一个_m128i的寄存器，将a和b中对应位置的16bit有符号整数分别相加。
	_mm_adds_epi16(_m128i a,_m128i b)	返回一个_m128i的寄存器，将a和b中对应位置的16bit有符号整数分别相加，当计算结果溢出时将其置为边界值
	_mm_madd_epi16(_m128i a,_m128i b)	返回一个_m128i的寄存器，它含有4个有符号32bit的整数，分别满足： r0 := (a0 * b0) + (a1 * b1), r1 := (a2 * b2) + (a3 * b3), r2 := (a4 * b4) + (a5 * b5), r3 := (a6 * b6) + (a7 * b7)
	_mm_mulhi_epi16(_m128i a,_m128i b)	返回一个_128i的寄存器，它含8个有符号16bit的整数，分别为a和b对应位置的16bit有符号整数相乘结果的高16bit数据，即,r0 := (a0 * b0)[31:16],r1 := (a1 * b1)[31:16],…,r7 := (a7 * b7)[31:16]
整型移位指令	_mm_slli_epi16(_m128i a, int count)	返回一个_m128i的寄存器，将寄存器a中的8个16bit整数按照count进行相同的逻辑左移
整型移位指令	_mm_sll_epi16(_m128i a,_m128i count)	返回一个_m128i的寄存器，将寄存器a中的8个16bit整数按照count寄存器中对应位置的整数进行逻辑左移
整型逻辑指令	_mm_and_si128(_m128i a,_m128ib)	返回一个_m128i的寄存器，将寄存器a和寄存器b的对应位进行按位与运算。
整型逻辑指令	_mm_xor_si128(_m128i a,_m128ib)	返回一个_m128i的寄存器，将寄存器a和寄存器b的对应位进行按位异或运算。
整型混杂指令	_mm_shuffle_epi32(_m128i a,int imm)	返回一个_m128i的寄存器，它是将a中128i数据以32bit为单位重新排列得到的，imm为一个四元组，表示重新排列的顺序。
	_mm_unpacklo_epi16(_m128i a,_m128i b)	返回一个_m128i的寄存器，它将寄存器a和寄存器b的低64bit数以16bit为单位交织在一块。
	_mm_extract_epi16(_m128i a,int imm)	返回一个16bit整数，根据imm从a中8个16bit数中选取对应编号的数。
整型读取存储指令	_mm_load_si128(_m128i *p)	返回一个_m128i的寄存器，它将p指向的数据读到指定寄存器中，实际使用时，p一般是通过类型转换得到的。

4. SSE实现复数乘法

数学原理

A = a + bj
B = c + dj
A * B = (a + bj)*(c ＋ dj) = (ac-bd) + (bc+ad)j

源码

\#include <stdio.h>
\#include <stdlib.h>
\#include "emmintrin.h"
\#include "tmmintrin.h"
\#include <stdint.h>
int16_t conjugate[8] __attribute__((aligned(16))) = {1,-1,1,-1,1,-1,1,-1};
void two_int_complex_multiply(int16_t *inputa, int16_t *inputb, int length, int16_t *output, int16_t shift)
{
    /*
     *    inputa: RE(16 bit),IM(16 bit); RE(16 bit),IM(16 bit)
     *    inputb: RE(16 bit),IM(16 bit); RE(16 bit),IM(16 bit)
     *    output: RE(16 bit),IM(16 bit); RE(16 bit),IM(16 bit)
     */
    __m128i *a128 = (__m128i *) inputa;
    __m128i *b128 = (__m128i *) inputb;
    __m128i *o128 = (__m128i *) output;
    __m128i mmtmp_re, mmtmp_im, mmtmp_lo, mmtmp_hi;
    int index;
    for(index = 0; index < (length >> 2); index++)
    {
        mmtmp_re = _mm_sign_epi16(*b128, *(__m128i*)&conjugate[0]);
        mmtmp_re = _mm_madd_epi16(*a128, mmtmp_re);
        mmtmp_im = _mm_shufflelo_epi16(*b128, _MM_SHUFFLE(2,3,0,1));
        mmtmp_im = _mm_shufflehi_epi16(mmtmp_im, _MM_SHUFFLE(2,3,0,1));
        mmtmp_im = _mm_madd_epi16(*a128, mmtmp_im);
        mmtmp_re = _mm_srai_epi32(mmtmp_re, shift);
        mmtmp_im = _mm_srai_epi32(mmtmp_im, shift);
        mmtmp_lo = _mm_unpacklo_epi32(mmtmp_re,mmtmp_im);
        mmtmp_hi = _mm_unpackhi_epi32(mmtmp_re,mmtmp_im);
        o128[0]  = _mm_packs_epi32(mmtmp_lo,mmtmp_hi);
        a128 ++;
        b128 ++;
        o128 ++;
    }
}
int main(int argc, char *argv[])
{
    int shift = 9;
    int data_length = 16;
    int *inputa = malloc( data_length * sizeof(int));
    int *inputb = malloc( data_length * sizeof(int));
    int *output = malloc( data_length * sizeof(int));
    int index;
    for(index = 0; index < data_length; index++)
    {
        ((int16_t *)&inputa[index])[0] = rand() % 512;
        ((int16_t *)&inputa[index])[1] = rand() % 512;
        ((int16_t *)&inputb[index])[0] = rand() % 512;
        ((int16_t *)&inputb[index])[1] = rand() % 512;
    }
    puts("Input data A:");
    for(index = 0; index < data_length; index++)
    {
        if(index % 10 == 9)
            puts("");
        printf("(%d %d) ",((int16_t *)&inputa[index])[0],((int16_t *)&inputa[index])[1]);
    }
    puts("\n");
    puts("Input data B:");
    for(index = 0; index < data_length; index++)
    {
        if(index % 10 == 9)
            puts("");
        printf("(%d %d) ",((int16_t *)&inputb[index])[0],((int16_t *)&inputb[index])[1]);
    }
    puts("\n");
    two_int_complex_multiply((int16_t *)inputa, (int16_t *)inputb, data_length, (int16_t *)output, shift);
    puts("output data:");
    for(index = 0; index < data_length; index++)
    {
        if(index % 10 == 9)
            puts("");
        printf("(%d %d) ",((int16_t *)&output[index])[0],((int16_t *)&output[index])[1]);
    }
    puts("");
    return 0;
}

编译
gcc two_int_complex_multiply.c -o two_int_complex_multiply -msse4.1

5. SSE实现复数共轭乘法

数学原理

A = a + bj
B = c + dj
A * conj(B) = (a + bj)*(c -dj) = (ac+bd) + (bc-ad)j

源码

\#include <stdio.h>
\#include <stdlib.h>
\#include "emmintrin.h"
\#include "tmmintrin.h"
\#include <stdint.h>
int16_t conjugate[8] __attribute__((aligned(16))) = {1,-1,1,-1,1,-1,1,-1};
void two_int_complex_conj_multiply(int16_t *inputa, int16_t *inputb, int length, int16_t *output, int16_t shift)
{
    /*
     *inputa: RE(16 bit),IM(16 bit); RE(16 bit),IM(16 bit)
     *inputb: RE(16 bit),IM(16 bit); RE(16 bit),IM(16 bit)
     *output: RE(16 bit),IM(16 bit); RE(16 bit),IM(16 bit)
     */
    __m128i *a128 = (__m128i *) inputa;
    __m128i *b128 = (__m128i *) inputb;
    __m128i *o128 = (__m128i *) output;
    __m128i mmtmp_re, mmtmp_im, mmtmp_lo, mmtmp_hi;
    int index;
    for(index = 0; index < (length >> 2); index++)
    {
        mmtmp_re = _mm_madd_epi16(*a128, *b128);
        *b128    = _mm_sign_epi16(*b128, *(__m128i*)&conjugate[0]);
        mmtmp_im = _mm_shufflelo_epi16(*b128, _MM_SHUFFLE(2,3,0,1));
        mmtmp_im = _mm_shufflehi_epi16(mmtmp_im, _MM_SHUFFLE(2,3,0,1));
        mmtmp_im = _mm_madd_epi16(*a128, mmtmp_im);
        mmtmp_re = _mm_srai_epi32(mmtmp_re, shift);
        mmtmp_im = _mm_srai_epi32(mmtmp_im, shift);
        mmtmp_lo = _mm_unpacklo_epi32(mmtmp_re,mmtmp_im);
        mmtmp_hi = _mm_unpackhi_epi32(mmtmp_re,mmtmp_im);
        o128[0]  = _mm_packs_epi32(mmtmp_lo,mmtmp_hi);
        a128 ++;
        b128 ++;
        o128 ++;
    }
}
int main(int argc, char *argv[])
{
    int shift = 9;
    int data_length = 16;
    int *inputa = malloc( data_length * sizeof(int));
    int *inputb = malloc( data_length * sizeof(int));
    int *output = malloc( data_length * sizeof(int));
    int index;
    for(index = 0; index < data_length; index++)
    {
        ((int16_t *)&inputa[index])[0] = rand() % 512;
        ((int16_t *)&inputa[index])[1] = rand() % 512;
        ((int16_t *)&inputb[index])[0] = rand() % 512;
        ((int16_t *)&inputb[index])[1] = rand() % 512;
    }
    puts("Input data A:");
    for(index = 0; index < data_length; index++)
    {
        if(index % 10 == 9)
            puts("");
        printf("(%d %d) ",((int16_t *)&inputa[index])[0],((int16_t *)&inputa[index])[1]);
    }
    puts("\n");
    puts("Input data B:");
    for(index = 0; index < data_length; index++)
    {
        if(index % 10 == 9)
            puts("");
        printf("(%d %d) ",((int16_t *)&inputb[index])[0],((int16_t *)&inputb[index])[1]);
    }
    puts("\n");
    two_int_complex_conj_multiply((int16_t *)inputa, (int16_t *)inputb, data_length, (int16_t *)output, shift);
    puts("output data:");
    for(index = 0; index < data_length; index++)
    {
        if(index % 10 == 9)
            puts("");
        printf("(%d %d) ",((int16_t *)&output[index])[0],((int16_t *)&output[index])[1]);
    }   

    puts("");
    return 0;
}

编译
gcc two_int_complex_conj_multiply.c -o two_int_complex_conj_multiply -msse4.1

光光学长

关注

8
点赞
踩
15

收藏

觉得还不错? 一键收藏
0
评论
Ubuntu SSE指令集编程实例---复数乘法与共轭乘法

Ubuntu SSE库编程实例---复数乘法与共轭乘法1. SSE指令集介绍2. SSE指令集原理3. 常用SSE指令集4. SSE实现复数乘法5. SSE实现复数共轭乘法
复制链接

扫一扫