基于鲲鹏数学库的高性能数学计算加速方法研究与实践

First Post:

Last Update:

Word Count:
10.1k

Read Time:
48 min

基于鲲鹏数学库的高性能数学计算加速方法研究与实践

1 实验目的

1. 验证鲲鹏数学库的性能优化效果

  • 核心目标:通过一系列数学计算实验,系统验证基于鲲鹏架构的 Kunpeng Math Library (KML) 相较传统实现方式的性能提升程度。
  • 细化目标
    • 验证 KML 在不同数学计算任务中的加速效果,包括基础数学运算(如三角函数计算)、矩阵计算(如矩阵乘法)和高级线性代数问题(如特征值分解)。
    • 比较 KML 在单线程和多线程模式下的性能表现,探索其并行能力。

2. 探索鲲鹏架构的硬件适配性

  • 硬件架构支持验证
    • 检验 KML 对鲲鹏处理器的硬件特性(如 SIMD 指令集、多核并行)和内存架构(如 L1/L2 缓存优化)的适配程度。
    • 探索基于 ARM 架构优化的数学库在高性能计算(HPC)场景中的实际表现。
  • 平台优势探索
    • 评估 KML 在典型计算任务(如矩阵运算、快速傅里叶变换)中对其他平台(如 x86 架构)的相对性能优势,为基于鲲鹏的数学计算方案提供参考。

3. 优化实际科学计算任务

  • 面向工程计算和科学研究
    • 通过 KML 的高性能计算加速,解决实际工程和科研中的计算瓶颈问题。例如,优化数据分析、机器学习训练中的矩阵运算,加速信号处理中的傅里叶变换等。
  • 扩展应用场景
    • 在本实验中探索 KML 在不同任务中的适配性,为其在更多计算场景(如气象模拟、分子动力学模拟等)提供基础。

4. 提升鲲鹏生态的实践经验

  • 用户指导与案例积累
    • 总结使用 KML 的安装、配置、调试方法,以及在各类任务中的具体实现步骤,为开发者提供实践指南。
    • 为鲲鹏生态的推广提供优化实践案例,进一步推动基于 ARM 架构的高性能计算技术在行业中的普及。
  • 性能优化经验总结
    • 通过对比实验,提炼出优化数学计算性能的经验,如向量化处理、多线程并行优化、内存管理优化等,为后续开发高效的数学计算程序提供理论和实践依据。

5. 评估性能与准确性的平衡

  • 性能测试
    • 在数学运算中,尤其是高级运算(如特征值分解、快速傅里叶变换),性能与准确性可能存在一定平衡点。实验目标之一是分析 KML 在保持高精度的同时实现性能提升的能力。
  • 结果一致性验证
    • 确保使用 KML 计算的结果在科学研究或工程应用中具有可接受的准确性,并验证其与传统实现方式的一致性。

6. 推动数学计算领域的技术进步

  • 实践新算法
    • 利用鲲鹏架构和 KML 优化数学运算的能力,探索传统算法在新型硬件架构下的优化潜力。
  • 高性能计算的技术转化
    • 将 KML 的优势推广到实际的工业和学术计算任务中,为数学计算领域的技术进步提供支持。

2 实验设备

  • 华为鲲鹏云服务器;

  • 具备网络连接的个人电脑。

3 实验原理

1. 数学库优化的核心思想

鲲鹏数学库(KML)的设计理念是通过充分利用硬件特性提升数学运算的性能。以下是其核心优化思想:

  • 指令级优化
    • KML 借助 ARMv8 架构的 SIMD(Single Instruction Multiple Data)指令集,能够在单次指令操作中对多个数据元素进行并行处理,从而显著提升运算速度。
    • 使用如 Neon 指令等优化数学函数(如三角函数、对数函数等)的批量处理能力。
  • 多线程并行
    • 通过 pthread 和 OpenMP 等多线程技术,KML 可高效分解计算任务,充分利用多核处理器的计算资源。
  • 内存访问优化
    • 避免频繁的内存读取与写入,通过对内存对齐、缓存友好型算法等优化策略,减少数据传输瓶颈。

2. 硬件架构与性能优化的关系

鲲鹏服务器基于 ARM 架构,具备以下硬件特性,KML 在设计时针对这些特性进行了针对性优化:

  • 高密度核心: 鲲鹏处理器具有高核数的特点,非常适合高并发计算任务,KML 可在矩阵运算、FFT 中有效分摊任务到多个核心,降低单核压力。
  • 宽向量寄存器: ARM 的 SIMD 技术利用宽向量寄存器,使得多数据并行计算成为可能。例如,KML_SVML 的优化实现可一次性处理多个向量元素。
  • 内存子系统优化: 鲲鹏架构的 L1 和 L2 缓存对矩阵计算中常见的密集存储访问模式进行了专门优化。KML 的矩阵运算函数充分利用了这一特性,通过块状处理减少了内存访问延迟。

3. 数学计算的多层次优化方法

KML 的优化覆盖了从低级指令到高级数学库的多个层次:

  • 基础数学函数优化: 优化基本数学函数(如三角函数、指数函数)的实现,通过硬件寄存器和流水线指令集实现批量处理。
  • BLAS 和 LAPACK 优化
    • BLAS:优化矩阵-向量运算、矩阵-矩阵运算等基础线性代数运算,广泛用于科学计算中。
    • LAPACK:在更高级别的线性代数操作(如特征值问题、奇异值分解)中,KML 将关键计算步骤分解为高效的 BLAS 调用,并结合硬件特性进行细粒度优化。
  • 快速傅里叶变换(FFT): FFT 算法中涉及复杂的递归和循环结构,KML 通过流水线技术和分块计算优化了时间复杂度,同时减少了内存访问次数。

4. 实验的对比设计与性能评估策略

为了全面评估 KML 的性能,实验采用了“传统实现”和“优化实现”的对比方式:

  • 传统实现:直接使用标准数学库(如 math.h)或手动实现算法,模拟常规的计算方式。
  • 优化实现:基于 KML 提供的库函数完成相同计算任务。
  • 性能评估
    • 时间复杂度:记录不同方法的执行时间,量化性能提升。
    • 准确性:验证 KML 的计算结果是否与传统实现一致,确保在优化性能的同时保证结果的可靠性。
    • 资源利用率:观察 CPU、内存的使用情况,分析 KML 如何在硬件资源利用上占据优势。

5. 典型优化场景

  • 矩阵运算: 矩阵计算(如矩阵-向量乘法、矩阵分解)是科学计算中的核心任务。KML 针对矩阵密集型计算进行了特殊优化,减少了循环嵌套带来的性能瓶颈。
  • 高维向量运算: 例如三角函数计算,普通实现逐元素计算效率低下,而 KML 利用向量化技术和批量计算方法,可显著加速处理大规模数据。
  • 特征值问题: 对称矩阵特征值与特征向量计算涉及复杂的迭代算法,KML 在 LAPACK 的基础上进一步优化了矩阵操作的并行性。

6. KML 的接口设计与易用性

  • 动态库链接: 用户可通过简单的编译选项链接 KML 提供的动态库,从而快速完成性能优化。
  • 模块化功能: 不同的子模块(如 BLAS、VML、SPBLAS 等)满足了从基础运算到高级线性代数的多种需求,提供了良好的兼容性和可扩展性。

4 实验任务操作指导

4.1安装KML

下面介绍如何安装KML

首先使用远程登录工具,登录到鲲鹏 ECS 服务器上

下载 WinSCP 客户端并安装。

启动WinSCP,启动后界面如下:

看不见是正常的,别担心

填写说明:

  • 协议:选填 SFTP 或者 SCP 均可。

  • 主机名:云服务器的公网 IP。登录管理控制台即可查看对应云服务器的公网 IP。

  • 端口:默认 22。

  • 用户名:云服务器的用户名。

  • 使用“SSH密钥方式”登录云服务器时:

    • 如果是“CoreOS”的公共镜像,用户名为“core”。

    • 如果是“非CoreOS”的公共镜像,用户名为“root”。

  • 使用“密码方式”登录云服务器,公共镜像(包括CoreOS)的用户名为:root。

  • 密码:购买云服务器设置的密码或通过密钥方式转化后的密码。

  • 单击“登录”,进入 “WinSCP” 文件传输界面。

  • 登录成功之后,您可以选择左侧本地计算机的文件,拖拽到右侧的远程云服务器,完成文件上传到云服务器。

具体操作可以参考:

本地Windows主机使用WinSCP上传文件到Linux云服务器 https://support.huaweicloud.com/ecs_faq/zh-cn_topic_0166284971.html

然后到https://www.hikunpeng.com/developer/boostkit/library/detail?subtab=%E6%95%B0%E5%AD%A6%E5%BA%93获取KML软件包(GCC版本)

看不见是正常的,别担心

看不见是正常的,别担心

下载软件包后解压得到此文件:

看不见是正常的,别担心

再通过“本地Windows主机使用WinSCP上传文件到Linux云服务器https://support.huaweicloud.com/ecs_faq/zh-cn_topic_0166284971.html”中所述的远程登录将该文件上传到云服务器

接下来安装KML。

步骤1 登录云服务器,进入刚刚上传文件所在的目录,输入

1
rpm -ivh kml-xxxx.aarch64.rpm

安装软件包,其中命令中涉及的xxxx代表版本号,下图示中的版本号是2.4.0-1

看不见是正常的,别担心

步骤 2 安装后验证

执行source命令或重新登录终端让环境变量生效。

1
source /etc/profile

看不见是正常的,别担心

查看环境变量LD_LIBRARY_PATH是否包含KML的安装路径“/usr/local/kml/lib”。

1
env | grep LD_LIBRARY_PATH

如果变量包含安装路径,说明安装成功。

看不见是正常的,别担心

安装成功后在安装路径(默认路径是“/usr/local/kml”)下生成相应文件,其中,include文件夹包含子库的头文件,lib文件夹包含了数学库的动态库文件。

看不见是正常的,别担心

使用时,请在GCC编译选项中添加动态库所在路径,链接需要使用的动态库文件,添加编译选项后用ldd命令检查程序依赖库是否准确链接。

若需要使用KML_BLAS请添加如下代码,此处对官方给出的代码进行适当修改以正常使用:

1
2
3
4
单线程不加锁版本:-L /usr/local/kml/lib/kblas/nolocking -lkblas -I /usr/local/kml/include -pthread
单线程加锁版本:-L /usr/local/kml/lib/kblas/locking -lkblas -I /usr/local/kml/include -pthread
pthread实现多线程版本:-L /usr/local/kml/lib/kblas/pthread -lkblas -I /usr/local/kml/include -pthread
OpenMP实现多线程版本:-L /usr/local/kml/lib/kblas/omp -lkblas -I /usr/local/kml/include -pthread

若需要使用KML_VML请添加:

1
2
单线程版本:-L /usr/local/kml/lib/kvml/single -lkvml -L /usr/local/kml/lib -lkm -I /usr/local/kml/include -lm
多线程版本:-L /usr/local/kml/lib/kvml/multi -lkvml -L /usr/local/kml/lib -lkm -I /usr/local/kml/include -lm -fopenmp

若需要使用KML_SPBLAS请添加:

1
2
单线程版本:-L /usr/local/kml/lib/kspblas/single -lkspblas -I /usr/local/kml/include -pthread
多线程版本:-L /usr/local/kml/lib/kspblas/multi -lkspblas -I /usr/local/kml/include -pthread

若需要使用KML_FFT请添加:

1
2
单精度版本:-L /usr/local/kml/lib -lkfftf -I /usr/local/kml/include -pthread
双精度版本:-L /usr/local/kml/lib -lkfft -I /usr/local/kml/include -pthread

若需要使用KML_MATH请添加:

1
2
高性能版本:-L /usr/local/kml/lib -lkm -lm -I /usr/local/kml/include -pthread
高精度版本:-L /usr/local/kml/lib -lkm_l9 -lm -I /usr/local/kml/include -pthread

若需要使用KML_SVML请添加:

1
-L /usr/local/kml/lib -lksvml -lm -I /usr/local/kml/include -pthread

若需要使用KML_VSL请添加:

1
-L /usr/local/kml/lib -lkvsl -I /usr/local/kml/include -pthread -lm

若需要使用KML_LAPACK:

先生成完整的LAPACK,然后添加:

1
2
3
4
export KML_LAPACK_ROOT=/usr/local/kml/lib
export ADAPT_ROOT=/home/lapack_adapt
export KML_BLAS_ROOT=/usr/local/kml/lib/kblas/omp
gcc test.c -o test -fopenmp -I $KML_LAPACK_ROOT/include/kml-0.3.0 -L /usr/local/kml/lib -lklapack -L $ADAPT_ROOT -l:liblapack_adapt.a -L $KML_BLAS_ROOT -lkblas -lgfortran -lm -lkservice -I /usr/local/kml/include

若需要使用KML_IPL请添加:

1
-L /usr/local/kml/lib -lkipl -lklapack_full -L /usr/local/kml/lib/kblas/pthread -lkblas -lm -I /usr/local/kml/include -pthread

若需要使用KML_SCALAPACK :
先生成完整的SCALAPACK,然后添加:

1
2
3
4
5
6
7
8
9
# 动态
gcc test.c -o app -fopenmp -I /usr/local/kml/lib/kblas/omp/include/kml-0.3.0 -L /usr/local/kml/lib -l:libkscalapack.a -L /home/scalapack_adapt -l:libscalapack_adapt.a -L /usr/local/kml/lib/kblas/omp -l kblas -L /usr/local/kml/lib -l:libkservice.a -L /home/lapack_adapt -l:liblapack_adapt.a -lm -I /usr/local/kml/include

# 静态
export KML_LAPACK_ROOT=/usr/local/kml/lib
export ADAPT_ROOT=/home/scalapack_adapt
export KML_BLAS_ROOT=/usr/local/kml/lib/kblas/omp

gcc test.c -o app -fopenmp -I $KML_LAPACK_ROOT/include/kml-0.3.0 -L /usr/local/kml/lib -l:libkscalapack.a -L $ADAPT_ROOT -l:libscalapack_adapt.a -L $KML_BLAS_ROOT -L /home/lapack_adapt -l:liblapack_adapt.a -l:libkservice.a -lm -I /usr/local/kml/include

下面给出两个基础的测试程序用于测试是否已经成功安装

使用mkdir创建文件夹FFTTEST,使用cd FFTTEST进入,vim test.c创建测试文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include <stdio.h>
#include "km.h"

int main()
{
double pi = acos(-1);
// typical usage
double a = pi/6, b = 1.0, c = -3*pi/4, d = pi/3;
// print result
printf("sin(pi/6) = %.15f\n", sin(a));
printf("sin(1.0) = %.15f\n", sin(b));
printf("sin(-3*pi/4) = %.15f\n", sin(c));
/*
* sin(pi/6) = 0.500000000000000
* sin(1.0) = 0.841470984807897
* sin(-3*pi/4) = -0.707106781186548
* sin(pi/3) = 0.866025403784438
* */
return 0;
}

利用

1
gcc test.c -I /usr/local/kml/include -L /usr/local/kml/lib -lkm_l9 -lm

编译,生成a.out

看不见是正常的,别担心

输入

1
./a.out

运行,结果如下

看不见是正常的,别担心

接下来尝试用牛顿迭代法求解非线性方程的根

依次执行命令 mkdir NEWTON、cd NEWTON 创建并进入到 NEWTON 目录。

创建 sum.c 文件,编写内容如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
//牛顿迭代法求解非线性方程的根
#include <stdio.h>
#include "km.h"
double f(double x)
{
return x * x - 2; // 函数 f(x) = x^2 - 2
}
double df(double x)
{
return 2 * x; // 函数的导数 f'(x) = 2x
}
double newton(double initial_guess, double tolerance, int max_iterations)
{
double x = initial_guess;
int iteration = 0;

while (iteration < max_iterations) {
double fx = f(x);
if (fabs(fx) < tolerance) {
return x; // 找到根,返回当前值
}
double dfx = df(x);
if (dfx == 0) {
printf("Error: Derivative is zero. No solution found.\n");
return x; // 导数为零,无法继续
}
x = x - fx / dfx; // 牛顿迭代公式
iteration++;
}
printf("Max iterations reached. Last approximation: %f\n", x);
return x; // 返回最后的近似值
}
int main() {
double initial_guess = 1.0; // 初始猜测
double tolerance = 1e-7; // 容忍度
int max_iterations = 100; // 最大迭代次数
double root;
printf("With KML:\n");
root = newton(initial_guess, tolerance, max_iterations);
printf("Root found: %f\n", root);
return 0;
}

输入

1
2
gcc NEWTON_KML.c -I /usr/local/kml/include -L /usr/local/kml/lib -lkm -lm -o NEWTON_KML
./NEWTON_KML

运行之,结构如下:

看不见是正常的,别担心

至此说明kml安装成功

4.2 矩阵-向量乘加运算对比实验

该代码主要进行了如下工作:初始化规模为1000*300的矩阵A,长度为300的向量x,长度为1000的向量y1和y2。分别用两种方法求矩阵-向量乘加运算,即y=alpha*A*x+beta*y:

(1)按照矩阵-向量的乘加规则实现算法求解

(2)调用KML_BLAS提供的函数cblas_dgemv求解 分别用计时器记录两种方法消耗的时间,对比KML_BLAS与手动实现矩阵乘加的性能。

编写代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#include<stdio.h>
#include<time.h>
#include "kblas.h"
#define M 1000
#define N 300
int main()
{
double A[M*N]={0};
double x[N]={0};
double y1[M]={0};
double y2[M]={0};
for (int i=0;i<M*N;i++)
{
A[i]=i;
}
for (int i=0;i<N;i++)
{
x[i]=i;
}
for (int i=0;i<M;i++)
{
y1[i]=1;
y2[i]=1;
}
double alpha=1.2;
double beta=2.5;
struct timespec t1,t2; //定义初始与结束时间
printf("Without KML:\n");
clock_gettime(CLOCK_MONOTONIC,&t1); //计算开始时间。
for(int i=0;i<M;i++)
{
double tmp=0;
for(int j=0;j<N;j++)
{
tmp+=(A[i*N+j]*x[j]);
}
y1[i]=alpha*tmp+beta*y1[i];
}
clock_gettime(CLOCK_MONOTONIC,&t2); //计算结束时间。
printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec);
printf("With KML:\n");
clock_gettime(CLOCK_MONOTONIC,&t1); //计算开始时间。
cblas_dgemv(CblasColMajor,CblasNoTrans, M, N, alpha, A, M, x, 1, beta, y2, 1);
clock_gettime(CLOCK_MONOTONIC,&t2); //计算结束时间。
printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec);

return 0;

}

输入:

1
gcc kblas.c -o kblas -L /usr/local/kml/lib/kblas/nolocking -lkblas -I /usr/local/kml/include

运行之,得到运行结果如下:

看不见是正常的,别担心

可知KML_BLAS提供的函数cblas_dgemv对矩阵-向量乘加运算的优化效果十分显著

4.3 调用KML_VML库计算100000个数的sin值和普通循环对比体现优化

使用多线程版本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#include<stdio.h>
#include<time.h>
#include<math.h>
#include"kvml.h"
#define LEN 100000

int main()
{
   double src[LEN]={0};
   double dst1[LEN]={0};
   double dst2[LEN]={0};
   for(int i=0;i<LEN;i++)
       src[i]=i;
    struct timespec t1,t2;    //定义初始与结束时间

    printf("Without KML_VML ");
    clock_gettime(CLOCK_MONOTONIC,&t1);      //计算开始时间。
    for(int i=0;i<LEN;i++)
       dst1[i]=sin(src[i]);
    clock_gettime(CLOCK_MONOTONIC,&t2);    //计算结束时间。
    printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec);  //得出目标代码段的执行时间。

    printf("With KML_VML ");
    clock_gettime(CLOCK_MONOTONIC,&t1);      //计算开始时间。
    vdsin(LEN,src,dst2);
    clock_gettime(CLOCK_MONOTONIC,&t2);    //计算结束时间。
    printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec);  //得出目标代码段的执行时间。

return 0;

}

输入

1
2
gcc vsin.c -o vsin -L /usr/local/kml/lib/kvml/multi -lkvml -L /usr/local/kml/lib -lkm -I /usr/local/kml/include -lm -fopenmp
./vsin

运行结果如下所示:

看不见是正常的,别担心

可知KML_VML库对三角函数运算的优化效果十分显著

4.4 体现LAPACK库性能的对比实验

首先要生成完整的LAPACK用到的脚本

接下来的程序需要脚本才能正确运行

在/home目录下创建sh.sh文件,接着编写如下代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
set -eE

echo "LAPACK_SRC_DIR         ${LAPACK_SRC_DIR:-<undefined>}"
echo "LAPACK_TGZ             ${LAPACK_TGZ:=/home/lapack-3.12.0.tar.gz}"
echo "LIBKLAPACK_A           ${LIBKLAPACK_A:=/usr/local/kml/lib/libklapack.a}"
echo "LIBKSERVICE_A          ${LIBKSERVICE_A:=${LIBKLAPACK_A/klapack/kservice}}"
echo "ADAPT_DIR              ${ADAPT_DIR:=./lapack_adapt}"
echo "CMAKE_BUILD_TYPE       ${CMAKE_BUILD_TYPE:=Release}"
echo "LIBLAPACK_ADAPT_A      ${LIBLAPACK_ADAPT_A:=liblapack_adapt.a}"
echo "LIBKLAPACK_FULL_SO     ${LIBKLAPACK_FULL_SO:=libklapack_full.so}"
echo "CC                     ${CC:=gcc}"
echo "FC                     ${FC:=gfortran}"

mkdir -p ${ADAPT_DIR}ZQ
cd ${ADAPT_DIR}

# build netlib lapack
if [ ! -r "${LAPACK_SRC_DIR}/CMakeLists.txt" ]; then
    mkdir -p netlib
    ( cd netlib ; tar xzpf ${LAPACK_TGZ} )
    LAPACK_SRC_DIR=$(cd netlib/l* ; pwd)
fi

mkdir -p build
cmake_flags=(
    -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
    -DCMAKE_POSITION_INDEPENDENT_CODE=ON
    -DCMAKE_C_COMPILER=${CC}
    -DCMAKE_Fortran_COMPILER=${FC}
    -DCMAKE_RULE_MESSAGES=off
    -DBUILD_DEPRECATED=on
    -DBUILD_TESTING=off
)
( cd build ; cmake ${cmake_flags[*]} ${LAPACK_SRC_DIR} )
( cd build ; make -j )

cp build/lib/liblapack.a ${LIBLAPACK_ADAPT_A}

# get symbols defined both in klapack and netlib lapack
nm -g ${LIBLAPACK_ADAPT_A} | grep 'T ' | grep -oP '\K\w+(?=_$)' | sort | uniq > netlib.sym
nm -g ${LIBKLAPACK_A} | grep 'T ' | grep -oP '\K\w+(?=_$)' | sort | uniq > klapack.sym
comm -12 klapack.sym netlib.sym > comm.sym

# update symbols name of ${LIBLAPACK_ADAPT_A}
while read sym; do
    (
        if ! nm ${LIBLAPACK_ADAPT_A} | grep -qe " T ${sym}_\$"; then
            continue
        fi
        ar x ${LIBLAPACK_ADAPT_A} ${sym}.f.o
        mv ${sym}.f.o ${sym}_netlib.f.o

        objcopy --redefine-sym ${sym}_=${sym}_netlib_ ${sym}_netlib.f.o
    ) &
done < comm.sym
wait
ar d ${LIBLAPACK_ADAPT_A} $(sed -ne 's/$/.f.o/p' comm.sym)
ar d ${LIBLAPACK_ADAPT_A} xerbla.f.o
ar ru ${LIBLAPACK_ADAPT_A} *_netlib.f.o
rm *_netlib.f.o

在/home目录下输入

1
sh sh.sh

运行之,得到如下结果:

看不见是正常的,别担心

看不见是正常的,别担心

至此LAPACK用到的脚本已经生成完毕

接下来先展示如何使用一般c语言代码计算实对称矩阵的特征值与特征向量

编写eigenNoOpt.c文件,文件内容如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

float** Matrix_Jac_Eig(float **arrayint n, float *eig);
int Matrix_Free(float **tmp, int m, int n);
int main(void)
{
    int n;
    printf("请输入矩阵维度:\n");
    scanf("%d", &n);
    float **array = (float **)malloc(n * sizeof(float *));
    if (array == NULL)
    {
        printf("error :申请数组内存空间失败\n");
        return -1;
    }
    for (int i = 0; i < n; i++)
    {
        array[i] = (float *)malloc(n * sizeof(float));
        if (array[i] == NULL)
        {
            printf("error :申请数组子内存空间失败\n");
            return -1;
        }
    }
    printf("请输入矩阵元素:\n");
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            scanf("%f", &array[i][j]);
        }
    }
    float *eig = (float *)malloc(n * sizeof(float));
    
    struct timespec t1,t2;

    printf("Without KML:\n");
    clock_gettime(CLOCK_MONOTONIC,&t1);      //计算开始时间。
    float **Result = Matrix_Jac_Eig(array, n, eig);
    clock_gettime(CLOCK_MONOTONIC,&t2);    //计算结束时间。
    printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec);  

    printf("特征矩阵元素:\n");
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            printf("%f ", Result[i][j]);
        }
        printf("\n");
    }
    printf("特征根:\n");
    for (int i = 0; i < n; i++)
    {
        printf("%f \n", eig[i]);
    }
    Matrix_Free(Result, n, n);
    free(eig);
    eig = NULL;
    return 0;
}
float** Matrix_Jac_Eig(float **arrayint n, float *eig)
{
    int i, j, flag, k;
    flag = 0;
    k = 0;
    float sum = 0;
    float **temp_mat = (float **)malloc(n * sizeof(float *));
    for (i = 0; i < n; i++)
    {
        temp_mat[i] = (float *)malloc(n * sizeof(float));
    }
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            temp_mat[i][j] = array[i][j];
        }
    }
    //判断是否为对称矩阵
    for (i = 0; i < n; i++)
    {
        for (j = i; j < n; j++)
        {
            if (array[i][j] != array[j][i])
            {
                flag = 1;
                break;
            }
        }
    }
    if (flag == 1)
    {
        printf("error in Matrix_Eig: 输入并非是对称矩阵:\n");
        return NULL;
    }
    else
    {
        //开始执行算法
        int p, q;
        float thresh = 0.0000000001;
        float max = array[0][1];
        float tan_angle, sin_angle, cos_angle;
        float **result = (float **)malloc(n * sizeof(float *));
        if (result == NULL)
        {
            printf("error in Matrix_Eig:申请空间失败\n");
            return NULL;
        }
        float **result_temp = (float **)malloc(n * sizeof(float *));
        if (result_temp == NULL)
        {
            printf("error in Matrix_Eig:申请空间失败\n");
            return NULL;
        }
        float **rot = (float **)malloc(n * sizeof(float *));
        if (rot == NULL)
        {
            printf("error in Matrix_Eig:申请空间失败\n");
            return NULL;
        }
        float **mat = (float **)malloc(n * sizeof(float *));
        if (mat == NULL)
        {
            printf("error in Matrix_Eig:申请空间失败\n");
            return NULL;
        }
        for (i = 0; i < n; i++)
        {
            result[i] = (float *)malloc(n * sizeof(float));
            if (result[i] == NULL)
            {
                printf("error in Matrix_Eig:申请子空间失败\n");
                return NULL;
            }
            result_temp[i] = (float *)malloc(n * sizeof(float));
            if (result_temp[i] == NULL)
            {
                printf("error in Matrix_Eig:申请子空间失败\n");
                return NULL;
            }
            rot[i] = (float *)malloc(n * sizeof(float));
            if (rot[i] == NULL)
            {
                printf("error in Matrix_Eig:申请子空间失败\n");
                return NULL;
            }
            mat[i] = (float *)malloc(n * sizeof(float));
            if (mat[i] == NULL)
            {
                printf("error in Matrix_Eig:申请子空间失败\n");
                return NULL;
            }
        }
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                if (i == j)
                {
                    result[i][j] = 1;
                }
                else
                {
                    result[i][j] = 0;
                }
            }
        }
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                if (i == j)
                {
                    mat[i][j] = 1;
                }
                else
                {
                    mat[i][j] = 0;
                }
            }
        }
        max = array[0][1];
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                if (i == j)
                {
                    continue;
                }
                else
                {
                    if (fabs(array[i][j]) >= fabs(max))
                    {
                        max = array[i][j];
                        p = i;
                        q = j;
                    }
                    else
                    {
                        continue;
                    }
                }
            }
        }
        while (fabs(max) > thresh)
        {
            if (fabs(max) < thresh)
            {
                break;
            }
            tan_angle = -2 * array[p][q] / (array[q][q] - array[p][p]);
            sin_angle = sin(0.5*atan(tan_angle));
            cos_angle = cos(0.5*atan(tan_angle));
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {

                    if (i == j)
                    {
                        mat[i][j] = 1;
                    }
                    else
                    {
                        mat[i][j] = 0;
                    }
                }
            }
            mat[p][p] = cos_angle;
            mat[q][q] = cos_angle;
            mat[q][p] = sin_angle;
            mat[p][q] = -sin_angle;
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    rot[i][j] = array[i][j];
                }
            }
            for (j = 0; j < n; j++)
            {
                rot[p][j] = cos_angle*array[p][j] + sin_angle*array[q][j];
                rot[q][j] = -sin_angle*array[p][j] + cos_angle*array[q][j];
                rot[j][p] = cos_angle*array[j][p] + sin_angle*array[j][q];
                rot[j][q] = -sin_angle*array[j][p] + cos_angle*array[j][q];
            }
            rot[p][p] = array[p][p] * cos_angle*cos_angle +
                array[q][q] * sin_angle*sin_angle +
                2 * array[p][q] * cos_angle*sin_angle;
            rot[q][q] = array[q][q] * cos_angle*cos_angle +
                array[p][p] * sin_angle*sin_angle -
                2 * array[p][q] * cos_angle*sin_angle;
            rot[p][q] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +
                array[p][q] * (2 * cos_angle*cos_angle - 1);
            rot[q][p] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +
                array[p][q] * (2 * cos_angle*cos_angle - 1);
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    array[i][j] = rot[i][j];
                }
            }
            max = array[0][1];
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    if (i == j)
                    {
                        continue;
                    }
                    else
                    {
                        if (fabs(array[i][j]) >= fabs(max))
                        {
                            max = array[i][j];
                            p = i;
                            q = j;
                        }
                        else
                        {
                            continue;
                        }
                    }
                }
            }
            for (i = 0; i < n; i++)
            {
                eig[i] = array[i][i];
            }
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    sum = 0;
                    for (k = 0; k < n; k++)
                    {
                        sum = sum + result[i][k] * mat[k][j];
                    }
                    result_temp[i][j] = sum;
                }
            }
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    result[i][j] = result_temp[i][j];
                }
            }
        }
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                array[i][j] = temp_mat[i][j];
            }
        }
        Matrix_Free(result_temp, n, n);
        Matrix_Free(rot, n, n);
        Matrix_Free(mat, n, n);
        Matrix_Free(temp_mat, n, n);
        return result;
    }
}
int Matrix_Free(float **tmp, int m, int n)
{
    int i, j;
    if (tmp == NULL)
    {
        return(1);
    }
    for (i = 0; i < m; i++)
    {
        if (tmp[i] != NULL)
        {
            free(tmp[i]);
            tmp[i] = NULL;
        }
    }
    if (tmp != NULL)
    {
        free(tmp);
        tmp = NULL;
    }
    return(0);
}

输入如下代码对其进行编译并运行:

1
2
gcc eigenNoOpt.c -o eigenNoOpt -lm
./eigenNoOpt

得到运行结果如下:

看不见是正常的,别担心

对其进行改进,并作体现LAPACK库性能的对比实验如下:

将传统C语言算法的代码改为如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

float** Matrix_Jac_Eig(float **arrayint n, float *eig);
int Matrix_Free(float **tmp, int m, int n);
int main(void)
{
    int n;
    printf("请输入矩阵维度:\n");
    scanf("%d", &n);
    float **array = (float **)malloc(n * sizeof(float *));
    if (array == NULL)
    {
        printf("error :申请数组内存空间失败\n");
        return -1;
    }
    for (int i = 0; i < n; i++)
    {
        array[i] = (float *)malloc(n * sizeof(float));
        if (array[i] == NULL)
        {
            printf("error :申请数组子内存空间失败\n");
            return -1;
        }
    }
    printf("请输入矩阵元素:\n");
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            scanf("%f", &array[i][j]);
        }
    }
    float *eig = (float *)malloc(n * sizeof(float));
    
    struct timespec t1,t2;

    printf("Without KML:\n");
    clock_gettime(CLOCK_MONOTONIC,&t1);      //计算开始时间。
    float **Result = Matrix_Jac_Eig(array, n, eig);
    clock_gettime(CLOCK_MONOTONIC,&t2);    //计算结束时间。
    printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec);  
    Matrix_Free(Result, n, n);
    free(eig);
    eig = NULL;
    return 0;
}
float** Matrix_Jac_Eig(float **arrayint n, float *eig)
{
    int i, j, flag, k;
    flag = 0;
    k = 0;
    float sum = 0;
    float **temp_mat = (float **)malloc(n * sizeof(float *));
    for (i = 0; i < n; i++)
    {
        temp_mat[i] = (float *)malloc(n * sizeof(float));
    }
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            temp_mat[i][j] = array[i][j];
        }
    }
    //判断是否为对称矩阵
    for (i = 0; i < n; i++)
    {
        for (j = i; j < n; j++)
        {
            if (array[i][j] != array[j][i])
            {
                flag = 1;
                break;
            }
        }
    }
    if (flag == 1)
    {
        printf("error in Matrix_Eig: 输入并非是对称矩阵:\n");
        return NULL;
    }
    else
    {
        //开始执行算法
        int p, q;
        float thresh = 0.0000000001;
        float max = array[0][1];
        float tan_angle, sin_angle, cos_angle;
        float **result = (float **)malloc(n * sizeof(float *));
        if (result == NULL)
        {
            printf("error in Matrix_Eig:申请空间失败\n");
            return NULL;
        }
        float **result_temp = (float **)malloc(n * sizeof(float *));
        if (result_temp == NULL)
        {
            printf("error in Matrix_Eig:申请空间失败\n");
            return NULL;
        }
        float **rot = (float **)malloc(n * sizeof(float *));
        if (rot == NULL)
        {
            printf("error in Matrix_Eig:申请空间失败\n");
            return NULL;
        }
        float **mat = (float **)malloc(n * sizeof(float *));
        if (mat == NULL)
        {
            printf("error in Matrix_Eig:申请空间失败\n");
            return NULL;
        }
        for (i = 0; i < n; i++)
        {
            result[i] = (float *)malloc(n * sizeof(float));
            if (result[i] == NULL)
            {
                printf("error in Matrix_Eig:申请子空间失败\n");
                return NULL;
            }
            result_temp[i] = (float *)malloc(n * sizeof(float));
            if (result_temp[i] == NULL)
            {
                printf("error in Matrix_Eig:申请子空间失败\n");
                return NULL;
            }
            rot[i] = (float *)malloc(n * sizeof(float));
            if (rot[i] == NULL)
            {
                printf("error in Matrix_Eig:申请子空间失败\n");
                return NULL;
            }
            mat[i] = (float *)malloc(n * sizeof(float));
            if (mat[i] == NULL)
            {
                printf("error in Matrix_Eig:申请子空间失败\n");
                return NULL;
            }
        }
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                if (i == j)
                {
                    result[i][j] = 1;
                }
                else
                {
                    result[i][j] = 0;
                }
            }
        }
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                if (i == j)
                {
                    mat[i][j] = 1;
                }
                else
                {
                    mat[i][j] = 0;
                }
            }
        }
        max = array[0][1];
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                if (i == j)
                {
                    continue;
                }
                else
                {
                    if (fabs(array[i][j]) >= fabs(max))
                    {
                        max = array[i][j];
                        p = i;
                        q = j;
                    }
                    else
                    {
                        continue;
                    }
                }
            }
        }
        while (fabs(max) > thresh)
        {
            if (fabs(max) < thresh)
            {
                break;
            }
            tan_angle = -2 * array[p][q] / (array[q][q] - array[p][p]);
            sin_angle = sin(0.5*atan(tan_angle));
            cos_angle = cos(0.5*atan(tan_angle));
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {

                    if (i == j)
                    {
                        mat[i][j] = 1;
                    }
                    else
                    {
                        mat[i][j] = 0;
                    }
                }
            }
            mat[p][p] = cos_angle;
            mat[q][q] = cos_angle;
            mat[q][p] = sin_angle;
            mat[p][q] = -sin_angle;
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    rot[i][j] = array[i][j];
                }
            }
            for (j = 0; j < n; j++)
            {
                rot[p][j] = cos_angle*array[p][j] + sin_angle*array[q][j];
                rot[q][j] = -sin_angle*array[p][j] + cos_angle*array[q][j];
                rot[j][p] = cos_angle*array[j][p] + sin_angle*array[j][q];
                rot[j][q] = -sin_angle*array[j][p] + cos_angle*array[j][q];
            }
            rot[p][p] = array[p][p] * cos_angle*cos_angle +
                array[q][q] * sin_angle*sin_angle +
                2 * array[p][q] * cos_angle*sin_angle;
            rot[q][q] = array[q][q] * cos_angle*cos_angle +
                array[p][p] * sin_angle*sin_angle -
                2 * array[p][q] * cos_angle*sin_angle;
            rot[p][q] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +
                array[p][q] * (2 * cos_angle*cos_angle - 1);
            rot[q][p] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +
                array[p][q] * (2 * cos_angle*cos_angle - 1);
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    array[i][j] = rot[i][j];
                }
            }
            max = array[0][1];
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    if (i == j)
                    {
                        continue;
                    }
                    else
                    {
                        if (fabs(array[i][j]) >= fabs(max))
                        {
                            max = array[i][j];
                            p = i;
                            q = j;
                        }
                        else
                        {
                            continue;
                        }
                    }
                }
            }
            for (i = 0; i < n; i++)
            {
                eig[i] = array[i][i];
            }
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    sum = 0;
                    for (k = 0; k < n; k++)
                    {
                        sum = sum + result[i][k] * mat[k][j];
                    }
                    result_temp[i][j] = sum;
                }
            }
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    result[i][j] = result_temp[i][j];
                }
            }
        }
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                array[i][j] = temp_mat[i][j];
            }
        }
        Matrix_Free(result_temp, n, n);
        Matrix_Free(rot, n, n);
        Matrix_Free(mat, n, n);
        Matrix_Free(temp_mat, n, n);
        return result;
    }
}
int Matrix_Free(float **tmp, int m, int n)
{
    int i, j;
    if (tmp == NULL)
    {
        return(1);
    }
    for (i = 0; i < m; i++)
    {
        if (tmp[i] != NULL)
        {
            free(tmp[i]);
            tmp[i] = NULL;
        }
    }
    if (tmp != NULL)
    {
        free(tmp);
        tmp = NULL;
    }
    return(0);
}

将调用KML库的那段代码中改为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include<stdio.h>
#include <stdlib.h>
#include<time.h>
#include<stdlib.h>
#include "klapack.h"

int main()
{  
    char jobz = 'V'
    char uplo = 'L'
    int n = 10
    int lda = 10
    int info = 0
    double w[10]; 
    double *work = NULL
    double qwork; 
    int lwork = -1
    int *iwork = NULL
    int qiwork; 
    int liwork = -1

    double a[] =
     {  
1.232.753.104.565.927.018.409.886.344.59
2.751.562.493.314.475.906.557.108.225.78,  
3.102.491.885.396.237.683.794.126.447.85,  
4.563.315.398.944.582.536.837.461.563.21,  
5.924.476.234.589.768.905.123.982.636.34,  
7.015.907.682.538.906.754.801.943.552.88,  
8.406.553.796.835.124.801.649.205.904.75,  
9.887.104.127.463.981.949.207.442.386.71,  
6.348.226.441.562.633.555.902.380.762.13,  
4.595.787.853.216.342.884.756.712.139.32    }; 
    struct timespec t1,t2;    //定义初始与结束时间
    /* Query optimal work size */ 
    dsyevd_(&jobz, &uplo, &n, a, &lda, w, &qwork, &lwork, &qiwork, &liwork, &info); 
    if (info != 0
    { 
        return -1
    } 
    lwork = (int)qwork; 
    work = (double *)malloc(sizeof(double) * lwork); 
    liwork = (int)qiwork; 
    iwork = (int *)malloc(sizeof(int) * liwork); 
    /* Calculate eigenvalues and eigenvectors */ 
    printf("With KML:\n");
    clock_gettime(CLOCK_MONOTONIC,&t1);      //计算开始时间。
    dsyevd_(&jobz, &uplo, &n, a, &lda, w, work, &lwork, iwork, &liwork, &info); 
    clock_gettime(CLOCK_MONOTONIC,&t2);    //计算结束时间。  
    printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec);  

    free(work); 
    free(iwork); 

     
     return 0;
}

编译并运行,得到结果如下:

看不见是正常的,别担心

可知LAPACK库对运算的优化效果十分显著

4.5 使用KML库函数实现R2R变换优化

首先是使用一般C语言实现R2R变换

编写without.c文件,文件内容如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

void r2r_fft(double *in, double *out, int n) 
{
    int k, m, j;
    double wr, wi, arg, c, s;
    double *temp = (double*)malloc(n * sizeof(double));
    
    // 实数的 R2R FFT 变换实现
    for (k = 0; k < n; k++) {
        temp[k] = in[k];
    }

    // R2R FFT 变换
    for (m = 1; m <= n / 2; m *= 2) {
        wr = cos(M_PI / m);
        wi = sin(M_PI / m);
        for (k = 0; k < n; k += 2 * m) {
            for (j = 0; j < m; j++) {
                arg = j * (2 * M_PI / (2 * m));
                c = cos(arg);
                s = sin(arg);
                double u1 = temp[k + j];
                double u2 = temp[k + j + m];
                double v1 = temp[k + j + m + 1];
                out[k + j] = u1 + u2 * c - v1 * s;
                out[k + j + 1] = u1 * s + u2 * v1 * c;
            }
        }
    }

    free(temp);
}

int main() {
    int rank = 2
    int *n; 
    struct timespec t1t2; 

    n = (int*)malloc(sizeof(int) * rank); 
    n[0] = 128;  // 修改数据规模
    n[1] = 128;  // 修改数据规模

    double *in = (double*)malloc(sizeof(double) * n[0] * n[1]); 
    for (int i = 0; i < n[0] * n[1]; i++) { 
        in[i] = (double)(i % 10);  // 用一些基本的数值初始化输入数据
    } 

    double *out = (double*)malloc(sizeof(double) * n[0] * n[1]);

    clock_gettime(CLOCK_MONOTONIC, &t1);      // 计算开始时间
    r2r_fft(in, out, n[0] * n[1]);            // 执行 R2R FFT
    clock_gettime(CLOCK_MONOTONIC, &t2);      // 计算结束时间

    printf("Without KML:\n");
    printf("Time: %11u ns\n", t2.tv_nsec - t1.tv_nsec);  

    free(n);
    free(in); 
    free(out);

    return 0;
}

使用如下代码编译并运行

1
2
gcc without.c -o without -lm
./without

得到结果如下:

看不见是正常的,别担心

接着再使用KML_FFT库函数进行优化:

编写FFT.c文件,内容如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <stdio.h>
#include <time.h>
#include <math.h>
#include "kfft.h"

int main() {
    int rank = 2
    int *n; 
    struct timespec t1t2;    // 定义初始与结束时间

    n = (int*)kml_fft_malloc(sizeof(int) * rank); 
    n[0] = 128;  // 修改数据规模
    n[1] = 128;  // 修改数据规模

    double *in = (double*)kml_fft_malloc(sizeof(double) * n[0] * n[1]); 
    for (int i = 0; i < n[0] * n[1]; i++) { 
        in[i] = (double)(i % 10);  // 用一些基本的数值初始化输入数据
    } 

    double *out = (double*)kml_fft_malloc(sizeof(double) * n[0] * n[1]); 
    kml_fft_r2r_kind *kind = (kml_fft_r2r_kind*)kml_fft_malloc(sizeof(kml_fft_r2r_kind) * rank); 
    kind[0] = KML_FFT_DHT; 
    kind[1] = KML_FFT_DHT; 

    kml_fft_plan plan; 
    clock_gettime(CLOCK_MONOTONIC, &t1);      // 计算开始时间
    plan = kml_fft_plan_r2r(rank, n, in, out, kind, KML_FFT_ESTIMATE); 
    kml_fft_execute_r2r(plan, in, out); 
    clock_gettime(CLOCK_MONOTONIC, &t2);      // 计算结束时间 

    printf("With KML:\n");
    printf("Time: %11u ns\n", t2.tv_nsec - t1.tv_nsec);  

    printf("\n");

    kml_fft_destroy_plan(plan); 
    kml_fft_free(n); 
    kml_fft_free(kind); 
    kml_fft_free(in); 
    kml_fft_free(out); 

    return 0;
}

使用如下代码编译并运行

1
2
gcc FFT.c -o FFT -L /usr/local/kml/lib -lkfft -I /usr/local/kml/include -pthread
./FFT

得到运行结果如下:

看不见是正常的,别担心

可知KML库对R2R FFT运算的优化效果十分显著

5.实验结果分析

1. 牛顿迭代法求解非线性方程的根

  • 实现了牛顿迭代法,通过简单的非线性方程(如 $f(x)=x2−2$)验证计算结果的准确性。
  • 传统实现和 KML 的计算结果一致,说明 KML 提供的数学库可以准确求解非线性问题。

2. 矩阵-向量乘加运算性能对比

  • 在矩阵规模为 1000×300 的情况下,KML 的矩阵乘加操作(cblas_dgemv)相比手动实现加速显著。
  • 时间数据
    • 手动实现耗时较长,主要由于逐元素计算导致的循环开销。
    • 调用 KML 提供的 BLAS 库后,耗时明显减少,充分利用了向量化和并行化。
  • 结论:KML 在矩阵计算场景中展现了出色的性能提升。

3. KML_VML 向量数学运算优化

  • 比较普通循环与 KML 的矢量数学库处理大规模向量(如 100,000 元素)的性能。
  • 时间数据
    • 传统实现逐元素计算,耗时较长。
    • KML 利用硬件寄存器和 SIMD 指令对数据进行批量处理,极大缩短了运算时间。
  • 分析:适合高频调用数学函数的场景,例如模拟计算和信号处理。

4. 快速傅里叶变换(FFT)优化

  • 对比一般 FFT 算法和 KML_FFT 实现的性能。
  • 时间数据
    • 手动实现 FFT 的时间复杂度较高。
    • KML 的 kml_fft_plan_r2r 方法不仅加速了计算,还简化了实现过程。
  • 结论:KML 的 FFT 模块对高维变换计算尤为高效。

5. LAPACK 在实对称矩阵特征值计算中的优化

  • 测试任务:计算 $10\times10$ 实对称矩阵的特征值和特征向量。
  • 结果分析
    • 传统实现(如 Jacobi 方法)存在显著计算瓶颈。
    • 调用 KML 提供的 LAPACK 接口后,计算时间大幅缩短,进一步展示了 KML 在线性代数中的优化潜力。

通过上述实验,验证了鲲鹏数学库(KML)的性能和适用性。无论是基础数学运算、矩阵计算,还是更复杂的快速变换与特征值问题,KML 的表现均优于传统实现。同时,实验中总结出的优化方法对未来的高性能计算实践提供了有力支持。

6 思考题

使用KML_SVML进行短向量的数学运算优化

KML_SVML是短向量的数学运算,包括幂函数、三角函数、指数函数、双曲函数、对数函数等。 KML_SVML通过Neon指令优化、内联汇编等方法,对输入向量进行批量处理,充分利用了鲲鹏架构下的寄存器特点,实现了在鲲鹏服务器上的性能提升。

请在下面代码空缺处将调用KML_SVML库函数对短向量的数学运算进行优化的代码补充完整

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include<stdio.h>
#include<time.h>
#include<math.h>
#include"ksvml.h"
#define LEN 100000

int main()
{
double src[LEN]={0};
double dst1[LEN]={0};
double dst2[LEN]={0};
float32x4_t src2;
float32x4_t dst;

struct timespec t1,t2; //定义初始与结束时间
printf("Without KML_SVML ");
clock_gettime(CLOCK_MONOTONIC,&t1); //计算开始时间。
for(int i=0;i<LEN;i++)
src[i]=i;
for(int i=0;i<LEN;i++)
dst1[i]=sin(src[i]);
clock_gettime(CLOCK_MONOTONIC,&t2); //计算结束时间。
printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec); //得出目标代码段的执行时间。

printf("With KML_SVML ");
clock_gettime(CLOCK_MONOTONIC,&t1); //计算开始时间。
for(int i=0;i<LEN;i+=4)
{

//请在此处补充调用KML_SVML库函数对短向量的数学运算进行优化的代码

}
clock_gettime(CLOCK_MONOTONIC,&t2); //计算结束时间。
printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec); //得出目标代码段的执行时间。

return 0;

}

答案如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#include<stdio.h>
#include<time.h>
#include<math.h>
#include"ksvml.h"
#define LEN 100000

int main()
{
double src[LEN]={0};
double dst1[LEN]={0};
double dst2[LEN]={0};
float32x4_t src2;
float32x4_t dst;

struct timespec t1,t2; //定义初始与结束时间
printf("Without KML_SVML ");
clock_gettime(CLOCK_MONOTONIC,&t1); //计算开始时间。
for(int i=0;i<LEN;i++)
src[i]=i;
for(int i=0;i<LEN;i++)
dst1[i]=sin(src[i]);
clock_gettime(CLOCK_MONOTONIC,&t2); //计算结束时间。
printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec); //得出目标代码段的执行时间。

printf("With KML_SVML ");
clock_gettime(CLOCK_MONOTONIC,&t1); //计算开始时间。
for(int i=0;i<LEN;i+=4)
{
src2[0]=i;
src2[1]=i+1;
src2[2]=i+2;
src2[3]=i+3;
dst = svml128_sin_f32(src2);
}
clock_gettime(CLOCK_MONOTONIC,&t2); //计算结束时间。
printf("Time:%11u ns\n",t2.tv_nsec-t1.tv_nsec); //得出目标代码段的执行时间。

return 0;

}

输入

1
2
gcc svsin.c -o svsin -I /usr/local/kml/include -L /usr/local/kml/lib -lksvml -lm
./svsin

编译并运行得到运行结果如下:

看不见是正常的,别担心

通过对比可知使用KML_SVML进行短向量的数学运算优化的优化效果十分显著