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

CUDA-CUBLAS：解决许多3x3密集线性系统的问题

分类：编程问答

如何解决CUDA-CUBLAS：解决许多3x3密集线性系统的问题

我正在尝试使用CUDA 10.1解决大约1200000线性系统（3x3，Ax = B），尤其是使用CUBLAS库。我从this (useful!) post那里得到了提示，并在统一内存版本中重新编写了建议的代码。该算法首先使用cublasgetrfBatched（）执行LU分解，然后连续两次调用cublastrsm（）来求解上三角三角形系统或下三角三角形系统。该代码附在下面。它最多可以正确处理大约10000个矩阵，在这种情况下，执行LU分解（在NVIDIA GeForce 930MX上）大约需要570毫秒，而求解系统大约需要311毫秒。

我的问题是：

  1. 重载问题：为超过1万个矩阵分配内存时会崩溃。为什么？如何改善我的代码以解决整个120万个矩阵？

  2. 时间问题：我的目标是在不到1秒的时间内解决所有这些系统。我目前是否遵循正确的方法？否则有其他建议吗？

  3. 是否可以使用和/或有用，如果可以的话，使用批次1万个矩阵的“流”吗？

代码： 

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <vector>
#include <ctime>
#include <ratio>
#include <chrono>
#include <random>
#include <time.h>
#include <math.h>

// CUDA
#include <cuda.h>
#include <cuda_runtime.h>
#include "device_launch_parameters.h"
#include <cusolverDn.h>

//#include "Utilities.cuh"

using namespace std;
using namespace std::chrono;

/************************************/
/* COEFFICIENT REARRANGING FUNCTION */
/************************************/
void rearrange(double** vec,int* pivotArray,int N,int numMatrices) {
  for (int nm = 0; nm < numMatrices; nm++) {
    for (int i = 0; i < N; i++) {
      double temp = vec[nm][i];
      vec[nm][i] = vec[nm][pivotArray[N*i + nm] - 1];
      vec[nm][pivotArray[N * i + nm] - 1] = temp;
    }
  }
}


/************************************/
/* MAIN  */
/************************************/
int main() {

  const int N = 3; 
  const int numMatrices = 10000; // I want 1200000

  // random generator to fill matrices and coefficients
  random_device device;
  mt19937 generator(device());
  uniform_real_distribution<double> distribution(1.,5.);

  //ALLOCATE MEMORY - using unified memory
  double** h_A;
  cudamallocManaged(&h_A,sizeof(double*) * numMatrices);
  for (int nm = 0; nm < numMatrices; nm++) {
    cudamallocManaged(&(h_A[nm]),sizeof(double) * N * N);
  }

  double** h_b;
  cudamallocManaged(&h_b,sizeof(double*) * numMatrices);
  for (int nm = 0; nm < numMatrices; nm++) {
    cudamallocManaged(&(h_b[nm]),sizeof(double) * N );
  }
  cout << " memory allocated" << endl;

  // FILL MATRICES
  for (int nm = 0; nm < numMatrices; nm++) {
    for (int i = 0; i < N; i++) {
      for (int j = 0; j < N; j++) {
        h_A[nm][j * N + i] = distribution(generator);
      }
    }
  }
  cout << " Matrix filled " << endl;

  // FILL COEFFICIENTS
  for (int nm = 0; nm < numMatrices; nm++) {
    for (int i = 0; i < N; i++) {
      h_b[nm][i] = distribution(generator);
    }
  }
  cout << " Coeff. vector filled " << endl;
  cout << endl;

  // --- CUDA solver initialization
  cublasHandle_t cublas_handle;
  cublasCreate_v2(&cublas_handle);
  int* PivotArray;
  cudamallocManaged(&PivotArray,N * numMatrices * sizeof(int));
  int* infoArray;
  cudamallocManaged(&infoArray,numMatrices * sizeof(int));

  //CUBLAS LU SOLVER
  high_resolution_clock::time_point t1 = high_resolution_clock::Now();
  cublasDgetrfBatched(cublas_handle,N,h_A,PivotArray,infoArray,numMatrices);
  cudaDeviceSynchronize();
  high_resolution_clock::time_point t2 = high_resolution_clock::Now();
  duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
  cout << "It took " << time_span.count() * 1000. << " milliseconds." << endl;


  for (int i = 0; i < numMatrices; i++)
    if (infoArray[i] != 0) {
      fprintf(stderr,"Factorization of matrix %d Failed: Matrix may be singular\n",i);
    }

 // rearrange coefficient 
 // (temporarily on cpu,this step will be on a GPU Kernel as well)
  high_resolution_clock::time_point tA = high_resolution_clock::Now();
  rearrange(h_b,numMatrices);
  high_resolution_clock::time_point tB = high_resolution_clock::Now();
  duration<double> time_spanA = duration_cast<duration<double>>(tB - tA);
  cout << "rearrangement took " << time_spanA.count() * 1000. << " milliseconds." << endl;

//INVERT UPPER AND LOWER TRIANGULAR MATRICES 
  // --- Function solves the triangular linear system with multiple right-hand sides
  // --- Function overrides b as a result 
  const double alpha = 1.f;
  high_resolution_clock::time_point t3 = high_resolution_clock::Now();
  cublasDtrsmBatched(cublas_handle,CUBLAS_SIDE_LEFT,CUBLAS_FILL_MODE_LOWER,CUBLAS_OP_N,CUBLAS_DIAG_UNIT,1,&alpha,h_b,numMatrices);
  cublasDtrsmBatched(cublas_handle,CUBLAS_FILL_MODE_UPPER,CUBLAS_DIAG_NON_UNIT,numMatrices);
  cudaDeviceSynchronize();
  high_resolution_clock::time_point t4 = high_resolution_clock::Now();
  duration<double> time_span2 = duration_cast<duration<double>>(t4 - t3);
  cout << "second step took " << time_span2.count() * 1000. << " milliseconds." << endl;
  
  // --- Free resources
  if (h_A) cudaFree(h_A);
  if (h_b) cudaFree(h_b);
 
  cudaDeviceReset();

  return 0;
}

解决方法

重载问题：为超过1万个矩阵分配内存时会崩溃。为什么？如何改善我的代码以解决整个120万个矩阵？

我认为，代码中的最大问题是，在这些关键分配循环中，您对托管内存的使用效率非常低：

  //ALLOCATE MEMORY - using unified memory
  double** h_A;
  cudaMallocManaged(&h_A,sizeof(double*) * numMatrices);
  for (int nm = 0; nm < numMatrices; nm++) {
    cudaMallocManaged(&(h_A[nm]),sizeof(double) * N * N);
  }

  double** h_b;
  cudaMallocManaged(&h_b,sizeof(double*) * numMatrices);
  for (int nm = 0; nm < numMatrices; nm++) {
    cudaMallocManaged(&(h_b[nm]),sizeof(double) * N );
  }

问题在于，每次对cudaMallocManaged的调用都具有最小粒度。这意味着，如果您请求分配1个字节，则实际上会用掉大约4kbyte的内存（我相信这是linux的分配粒度。看起来您在Windows上，并且我相信Windows的分配粒度可能更大） 。另外，当您启动内核时，这会在托管内存子系统上造成巨大的低效率数据传输负载（内核将在您的cublas调用中启动）。

一种更好的方法是执行单个大分配，而不是循环分配，然后使用指针算法对该分配进行细分。代码看起来像这样：

  //ALLOCATE MEMORY - using unified memory
  double** h_A;
  cudaMallocManaged(&h_A,sizeof(double*) * numMatrices);
  cudaMallocManaged(&(h_A[0]),sizeof(double)*numMatrices*N*N);
  for (int nm = 1; nm < numMatrices; nm++) {
    h_A[nm] = h_A[nm-1]+ N * N;
  }

  double** h_b;
  cudaMallocManaged(&h_b,sizeof(double*) * numMatrices);
  cudaMallocManaged(&(h_b[0]),sizeof(double) * numMatrices * N);
  for (int nm = 1; nm < numMatrices; nm++) {
    h_b[nm] = h_b[nm-1] + N;
  }

这的另一个好处是分配过程的运行速度要快得多。

时间问题：我的目标是在不到1秒的时间内解决所有这些系统。我目前是否遵循正确的方法？否则有其他建议吗？

通过更改代码，我可以在1GB GPU（GeForce GT640）上成功运行，并且具有：

const int numMatrices = 1200000;

具有这样的输出：

$ ./t81
 memory allocated
 Matrix filled
 Coeff. vector filled

It took 70.3032 milliseconds.
rearrangement took 60.02 milliseconds.
second step took 156.067 milliseconds.

您的GPU可能会稍慢一些，但我认为总体时间应该很容易在不到1秒的时间内到达。

使用一批1万个矩阵的“流”是否可能和/或有用，如果可以的话，如何使用？

有了上述更改，我认为您不必为此担心。在这里，流不会帮助计算操作重叠。它们可以帮助复制/计算重叠（尽管在您的GPU上可能不多），但这在带有托管内存的Windows上很难构建。对于Windows使用，如果您想探索copy/compute overlap，我可能建议切换到主机和设备内存的普通CUDA分隔。

顺便说一句，您可能能够获得一组cublas调用，这些调用将通过使用直接反转来更快地完成工作。 CUBLAS具有批处理直接反演方法。对于线性方程组的解决方案，通常我不建议这样做，但是对于一组3x3或4x4的求逆，您可能需要考虑一些问题，您可以使用行列式方法轻松检查奇异性。 Here是一个例子。

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

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