微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦！

cuSOLVER-cusolverSpScsrlsvqr的设备版本比主机版本慢得多

分类：编程问答

如何解决cuSOLVER-cusolverSpScsrlsvqr的设备版本比主机版本慢得多

我有一些规则建立的稀疏3对角NxN矩阵A，想解决系统Ax=b。为此，我正在使用cusolverSpScsrlsvqr() from cuSolverSp module。大N的设备版本比cusolverSpScsrlsvqrHost()慢很多倍是否可以？例如。对于N = 2 ^ 14，设备的时间为174.1毫秒，主机为3.5毫秒。我正在使用RTX 2060。

代码：

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cusolverSp.h>
#include <cusparse_v2.h>

#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <chrono> 


using namespace std;

void checkCudaCusolverStatus(cusolverStatus_t status,char const* operation) {
    char const *str = "UNKNowN STATUS";
    switch (status) {
    case CUSOLVER_STATUS_SUCCESS:
        str = "CUSOLVER_STATUS_SUCCESS";
        break;
    case CUSOLVER_STATUS_NOT_INITIALIZED:
        str = "CUSOLVER_STATUS_NOT_INITIALIZED";
        break;
    case CUSOLVER_STATUS_ALLOC_Failed:
        str = "CUSOLVER_STATUS_ALLOC_Failed";
        break;
    case CUSOLVER_STATUS_INVALID_VALUE:
        str = "CUSOLVER_STATUS_INVALID_VALUE";
        break;
    case CUSOLVER_STATUS_ARCH_MISMATCH:
        str = "CUSOLVER_STATUS_ARCH_MISMATCH";
        break;
    case CUSOLVER_STATUS_MAPPING_ERROR:
        str = "CUSOLVER_STATUS_MAPPING_ERROR";
        break;
    case CUSOLVER_STATUS_EXECUTION_Failed:
        str = "CUSOLVER_STATUS_EXECUTION_Failed";
        break;
    case CUSOLVER_STATUS_INTERNAL_ERROR:
        str = "CUSOLVER_STATUS_INTERNAL_ERROR";
        break;
    case CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED:
        str = "CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED";
        break;
    case CUSOLVER_STATUS_ZERO_PIVOT:
        str = "CUSOLVER_STATUS_ZERO_PIVOT";
        break;
    }
    cout << left << setw(30) << operation << " " << str << endl;
}

__global__ void fillAB(float *aValues,int *aRowPtrs,int *aColIdxs,float *b,int const n) {
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    if (i >= n) return;
    if (i == 0) {
        float xn = 10 * (n + 1);
        aValues[n * 3] = xn;
        aRowPtrs[0] = 0;
        aRowPtrs[n + 1] = n * 3 + 1;
        aColIdxs[n * 3] = n;
        b[n] = xn * 2;
    }
    float xi = 10 * (i + 1);
    aValues[i * 3 + 0] = xi;
    aValues[i * 3 + 1] = xi + 5;
    aValues[i * 3 + 2] = xi - 5;
    aColIdxs[i * 3 + 0] = i;
    aColIdxs[i * 3 + 1] = i + 1;
    aColIdxs[i * 3 + 2] = i;
    aRowPtrs[i + 1] = 2 + (i * 3);
    b[i] = xi * 2;
}

int main() {
    int const n = (int)pow(2,14);  // <<<<<<<<<<<<<<<<<<<<<<<<<<<<< N HERE
    int const valCount = n * 3 - 2;
    float *const aValues = new float[valCount];
    int *const aRowPtrs = new int[n + 1];
    int *const aColIdxs = new int[valCount];
    float *const b = new float[n];
    float *const x = new float[n];

    float *dev_aValues;
    int *dev_aRowPtrs;
    int *dev_aColIdxs;
    float *dev_b;
    float *dev_x;
    int aValuesSize = sizeof(float) * valCount;
    int aRowPtRSSize = sizeof(int) * (n + 1);
    int aColIdxsSize = sizeof(int) * valCount;
    int bSize = sizeof(float) * n;
    int xSize = sizeof(float) * n;
    cudamalloc((void**)&dev_aValues,aValuesSize);
    cudamalloc((void**)&dev_aRowPtrs,aRowPtRSSize);
    cudamalloc((void**)&dev_aColIdxs,aColIdxsSize);
    cudamalloc((void**)&dev_b,bSize);
    cudamalloc((void**)&dev_x,xSize);
    fillAB<<<1024,(int)ceil(n / 1024.f)>>>(dev_aValues,dev_aRowPtrs,dev_aColIdxs,dev_b,n - 1);
    cudamemcpy(aValues,dev_aValues,aValuesSize,cudamemcpyDevicetoHost);
    cudamemcpy(aRowPtrs,aRowPtRSSize,cudamemcpyDevicetoHost);
    cudamemcpy(aColIdxs,aColIdxsSize,cudamemcpyDevicetoHost);
    cudamemcpy(b,bSize,cudamemcpyDevicetoHost);

    cusolverSpHandle_t handle;
    checkCudaCusolverStatus(cusolverSpCreate(&handle),"cusolverSpCreate");
    cusparseMatDescr_t aDescr;
    cusparseCreateMatDescr(&aDescr);
    cusparseSetMatIndexBase(aDescr,CUSPARSE_INDEX_BASE_ZERO);
    cusparseSetMatType(aDescr,CUSPARSE_MATRIX_TYPE_GENERAL);
    int singularity;
    cusolverStatus_t status;
    cusolverSpScsrlsvqr(handle,n,valCount,aDescr,0.1f,dev_x,&singularity);
    cudaDeviceSynchronize();
    auto t0 = chrono::high_resolution_clock::Now();
    for (int i = 0; i < 10; ++i)
        ////////////////////// CUSOLVER HERE //////////////////////
        status = cusolverSpScsrlsvqr(handle,&singularity);
        //status = cusolverSpScsrlsvqrHost(handle,aValues,aRowPtrs,aColIdxs,b,x,&singularity);
        ///////////////////////////////////////////////////////////
    cudaDeviceSynchronize();
    auto t1 = chrono::high_resolution_clock::Now();
    checkCudaCusolverStatus(status,"cusolverSpScsrlsvqr");
    checkCudaCusolverStatus(cusolverSpDestroy(handle),"cusolverSpDestroy");
    cout << "System solved: " << setw(20) << fixed << right << setprecision(3) << (t1 - t0).count() / 10.0 / 1000000 << " ms" << endl;

    cudamemcpy(x,xSize,cudamemcpyDevicetoHost);
    /*for (int i = 0; i < n; ++i) {
        cout << " " << x[i];
    }*/
    cudaFree(dev_aValues);
    cudaFree(dev_aRowPtrs);
    cudaFree(dev_aColIdxs);
    cudaFree(dev_b);
    cudaFree(dev_x);
    delete[] aValues;
    delete[] aRowPtrs;
    delete[] aColIdxs;
    delete[] b;
    delete[] x;
    cudaDeviceReset();
    return 0;
}

解决方法

我的猜测是这里的问题，它是一个三对角矩阵。我怀疑这可能消除某些并行性方面，这将对GPU cusolver例程有所帮助。除了我在这样的cusparse docs语句中阅读过之外，我对此声明没有任何理由：

例如，三对角矩阵不具有并行性。

确切地说，我不能说什么意思，除了它向我表明，对于一个三对角矩阵，可能需要采用另一种方法。 cusparse provides求解器专门用于三对角线情况。

如果使用其中之一，则在您的测试用例上，可以获得比您的特定主机cusolver示例更快的结果。这是一个示例：

$ cat t48.cu
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cusolverSp.h>
#include <cusparse_v2.h>

#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <chrono>
#include <cassert>
#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL

unsigned long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv,0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}
#ifdef USE_DOUBLE
#define START 3
#define TOL 0.000001
#define THR 0.00001
typedef double mt;
#else
#define START 0
#define TOL 0.01
#define THR 0.1
typedef float mt;
#endif


using namespace std;

void checkCudaCusolverStatus(cusolverStatus_t status,char const* operation) {
    char const *str = "UNKNOWN STATUS";
    switch (status) {
    case CUSOLVER_STATUS_SUCCESS:
        str = "CUSOLVER_STATUS_SUCCESS";
        break;
    case CUSOLVER_STATUS_NOT_INITIALIZED:
        str = "CUSOLVER_STATUS_NOT_INITIALIZED";
        break;
    case CUSOLVER_STATUS_ALLOC_FAILED:
        str = "CUSOLVER_STATUS_ALLOC_FAILED";
        break;
    case CUSOLVER_STATUS_INVALID_VALUE:
        str = "CUSOLVER_STATUS_INVALID_VALUE";
        break;
    case CUSOLVER_STATUS_ARCH_MISMATCH:
        str = "CUSOLVER_STATUS_ARCH_MISMATCH";
        break;
    case CUSOLVER_STATUS_MAPPING_ERROR:
        str = "CUSOLVER_STATUS_MAPPING_ERROR";
        break;
    case CUSOLVER_STATUS_EXECUTION_FAILED:
        str = "CUSOLVER_STATUS_EXECUTION_FAILED";
        break;
    case CUSOLVER_STATUS_INTERNAL_ERROR:
        str = "CUSOLVER_STATUS_INTERNAL_ERROR";
        break;
    case CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED:
        str = "CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED";
        break;
    case CUSOLVER_STATUS_ZERO_PIVOT:
        str = "CUSOLVER_STATUS_ZERO_PIVOT";
        break;
    }
    cout << left << setw(30) << operation << " " << str << endl;
}

__global__ void fillAB(mt *aValues,int *aRowPtrs,int *aColIdxs,mt *b,int const n) {
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    if (i >= n) return;
    if (i == 0) {
        mt xn = 10 * (n + 1);
        aValues[n * 3] = xn;
        aRowPtrs[0] = 0;
        aRowPtrs[n + 1] = n * 3 + 1;
        aColIdxs[n * 3] = n;
        b[n] = xn * 2;
    }
    mt xi = 10 * (i + 1);
    aValues[i * 3 + 0] = xi;
    aValues[i * 3 + 1] = xi + 5;
    aValues[i * 3 + 2] = xi - 5;
    aColIdxs[i * 3 + 0] = i;
    aColIdxs[i * 3 + 1] = i + 1;
    aColIdxs[i * 3 + 2] = i;
    aRowPtrs[i + 1] = 2 + (i * 3);
    b[i] = xi * 2;
}
__global__ void filld3(mt *d,mt *du,mt *dl,mt *aValues,mt *b2,const int n){
        int i = blockDim.x*blockIdx.x+threadIdx.x;
        if ((i > 0) && (i < n-1)){
                dl[i] = aValues[i*3 - 1];
                d[i] = aValues[i*3];
                du[i] = aValues[i*3+1];
        }
        if (i == 0){
                dl[0] = 0;
                d[0]  = aValues[0];
                du[0] = aValues[1];}
        if (i == (n-1)){
                dl[i] = aValues[i*3-1];
                d[i]  = aValues[i*3];
                du[i] = 0;}
        if (i < n) b2[i] = b[i];
}

int main() {
    int const n = (int)pow(2,14);  // <<<<<<<<<<<<<<<<<<<<<<<<<<<<< N HERE
    int const valCount = n * 3 - 2;
    mt *const aValues = new mt[valCount];
    int *const aRowPtrs = new int[n + 1];
    int *const aColIdxs = new int[valCount];
    mt *const b = new mt[n];
    mt *const x = new mt[n];
    mt *const x2= new mt[n];

    mt *dev_aValues;
    int *dev_aRowPtrs;
    int *dev_aColIdxs;
    mt *dev_b;
    mt *dev_x;
    mt *dev_b2,*dev_d,*dev_dl,*dev_du;
    int aValuesSize = sizeof(mt) * valCount;
    int aRowPtrsSize = sizeof(int) * (n + 1);
    int aColIdxsSize = sizeof(int) * valCount;
    int bSize = sizeof(mt) * n;
    int xSize = sizeof(mt) * n;
    cudaMalloc((void**)&dev_aValues,aValuesSize);
    cudaMalloc((void**)&dev_aRowPtrs,aRowPtrsSize);
    cudaMalloc((void**)&dev_aColIdxs,aColIdxsSize);
    cudaMalloc((void**)&dev_b,bSize);
    cudaMalloc((void**)&dev_x,xSize);
    cudaMalloc((void**)&dev_b2,bSize);
    cudaMalloc(&dev_d,n*sizeof(mt));
    cudaMalloc(&dev_dl,n*sizeof(mt));
    cudaMalloc(&dev_du,n*sizeof(mt));
    fillAB<<<1024,(int)ceil(n / 1024.f)>>>(dev_aValues,dev_aRowPtrs,dev_aColIdxs,dev_b,n - 1);
    filld3<<<(n+1023)/1024,1024>>>(dev_d,dev_du,dev_dl,dev_aValues,dev_b2,n);
    cudaMemcpy(aValues,aValuesSize,cudaMemcpyDeviceToHost);
    cudaMemcpy(aRowPtrs,aRowPtrsSize,cudaMemcpyDeviceToHost);
    cudaMemcpy(aColIdxs,aColIdxsSize,cudaMemcpyDeviceToHost);
    cudaMemcpy(b,bSize,cudaMemcpyDeviceToHost);

    cusolverSpHandle_t handle;
    checkCudaCusolverStatus(cusolverSpCreate(&handle),"cusolverSpCreate");
    cusparseMatDescr_t aDescr;
    cusparseCreateMatDescr(&aDescr);
    cusparseSetMatIndexBase(aDescr,CUSPARSE_INDEX_BASE_ZERO);
    cusparseSetMatType(aDescr,CUSPARSE_MATRIX_TYPE_GENERAL);
    int singularity;
    cusolverStatus_t status;
    unsigned long long dt = dtime_usec(0);
#ifdef USE_DOUBLE
    cusolverSpDcsrlsvqr(handle,n,valCount,aDescr,0.1f,dev_x,&singularity);
#else
    cusolverSpScsrlsvqr(handle,&singularity);
#endif
    cudaDeviceSynchronize();
    dt = dtime_usec(dt);
    std::cout << "time: " << dt/(float)USECPSEC << "s" << std::endl;
    auto t0 = chrono::high_resolution_clock::now();
    for (int i = 0; i < 10; ++i)
        ////////////////////// CUSOLVER HERE //////////////////////
#ifdef USE_DEVICE
#ifdef USE_DOUBLE
        status = cusolverSpDcsrlsvqr(handle,&singularity);
#else
        status = cusolverSpScsrlsvqr(handle,&singularity);
#endif
#else
#ifdef USE_DOUBLE
        status = cusolverSpDcsrlsvqrHost(handle,aValues,aRowPtrs,aColIdxs,b,x,&singularity);
#else
        status = cusolverSpScsrlsvqrHost(handle,&singularity);
#endif
#endif
    ///////////////////////////////////////////////////////////
    cudaDeviceSynchronize();
    auto t1 = chrono::high_resolution_clock::now();
    checkCudaCusolverStatus(status,"cusolverSpScsrlsvqr");
    checkCudaCusolverStatus(cusolverSpDestroy(handle),"cusolverSpDestroy");
    cout << "System solved: " << setw(20) << fixed << right << setprecision(6) << (t1 - t0).count() / 10.0 / 1000000 << " ms" << endl;

    cudaMemcpy(x,xSize,cudaMemcpyDeviceToHost);
    /*for (int i = 0; i < n; ++i) {
        cout << " " << x[i];
    }*/
    cusparseHandle_t csphandle;
    cusparseStatus_t  cstat = cusparseCreate(&csphandle);
    assert(cstat == CUSPARSE_STATUS_SUCCESS);
    size_t bufferSize;
#ifdef USE_DOUBLE
    cstat = cusparseDgtsv2_nopivot_bufferSizeExt(csphandle,1,dev_d,&bufferSize);
#else
    cstat = cusparseSgtsv2_nopivot_bufferSizeExt(csphandle,&bufferSize);
#endif
    assert(cstat == CUSPARSE_STATUS_SUCCESS);
    unsigned char *dev_buffer;
    cudaMalloc(&dev_buffer,bufferSize);
    dt = dtime_usec(0);
#ifdef USE_DOUBLE
    cstat = cusparseDgtsv2_nopivot(csphandle,(void *)dev_buffer);
#else
    cstat = cusparseSgtsv2_nopivot(csphandle,(void *)dev_buffer);
#endif
    if(cstat != CUSPARSE_STATUS_SUCCESS) {std::cout << "cusparse solve error: " << (int)cstat  << std::endl;}
    cudaDeviceSynchronize();
    dt = dtime_usec(dt);
    std::cout << "cusparse time: " << (dt*1000.f)/(float)USECPSEC << "ms" << std::endl;
    std::cout << cudaGetErrorString(cudaGetLastError()) << std::endl;
    cudaMemcpy(x2,cudaMemcpyDeviceToHost);
    for (int i = START; i < n; i++) if ((x[i] > THR) && (fabs((x[i] - x2[i])/x[i]) > TOL)) {std::cout << "mismatch at: " << i << " was: " << x2[i] << " should be: " << x[i] << std::endl; return 0;}

    for (int i = 0; i < 40; i++)
            std::cout << x2[i] << "    " << x[i] <<  std::endl;
    cudaFree(dev_aValues);
    cudaFree(dev_aRowPtrs);
    cudaFree(dev_aColIdxs);
    cudaFree(dev_b);
    cudaFree(dev_x);
    delete[] aValues;
    delete[] aRowPtrs;
    delete[] aColIdxs;
    delete[] b;
    delete[] x;
    cudaDeviceReset();
    return 0;
}
$ nvcc -o t48 t48.cu -lcusparse -lcusolver
$ ./t48
cusolverSpCreate               CUSOLVER_STATUS_SUCCESS
time: 0.202933s
cusolverSpScsrlsvqr            CUSOLVER_STATUS_SUCCESS
cusolverSpDestroy              CUSOLVER_STATUS_SUCCESS
System solved:             6.653404 ms
cusparse time: 0.089000ms
no error
-11243.155273    -11242.705078
7496.770508    7496.473145
-3747.185303    -3747.039551
0.685791    0.689445
2083.059570    2082.854004
-1892.308716    -1892.124756
306.474457    306.447662
1103.516846    1103.407104
-1271.085938    -1270.961060
334.883911    334.852417
711.941956    711.870911
-955.718140    -955.624390
321.378174    321.348175
506.580902    506.530060
-764.231689    -764.156799
300.298950    300.270935
382.335785    382.299347
-635.448120    -635.388000
279.217651    279.191864
300.164459    300.135559
-542.869019    -542.817688
259.955475    259.931702
242.390839    242.367310
-473.107758    -473.063080
242.806229    242.785751
199.916733    199.895645
-418.663696    -418.624115
227.637909    227.618652
167.604431    167.586487
-375.002411    -374.966827
214.208069    214.189835
142.353058    142.337738
-339.221130    -339.187653
202.273911    202.255341
122.184746    122.171494
-309.370209    -309.339600
191.615189    191.597580
105.783485    105.771858
-284.096802    -284.068604
182.047958    182.031158
$

注意：

  1. 这里没有正确性或适用性的主张。主要是您的代码，我做了一些修改以进行调查。
  2. 方法之间的结果不完全匹配，但是在float情况下，结果相差约1％。我认为部分原因是因为float的分辨率，但可能还有其他因素。没有进一步的研究，我将没有理由声称一个比另一个“更正确”。
  3. 我使用了nopivot的{​​{1}}变体，因为它似乎表明在2乘方幂情况下，它会更快，这就是您的情况。而且根据我的测试，它更快。
  4. 当我以2 ^ 12大小而不是2 ^ 14大小运行nopivot盒时，它的确在我的GPU（GTX 960）上运行得更快。 YMMV。
  5. 在研究各种事物时，我在代码中添加了其他各种垃圾，因此有点混乱。
  6. 再次，我真的无法解释这个cusolver案。围绕三对角线问题性质的推测就是这样-推测。不过，在我看来，如果挑剔的开发人员找到了为三对角线案例提供（单独）求解器的充分理由，那么可能会有一些合理的理由（例如，可以利用问题的某些方面，事先知道该信息时）。因此使用它们似乎是合理的，在这种情况下，它似乎运行得更快。

版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。

苹果市值2025年有望达4万亿美元
pythonJavaScriptjavaHTMLPHPreactjsC#AndroidCSSNode.jssqlrpython-3.xMysqLjQueryc++pandasFlutterangularIOSdjangolinuxswifttypescript路由器JSON路由器设置无线路由器h3c华三华三路由器设置华三路由器电脑软件教程arraysdocker软件图文教程Cvue.jslaravelspring-boot