本发明属于高性能计算，尤其涉及一种基于gpu并行加速的对称条带矩阵三对角化变换方法。
背景技术：
1、对称矩阵特征值分解是现在高性能数据计算的重要课题，涉及到人工智能、科学计算等多个领域。对称矩阵特征值分解目前常用的方法是三阶段分解算法，分为对称矩阵条带化阶段、对称条带化矩阵三对角化阶段和求特征值阶段，分别完成对称矩阵的条带化分解、将条带化对称矩阵分解为三对角化和使用qr迭代算法或者分治法求出最终的特征值三个任务。现有算法的性能瓶颈体现在前两个阶段，在常用算法包cuda cusolver中所花时间占95%以上。
2、对称条带化矩阵三对角化分解，其常用的方法为凸块追逐(bulge chasing)。但此方法在目前通常认为是内存瓶颈，大家认为其算法瓶颈主要存在于内存搬运上，是内存瓶颈问题不是计算瓶颈问题，所以到目前为止并没有在gpu上进行实现的版本。
技术实现思路
1、随着硬件加速技术的不断发展，英伟达gpu在人工智能、并行计算等领域取得了巨大成功。其simd(single instruction/multiple data单指令多数据)架构能够同时运行多个线程块，每个线程块中又可以同时运行多个线程，这样极高的并行性为多种算法的并行优化提供了巨大潜力。凸块追逐算法虽然在单趟凸块追逐过程中只能进行串行计算，但其多趟凸块追逐过程具有巨大并行潜力，本发明实现了一种新颖高效的gpu并行算法，对凸块追逐整个过程进行了并行实现。对比magma(美国田纳西大学开源的线性代数软件库)与cusolver（英伟达发布的线性代数软件库），分别具有7.5倍和10.1倍的性能提升。
2、为解决上述技术问题，本发明的一种基于gpu并行加速的对称条带矩阵三对角化变换方法的具体技术方案如下：
3、一种基于gpu并行加速的对称条带矩阵三对角化变换方法，针对n*n的对称条带化矩阵a，其条带化宽度为b，进行三条带化变换，包括以下步骤：
4、步骤1：在gpu上创建大小为2*b*n的全局数组ds[2*b*n]，用于步骤2中各线程块进行数据同步；设置变量i，j，分别表示待复制矩阵a中元素的行数和列数，将矩阵a中满足j≤i≤j+b的元素复制到数组ds中；
5、具体的，初始化j=0。将矩阵a中满足j≤i≤j+b的元素复制到数组ds中的具体步骤如下：
6、步骤1.1：初始化i=0，令待复制数据的个数count=min(n,j+b+1)-j，跳转到步骤1.2；
7、步骤1.2：判断i是否小于count，若是则令ds[i+j*2*b]=a[i+j+j*n]，i=i+1，跳转到步骤1.2，否则跳转到步骤1.3；
8、步骤1.3：令j=j+1,判断j是否小于n，若是则跳转到步骤1.1，否则跳转到步骤2。
9、步骤2：在gpu上创建n-2个线程块block，每个线程块分别处理从第1列到第n-2列开始的凸块追逐，每个线程块包括k*k个线程，k*k个线程一起处理单次凸块追逐；初始化数组com[n-2]，用于各线程块进行处理位置同步，用于记录每个线程块进行凸块追逐的当前位置oprow，初始化时com[n-2]={0}；
10、步骤3：每个线程块单独处理对应的凸块追逐；令每个线程块处理的当前位置为oprow，每个线程块的编号记为bindx，每个凸块追逐的householder变换过程涉及到3个部分，分别为b1、s和b2；b1的行数和列数分别为b1row和b1col，s的行数和列数分别为srow和scol， b2的行数和列数分别为b2row和b2col；具体的，在一次凸块追逐中，b1是矩阵a中行数为oprow到oprow+b1row-1，列数为oprow-b1col到oprow-1的部分；s是矩阵a中行数为oprow到oprow+srow-1，列数为oprow到oprow+scol-1的部分；b2是矩阵a中行数为oprow+srow到oprow+srow+b2row-1，列数为oprow到oprow+b2col-1的部分；在每个线程块中创建共享内存数组s_b1[b*b]、s_s[b*b]、s_b2[b*b]分别用于存放b1，s和b2三部分；令b2col=1，oprow=bindx+1；凸块追逐的具体过程如下：
11、步骤3.1：若当前线程块处理事件满足以下情况：
12、bindx≠0
13、且oprow+2*b＞com[binx-1]
14、则进行循环判断，否则跳转至3.2；其中com[binx-1] 表示数组中第binx个元素；
15、步骤3.2：
16、b1col=b2col
17、b1row=min(b，n- oprow)
18、b2col=scol =srow= b1row
19、b2row=min(b， n-srow-oprow)
20、复制数组ds中的相应部分到s_b1、s_s和s_b2中，跳转到步骤3.3;
21、步骤3.3：首先对s_b1的第1列进行householder变换，得到对应的householder向量u，令
22、分别对s_b1、s_s、s_b2三部分进行变换，其中，i为单位矩阵，大小为b1row *b1row ；跳转到步骤3.4；
23、步骤3.4：复制s_b1、s_s和s_b2到数组ds中的相应部分，跳转到步骤3.5；
24、步骤3.5：令
25、oprow=oprow+b
26、com[binx]= oprow
27、判断b2row是否小于2，若是，则进入步骤4，若否，则跳转到步骤3.1；其中com[binx] 表示数组中第binx+1个元素。
28、步骤4：等待所有线程块都完成自己的凸块追逐，将数组ds中的第1行复制到矩阵a的主对角线上，将数组ds中的第2行复制到矩阵a的上次对角线与下次对角线得到三条带化矩阵，完成整个凸块追逐过程。
29、具体的，所述步骤3.2具体如下：
30、步骤3.2.1：设置变量ds_b1、ds_s和ds_b2分别为待复制ds中b1、s和b2的对应位置；令ds_b1=ds+b1col + (oprow - b1col)*2*b、ds_s=ds+ oprow*2*b、ds_b2=ds+srow+oprow*2*b；令i，j分别为需要复制的数据的行数和列数，令j=0，跳转至步骤3.2.2；
31、步骤3.2.2：判断j是否小于b1col，若是则令i=0，跳转至3.2.3；否则令j=0，跳转至3.2.4；
32、步骤3.2.3：判断i是否小于b1row，若是则令s_b1[i+j*b]= ds_b1[i+(i-j)*b]，i++，跳转到3.2.3；否则令j=j+1，跳转至3.2.2；
33、步骤3.2.4：判断j是否小于scol，若是则令i=0跳转至3.2.5；否则令j=0，跳转至3.2.6；
34、步骤3.2.5：判断i是否小于srow，若是则令s_s[i+j*b]= ds_s[i+(i-j)*b]，i++，跳转到3.2.5；否则令j=j+1，跳转至3.2.4；
35、步骤3.2.6：判断j是否小于b2col，若是则令i=0，跳转至3.2.7；否则完成步骤3.2，跳转至步骤3.3；
36、步骤3.2.7：判断i是否小于b2row，若是则令s_b2[i+j*b]= ds_b2[i+(i-j)*b]，i++，跳转到3.2.7；否则令j=j+1，跳转至3.2.6；
37、具体的，所述步骤3.4恰好和步骤3.2相反，具体如下：
38、步骤3.4.1：设置变量ds_b1、ds_s和ds_b2分别为待复制的ds中b1、s和b2的相应位置；令ds_b1=ds+b1col + (oprow - b1col) *2*b、ds_s=ds+ oprow*2*b、ds_b2=ds+srow+oprow*2*b;令i，j分别为需要复制的数据的行数和列数，令j=0，跳转至步骤3.4.2；
39、步骤3.4.2：判断j是否小于b1col，若是则令i=0，跳转至3.4.3；否则令j=0，跳转至3.4.4；
40、步骤3.4.3：判断i是否小于b1row，若是则令ds_b1[i+(i-j)*b]=s_b1[i+j*b]，i++，跳转到3.4.3；否则令j=j+1，跳转至3.4.2；
41、步骤3.4.4：判断j是否小于scol，若是则令i=0，跳转至3.4.5；否则令j=0，跳转至3.4.6；
42、步骤3.4.5：判断i是否小于srow，若是则令ds_s[i+(i-j)*b]=s_s[i+j*b]，i++，跳转到3.4.5；否则令j=j+1，跳转至3.4.4；
43、步骤3.4.6：判断j是否小于b2col，若是则令i=0，跳转至3.4.7，否则完成步骤3.4，跳转至步骤3.5；
44、步骤3.4.7：判断i是否小于b2row，若是则令ds_b2[i+(i-j)*b]=s_b2[i+j*b]，i++，跳转到3.4.7；否则令j=j+1，跳转至3.4.6；
45、具体的，所述步骤3.3具体如下：
46、步骤3.3.1：将每个线程块中创建的k*k个线程划分为k个线程束，其编号为p(0≤p<k)；每个线程束中含有k个线程,其编号为q(0≤q<k)；跳转至3.3.2；
47、步骤3.3.2：令变量col为待进行householder变换的s_b1的列数，令col=p；跳转至步骤3.3.3；
48、步骤3.3.3：判断col是否小于b1col，若是则线程束p中的k个线程对s_b1的col列进行householder变换，令col=col+k，跳转至3.3.3；否则跳转至步骤3.3.4；
49、步骤3.3.4：令变量col为待进行householder变换的s_s的列数，令col=p；跳转至步骤3.3.5；
50、步骤3.3.5：判断col是否小于scol，若是则线程束p中的k个线程对s_s的col列进行householder变换，令col=col+k，跳转至3.3.5；否则跳转至步骤3.3.6；
51、步骤3.3.6：令变量row为待进行householder变换的s_s的行数，令row=p；跳转至步骤3.3.7；
52、步骤3.3.7：判断row是否小于b1row，若是则线程束p中的k个线程对s_s的row行进行householder变换，令row=row+k，跳转至3.3.7；否则跳转至步骤3.3.8；
53、步骤3.3.8：令变量col为待进行householder变换的s_b2的列数，令col=p；跳转至步骤3.3.9；
54、步骤3.3.9：判断col是否小于b2col，若是则线程束p中的k个线程对s_b2的col列进行householder变换，令col=col+k，跳转至3.3.9；否则跳转至步骤3.4；
55、具体的，所述步骤4具体如下：
56、步骤4.1：令i=0，为待复制数组ds的列数；
57、步骤4.2：判断i是否小于n，若是则跳转至步骤4.3；否则结束；
58、步骤4.3：判断i是否等于0且等于n-1，若不是则令a[i+i*n]=ds[i*2*b]，a[i+1+i*n]=ds[1+i*2*b]，a[i-1+i*n]=ds[1+i*2*b]，i=i+1跳转至步骤4.2；判断i是否等于n-1，若是令a[i+i*n]=ds[i*2*b]，a[i-1+i*n]=ds[1+i*2*b]，i=i+1，跳转至步骤4.2；判断i是否等于0，若是令a[i+i*n]=ds[i*2*b]，a[i+1+i*n]=ds[1+i*2*b]，i=i+1，跳转至步骤4.2。
59、本发明的有益效果如下：
60、本发明提供了一种基于gpu并行加速的对称条带矩阵三对角化变换方法，该方法首次在gpu上实现凸块追逐，并实现了并行加速，与现有的对称条带矩阵三对角化方法相比，能够大幅提升gpu处理性能，计算效率高。
1.一种基于gpu并行加速的对称条带矩阵三对角化变换方法，针对n*n的对称条带化矩阵a，其条带化宽度为b，进行三条带化变换，其特征在于，包括以下步骤：
2.根据权利要求1所述的一种基于gpu并行加速的对称条带矩阵三对角化变换方法，其特征在于，所述步骤1具体如下：
3.根据权利要求1所述的一种基于gpu并行加速的对称条带矩阵三对角化变换方法，其特征在于，所述步骤3具体如下：
4.根据权利要求3所述的一种基于gpu并行加速的对称条带矩阵三对角化变换方法，其特征在于，所述步骤3.2具体如下：
5.根据权利要求3所述的一种基于gpu并行加速的对称条带矩阵三对角化变换方法，其特征在于，所述步骤3.4恰好和步骤3.2相反，具体如下：
6.根据权利要求3所述的一种基于gpu并行加速的对称条带矩阵三对角化变换方法，其特征在于，所述步骤3.3具体如下：
7.根据权利要求1所述的一种基于gpu并行加速的对称条带矩阵三对角化变换方法，其特征在于，所述步骤4具体如下：步骤4.1：令i=0，为待复制数组ds的列数；
8.根据权利要求1所述的一种基于gpu并行加速的对称条带矩阵三对角化变换方法，其特征在于，所述gpu具有单指令多数据simd架构。
