cusolver庫使用(LU分解)

2021-07-29 02:00:39 字數 3362 閱讀 4820

由於網上關於使用cusolver庫的介紹比較少,今天就簡單地根據samples寫了乙個lu分解,來解決稠密矩陣ax=b的例子

需要安裝cuda

首先使用vs建立工程

1.導入庫

第一步 右鍵專案 新增依賴性 生成自定義

第二步 專案屬性 鏈結器 輸入 新增三個依賴項cudart.lib,cusolver.lib,cublas.lib。這是使用cusolver需要的庫。

下文**中呼叫了matrix.h,mat等庫來讀取matlab生成的.mat檔案。這還需要另外新增幾個lib,設定一下環境變數,這部分參考

linearsolverlu這個函式是來自於cuda 的sample cusolverdn_linearsolver。這裡包含了rq分解lu分解等稠密矩陣的gpu實現。
在lapack和 cblas, cusolver等庫中 都會有 lleading dimension 的概念 參考博文

#include "stdafx.h"

#include #include #include #include #include #include#include #include "cublas_v2.h"

#include "cusolverdn.h"

#include #includeint linearsolverlu(

cusolverdnhandle_t handle,

int n,

const double *acopy,

int lda,

const double *b,

double *x)

cudamemcpy(x, b, sizeof(double)*n, cudamemcpydevicetodevice);

cusolverdndgetrs(handle, cublas_op_n, n, 1, a, lda, ipiv, x, n, info);

cudadevicesynchronize();

if (info)

if (buffer)

if (a)

if (ipiv)

return 0;

}int main(int argc, char *argv)

cusolverdnhandle_t handle = null;

cublashandle_t cublashandle = null; // used in residual evaluation

cudastream_t stream = null;

int rowsa = n; // number of rows of a

int colsa = n; // number of columns of a

int nnza = m*n; // number of nonzeros of a

int basea = 0; // base index in csr format

int lda = n; // leading dimension in dense matrix

double *h_a = null;

double *h_x = null; // a copy of d_x

double *h_b = null; // b = ones(m,1)

double *d_a = null; // a copy of h_a

double *d_x = null; // x = a \ b

double *d_b = null; // a copy of h_b

h_a = (double*)malloc(sizeof(double)*lda*colsa);

h_x = (double*)malloc(sizeof(double)*colsa);

h_b = (double*)malloc(sizeof(double)*rowsa);

cusolverdncreate(&handle);

cublascreate(&cublashandle);

cudastreamcreate(&stream);

cusolverdnsetstream(handle, stream);

cublassetstream(cublashandle, stream);

cudamalloc((void **)&d_a, sizeof(double)*lda*colsa);

cudamalloc((void **)&d_x, sizeof(double)*colsa);

cudamalloc((void **)&d_b, sizeof(double)*rowsa);

cudamemcpy(d_a, a0, sizeof(double)*lda*colsa, cudamemcpyhosttodevice);

cudamemcpy(d_b, b0, sizeof(double)*rowsa, cudamemcpyhosttodevice);

clock_t start, finish;

double duration;

start = clock();

linearsolverlu(handle, colsa, d_a, lda, d_b, d_x);

finish = clock();

duration = (double)(finish - start) / clocks_per_sec;

cudamemcpy(h_x, d_x, sizeof(double)*colsa, cudamemcpydevicetohost);

double max_r = 0;

for (int i = 0; i < n; ++i)

printf("ax=b\n");

printf("rows=%d,cols=%d\n", rowsa, colsa);

printf("%f seconds\n", duration);

printf("max(x_gpu-xcpu)= %lf\n", max_r);

getchar();

matclose(a_file);

matclose(b_file);

cudafree(d_a);

cudafree(d_b);

cudafree(d_x);

mxfree(a0);

mxfree(b0);

mxfree(x);

free(h_a);

free(h_b);

free(h_x);

return 0;

}

結果如下。稠密矩陣a的行數和列數都為4384,gpu解決需要0.79s,而matlab這個過程需要1.2s。和matlab結果差為0.000003.

測試cpu為i7 4710hq,gpu為gtx880m。

LU分解和轉置

我們通過消元法解出 ax b 在消元過程中,我們要把a化成上三角形式,我們稱之為 u 而 e 就是描述 a 與 u 關係的矩陣,我們有 ea u 這是之前學過的內容,再熟悉不過了,但是我們的行交換矩陣 e 並不是乙個漂亮的形式,為啥這麼說呢!看個例幾 a left begin 2 4 2 4 10 ...

C C 使用Lu擴充套件動態庫

lu程式設計 c c 使用lu擴充套件動態庫 1 說明 lu32.h,相信你會找到並正確使用這幾個檔案。用c c 編譯器建立乙個控制台應用程式,複製本文的例子 直接編譯執行即可。2 關於lu擴充套件動態庫的使用 lu擴充套件動態庫只有乙個輸出函式 lu擴充套件動態庫唯一的輸出函式 hlu lu32....

c 矩陣求逆的lu分解實現

c 矩陣求逆的lu分解實現 初學c 嘗試利用基礎知識編寫矩陣求逆。但發現求解伴隨陣的過程非常複雜且難以實現,而我正好看到一篇求三角陣伴隨矩陣的文章,故嘗試程式設計實現。在這種方法下,計算量明顯減小,實現方法,思路適合初學者。參考文獻 三角形矩陣求伴隨矩陣的一種方法 曾月新 求逆矩陣思路 1.求矩陣的...