CUDA求解特征值和特征向量 使用cusolver库

         项目中遇到了求解复数Hermit矩阵的特征值分解问题（使用MUSIC方法进行DOA估计的GPU工程化实现），网上已经有使用GSL科学计算库（http://www.gnu.org/software/gsl/）完成特征值和特征向量求解问题，也可以使用QR算法和Jacobi算法等数值分析方法自己编写，现在想使用CUDA在GPU平台上并行完成特征值和特征向量求解问题。

        之前想先理解数值分析算法自己写并行程序，后来发现CUDA已经提供库函数完成此功能——cusolverDnCheevd()，通过例程和尝试使用，现将学习笔记记录如下：

系统环境：

win10 64位
VS2017
CUDA 10

一、函数说明

        参考文档 CUSOLVER_Library.pdf 中5.3.9节有关Hermit矩阵or对称矩阵特征值分解方法。GPU求解特征值和特征向量包括两个主要函数：缓存空间计算函数和特征分解函数。

首先，缓存空间计算函数用于提前计算特征分解所需的buffer大小保存在lwork变量中； 
然后，在显存中通过cudaMalloc函数分配lwork大小的空间供特征分解函数使用； 
最后，启动特征分解函数求解特征值和特征向量。 	
注：两个函数根据处理的矩阵是实数矩阵还是复数矩阵，以及要完成的是单精度计算还是双精度计算又可以分为四个具体函数，因而一共8个函数。

1. 缓存空间计算函数：

2.特征分解函数：

变量解释：

0. 实复数与单双精度：S 实数单精度D 实数双精度C 复数单精度Z复数双精度
1. 函数句柄：用于声明函数句柄，其中Dn是Dense缩写指非稀疏矩阵

cusolverDnHandle_t cusolverH = NULL;
cusolverDnCreate(&cusolverH);

2. 分解模式：分为只计算特征值模式 和 特征值特征向量模式

cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; //结果有特征值和特征向量
cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_NOVECTOR; //结果只有特征值

3. 存储模式：由于是方阵是对称或共轭的，只需要取用上三角或下三角即可，但习惯上我们还是建议输入的方阵全部存数据，而不能只存上下三角。

cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;//下三角
cublasFillMode_t uplo = CUBLAS_FILL_MODE_UPPER;//上三角

4. 方阵大小n：表示方阵行数列数，行列是一个数，内存存储
5. 方阵数据A：是n*n个数据，按照列优先方式存储，如果选择输出特征向量，该空间会存储归一化的特征向量，每一列对应w中保存的特征值。显存存储。
6. 子阵行数lda：有在大方阵中求解小方阵的特征分解的计算需求，由于矩阵采用列优先的方式存储，确定lda就能自动建立子阵。默认求解大方阵特征分解，只需使lda = n 即可。内存存储
   
7. 特征值结果w：计算的特征值会保存在该变量中，且是升序排列。显存存储。
8. 计算缓存空间work：在显存中分配的lwork大小的计算缓存空间。
9. 缓存空间大小lwork: 内存存储
10. 状态码devInfor：0表示计算正确，-i表示第i个参数有误，i表示第i个非对角线元素未收敛于0。

二、例程操作

1. 实数对称矩阵单精度特征分解例程：

#include <stdio.h>
#include <stdlib.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cusolverDn.h>

int main()
{
	cusolverDnHandle_t cusolverH = NULL;
	const int m = 3;
	const int lda = m;
	float *A = (float*)malloc(lda*m * sizeof(float));
	A[0] = 3.5;
	A[1] = 0.5;
	A[2] = 0;
	A[3] = 0.5;
	A[4] = 3.5;
	A[5] = 0;
	A[6] = 0;
	A[7] = 0;
	A[8] = 2;
	float W[m]; // eigenvalues最终保存结果
	int info_gpu = 0;//计算状态保存
	// 步骤1：声明句柄
	cusolverDnCreate(&cusolverH);
	// 步骤2：分配显存空间
	float *d_A = NULL; cudaMalloc((void**)&d_A, sizeof(float) * lda * m);//声明Hermite矩阵（与计算后的特征向量为同一空间）
	float *d_W = NULL; cudaMalloc((void**)&d_W, sizeof(float) * m);//声明特征值存储空间
	int *devInfo = NULL; cudaMalloc((void**)&devInfo, sizeof(int));//声明计算结果状态空间
	cudaMemcpy(d_A, A, sizeof(float) * lda * m, cudaMemcpyHostToDevice);//数据拷贝
	// 步骤3：申请计算缓存空间，并在显存中申请该空间
	float *d_work = NULL;
	int lwork = 0;
	cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvalues and eigenvectors.
	cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
	cusolverDnSsyevd_bufferSize(cusolverH, jobz, uplo, m, d_A, lda, d_W, &lwork);//计算evd计算所需存储空间,保存到lwork中
	cudaMalloc((void**)&d_work, sizeof(float)*lwork);
	// 步骤4：特征分解
	cusolverDnSsyevd(cusolverH, jobz, uplo, m, d_A, lda, d_W, d_work, lwork, devInfo);
	cudaDeviceSynchronize();
	//步骤5：数据读回
	cudaMemcpy(A, d_A, sizeof(float)*lda*m, cudaMemcpyDeviceToHost);
	cudaMemcpy(W, d_W, sizeof(float)*m, cudaMemcpyDeviceToHost);
	cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost);

	printf("%d\n", info_gpu);
	printf("eigenvalue = (matlab base-1), ascending order\n");
	for (int i = 0; i < m; i++) {
		printf("W[%d] = %E\n", i + 1, W[i]);
	}
	for (size_t i = 0; i < m; i++)
	{
		for (size_t j = 0; j < m; j++)
		{
			printf("%.4f\n", A[i + j*m]);
		}
	}

    return 0;
}

2. 复数Hermit矩阵双精度特征分解例程：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <stdlib.h>
#include <cusolverDn.h>

int main()
{
	cusolverDnHandle_t cusolverH = NULL;
	const int m = 4;
	const int lda = m;

	cuDoubleComplex *A = (cuDoubleComplex*)malloc(lda*m*sizeof(cuDoubleComplex));
	A[0] = make_cuDoubleComplex(1.9501e2,0);
	A[1] = make_cuDoubleComplex(0.2049e2, 0.1811e2);
	A[2] = make_cuDoubleComplex(0.5217e2, 0.3123e2);
	A[3] = make_cuDoubleComplex(0.2681e2, 0.3998e2);
	A[4] = make_cuDoubleComplex(0.2049e2, -0.1811e2);
	A[5] = make_cuDoubleComplex(1.8272e2, 0);
	A[6] = make_cuDoubleComplex(0.5115e2, -0.0987e2);
	A[7] = make_cuDoubleComplex(0.4155e2, -0.0435e2);
	A[8] = make_cuDoubleComplex(0.5217e2, -0.3123e2);
	A[9] = make_cuDoubleComplex(0.5115e2, 0.0987e2);
	A[10] = make_cuDoubleComplex(2.3984e2, 0);
	A[11] = make_cuDoubleComplex(0.4549e2, -0.0510e2);
	A[12] = make_cuDoubleComplex(0.2681e2, -0.3998e2);
	A[13] = make_cuDoubleComplex(0.4155e2, 0.0435e2);
	A[14] = make_cuDoubleComplex(0.4549e2, 0.0510e2);
	A[15] = make_cuDoubleComplex(2.2332e2, 0);

	//1.9501 + 0.0000i   0.2049 - 0.1811i   0.5217 - 0.3123i   0.2681 - 0.3998i
	//	0.2049 + 0.1811i   1.8272 + 0.0000i   0.5115 + 0.0987i   0.4155 + 0.0435i
	//	0.5217 + 0.3123i   0.5115 - 0.0987i   2.3984 + 0.0000i   0.4549 + 0.0510i
	//	0.2681 + 0.3998i   0.4155 - 0.0435i   0.4549 - 0.0510i   2.2332 + 0.0000i

	double W[m]; // eigenvalues最终保存结果
	int info_gpu = 0;//计算状态保存
	// step 1: create cusolver/cublas handle
	cusolverDnCreate(&cusolverH);
	// step 2: copy A and B to device
	cuDoubleComplex *d_A = NULL; cudaMalloc((void**)&d_A, sizeof(cuDoubleComplex) * lda * m);//声明Hermite矩阵（与计算后的特征向量为同一空间）
	double *d_W = NULL; cudaMalloc((void**)&d_W, sizeof(double) * m);//声明特征值存储空间
	int *devInfo = NULL;cudaMalloc((void**)&devInfo, sizeof(int));//声明计算结果状态空间

	cudaMemcpy(d_A, A, sizeof(cuDoubleComplex) * lda * m, cudaMemcpyHostToDevice);//数据拷贝
	// step 3: query working space of syevd
	cuDoubleComplex *d_work = NULL;
	int lwork = 0;
	cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvalues and eigenvectors.
	cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
	cusolverDnZheevd_bufferSize(cusolverH, jobz, uplo, m, d_A, lda, d_W, &lwork);//计算evd计算所需存储空间,保存到lwork中
	cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex)*lwork);

	// step 4: compute spectrum
	cusolverDnZheevd(cusolverH, jobz, uplo, m, d_A, lda, d_W, d_work, lwork, devInfo);
	cudaDeviceSynchronize();
	cudaMemcpy(A, d_A, sizeof(cuDoubleComplex)*lda*m,cudaMemcpyDeviceToHost);
	cudaMemcpy(W, d_W, sizeof(double)*m, cudaMemcpyDeviceToHost);
	cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost);

	printf("%d\n", info_gpu);
	printf("eigenvalue = (matlab base-1), ascending order\n");
	for (int i = 0; i < m; i++) {
		printf("W[%d] = %E\n", i + 1, W[i]);
	}
	for (size_t i = 0; i < 4; i++)
	{
		for (size_t j = 0; j < 4; j++)
		{
			printf("%.4f + %.4f j\n", A[i * 4 + j].x, A[i * 4 + j].y);
		}
	}
	cudaFree(d_A);
	cudaFree(d_W);
	cudaFree(devInfo);
	cudaFree(d_work);
	cusolverDnDestroy(cusolverH);

    return 0;
}

三、注意事项

除了在.cu文件中引入头文件**#include <cusolverDn.h>，也要记得在项目属性中添加依赖库cusolver.lib**