devbox
我家小花儿
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

一个使用cublasSgetrfBatched计算逆矩阵的示例_cublas 求逆

作者：我家小花儿 | 2024-03-15 23:51:39

赞

踩

cublas 求逆

  1. //mat_inv.cu
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <cublas_v2.h>
  5. #include<math.h>
  6.  
  7.  
  8. #define cudacall(call)                                                                                                          \
  9.     do                                                                                                                          \
  10.     {                                                                                                                           \
  11.         cudaError_t err = (call);                                                                                               \
  12.         if(cudaSuccess != err)                                                                                                  \
  13.         {                                                                                                                       \
  14.             fprintf(stderr,"CUDA Error:\nFile = %s\nLine = %d\nReason = %s\n", __FILE__, __LINE__, cudaGetErrorString(err));    \
  15.             cudaDeviceReset();                                                                                                  \
  16.             exit(EXIT_FAILURE);                                                                                                 \
  17.         }                                                                                                                       \
  18.     }                                                                                                                           \
  19.     while (0)
  20.  
  21. #define cublascall(call)                                                                                        \
  22.     do                                                                                                          \
  23.     {                                                                                                           \
  24.         cublasStatus_t status = (call);                                                                         \
  25.         if(CUBLAS_STATUS_SUCCESS != status)                                                                     \
  26.         {                                                                                                       \
  27.             fprintf(stderr,"CUBLAS Error:\nFile = %s\nLine = %d\nCode = %d\n", __FILE__, __LINE__, status);     \
  28.             cudaDeviceReset();                                                                                  \
  29.             exit(EXIT_FAILURE);                                                                                 \
  30.         }                                                                                                       \
  31.                                                                                                                 \
  32.     }                                                                                                           \
  33.     while(0)
  34.  
  35.  
  36. //
  37. //
  38.  
  39. void invert(float** src, float** dst, int n, int batchSize)
  40. {
  41.     	cublasHandle_t handle;
  42.     	cublascall(cublasCreate_v2(&handle));
  43.  
  44.     	int *P, *INFO;
  45.  
  46.     	cudacall(cudaMalloc(&P, n * batchSize * sizeof(int)));
  47.     	cudacall(cudaMalloc(&INFO,  batchSize * sizeof(int)));
  48.  
  49.     	int lda = n;
  50.  
  51.     	float **A = (float **)malloc(batchSize*sizeof(float *));
  52.     	float **A_d, *A_dflat;
  53.  
  54.     	cudacall(cudaMalloc(&A_d,batchSize*sizeof(float *)));
  55.     	cudacall(cudaMalloc(&A_dflat, n*n*batchSize*sizeof(float)));
  56.  
  57. 	A[0] = A_dflat;
  58.     	for (int i = 1; i < batchSize; i++)
  59.       		A[i] = A[i-1]+(n*n);
  60.  
  61.     	cudacall(cudaMemcpy(A_d,A,batchSize*sizeof(float *),cudaMemcpyHostToDevice));
  62.    
  63.  	for (int i = 0; i < batchSize; i++)
  64.       		cudacall(cudaMemcpy(A_dflat+(i*n*n), src[i], n*n*sizeof(float), cudaMemcpyHostToDevice));
  65.  
  66.  
  67.     	cublascall(cublasSgetrfBatched(handle,n,A_d,lda,P,INFO,batchSize));
  68.  
  69.  
  70.     	int INFOh[batchSize];
  71.     	cudacall(cudaMemcpy(INFOh,INFO,batchSize*sizeof(int),cudaMemcpyDeviceToHost));
  72.  
  73.     	for (int i = 0; i < batchSize; i++)
  74.       		if(INFOh[i]  != 0)
  75.       		{
  76.         		fprintf(stderr, "Factorization of matrix %d Failed: Matrix may be singular\n", i);
  77.         		cudaDeviceReset();
  78.         		exit(EXIT_FAILURE);
  79.       		}
  80.  
  81.     	float **C = (float **)malloc(batchSize*sizeof(float *));
  82.     	float **C_d, *C_dflat;
  83.  
  84.     	cudacall(cudaMalloc(&C_d,batchSize*sizeof(float *)));
  85.     	cudacall(cudaMalloc(&C_dflat, n*n*batchSize*sizeof(float)));
  86.     	C[0] = C_dflat;
  87.     	for (int i = 1; i < batchSize; i++)
  88.       		C[i] = C[i-1] + (n*n);
  89.     	cudacall(cudaMemcpy(C_d,C,batchSize*sizeof(float *),cudaMemcpyHostToDevice));
  90.     	cublascall(cublasSgetriBatched(handle,n,(const float **)A_d,lda,P,C_d,lda,INFO,batchSize));
  91.  
  92.     	cudacall(cudaMemcpy(INFOh,INFO,batchSize*sizeof(int),cudaMemcpyDeviceToHost));
  93.  
  94.     	for (int i = 0; i < batchSize; i++)
  95. 	      	if(INFOh[i] != 0)
  96.       		{
  97.         		fprintf(stderr, "Inversion of matrix %d Failed: Matrix may be singular\n", i);
  98.         		cudaDeviceReset();
  99.       			exit(EXIT_FAILURE);
  100.       		}
  101.     	for (int i = 0; i < batchSize; i++)
  102.       	cudacall(cudaMemcpy(dst[i], C_dflat + (i*n*n), n*n*sizeof(float), cudaMemcpyDeviceToHost));
  103.     	
  104. 	cudaFree(A_d); cudaFree(A_dflat); free(A);
  105. 	cudaFree(C_d); cudaFree(C_dflat); free(C);
  106.     	cudaFree(P); cudaFree(INFO); cublasDestroy_v2(handle);
  107. }
  108.  
  109.  
  110.  
  111. //
  112. //
  113.  
  114.  
  115.  
  116. __global__ void MatrixMulKernel(float *Md, float *Nd, float *Pd, int Width) 
  117. {
  118.     //2D Thread ID
  119.     int col = blockIdx.x*blockDim.x+threadIdx.x;
  120.     int row = blockIdx.y*blockDim.y+threadIdx.y;
  121.  
  122.     //Pvalue stores the Pd element that is computed by the thread
  123.     float Pvalue = 0;
  124. 	if(col<Width && row < Width)
  125. 	{
  126. 		    for(int k = 0; k < Width ; ++k) 
  127. 		    {
  128. 		        float Mdelement = Md[row*Width + k];
  129. 		        float Ndelement = Nd[k*Width + col];
  130. 		        Pvalue += (Mdelement*Ndelement);
  131. 		
  132. 		    }
  133. 		    Pd[row*Width + col] = Pvalue;
  134. 	}
  135. }
  136.  
  137.  
  138.  
  139. void mul(float* M,float* N,int Width)
  140. {
  141. 		
  142. 	float * P = (float *) malloc(Width*Width*sizeof(float));
  143. 	float *Md, *Nd, *Pd;
  144.  
  145.  
  146.  
  147. 	unsigned long int size = Width*Width*sizeof(float);
  148.   
  149.  
  150.     //Transfer M and N to device memory
  151.     	cudaMalloc((void**)&Md, size);
  152.     	cudaMemcpy(Md,M,size,cudaMemcpyHostToDevice);
  153.  
  154.     	cudaMalloc((void**)&Nd, size);
  155.     	cudaMemcpy(Nd,N,size,cudaMemcpyHostToDevice);
  156.  
  157.     	//Allocate P on the device
  158.     	cudaMalloc((void**)&Pd,size);
  159.  
  160.     	//Setup the execution configuration
  161.     	dim3 dimBlock(Width,Width);
  162.     	dim3 dimGrid(1,1);
  163.  
  164.  
  165. 	if (Width*Width > 1024)
  166. 	{
  167. 		//printf("\n\n enter inside if condi\n\n");
  168. 		
  169. 		dimGrid.x = (Width-1)/32+1;
  170.         	dimGrid.y = (Width-1)/32+1;
  171. 	
  172. 		dimBlock.x = 32;
  173. 	        dimBlock.y = 32;
  174.  
  175.  
  176.  
  177. 	}
  178.  
  179.      
  180.     //Launch the device computation threads!
  181. 	MatrixMulKernel<<<dimGrid,dimBlock>>>(Md,Nd,Pd,Width);
  183.     //Transfer P from device to host
  184.     	cudaMemcpy(P,Pd,size,cudaMemcpyDeviceToHost);
  185.  
  186.     //Free device matrices
  187.     	cudaFree(Md);
  188.     	cudaFree(Nd);
  189.     	cudaFree(Pd);
  190.  
  191. 	int i;
  192.  
  193. 	fprintf(stdout,"\n\n");
  194.  
  195. 	if(Width<11)
  196. 	{
  197.  
  198.  
  199. 		fprintf(stdout,"\n\nMatrix Multiplication, M x Inv(M) :\n\n");
  200. 		for(i = 0; i < Width*Width; i++)
  201.    		{
  202. 			if(P[i])
  203. 				fprintf(stdout,"%10f ",P[i]) ;
  204. 			else
  205. 				fprintf(stdout,"%9f ",P[i]) ;
  206. 	        	
  207.  
  208.  
  209.  
  210. 			if((i+1)%Width==0)
  211. 				fprintf(stdout,"\n");
  212. 		}
  213.   
  214.  
  215. 	}
  216. 	else
  217. 	{
  218. 		FILE *fp;	
  219. 	
  220. 		fp = fopen("Mat_Inv_out", "a");
  221.  
  222.  		if (!fp) 
  223. 		{
  224. 	    		fprintf(stderr, "Failed to open matAdata.\n");
  225. 	    		exit(1);
  226. 	  	}
  227. 		fprintf(fp,"\n\nMatrix Multiplication, M x Inv(M) :\n\n");
  228. 	 	for(i = 0; i < Width*Width; i++)
  229.    		{	if(P[i])
  230. 				fprintf(fp,"%10f ",P[i]) ;
  231. 			else
  232. 				fprintf(fp,"%9f ",P[i]) ;
  233. 			
  234. 		        if((i+1)%Width==0)
  235. 				fprintf(fp,"\n");
  236. 		}
  237.     		fclose(fp);
  238. 	}
  239.  
  240. 	
  241. 	//printf("\n Matrix multiplication completed !!\n\n"); 
  242. 	free(M);
  243. 	free(N);
  244. 	free(P);
  245.  
  246. }
  247.  
  248.  
  249.  
  250. //
  251. //
  252.  
  253.  
  254.  
  255. void fill(float* h, int w)
  256. { 
  257. 	
  258. 	unsigned int i, num;
  259. 	int divide;
  260. 	FILE *f;
  261.  
  262. 	f=fopen("/dev/urandom", "r");
  263. 	if (!f) {
  264.         	fprintf(stderr, "Failed open file\n");
  265.         	exit(1);
  266.     	}
  267. 	for(i=0; i< w*w; i++)
  268. 	{
  269. 		fread(&num, sizeof(unsigned int), 1, f);
  270. 		fread(&divide, sizeof(int), 1, f);
  271. 		h[i] = ((float)num)/((float)divide);
  272. 		//scanf("%f",&h[i]);
  273. 	}
  274. 	fclose(f);
  275. /*
  276. 	unsigned int i;
  277. 	srand((unsigned int)time(NULL));
  278. 	for(i=0; i< w*w; i++)
  279. 	{
  280. 		h[i] = ((float)rand()/(float)(RAND_MAX)) * 99;
  281. 		//scanf("%f",&h[i]);
  282. 	}
  283. 	
  284. */
  285.  
  286. } 
  287.  
  288.  
  289. //
  290. //
  291.  
  292.  
  293. void test_invert(int n )
  294. {
  295.     	
  296. 	//printf("Enter the order of the square matrix :");
  297. 	//scanf("%d",&n);
  298.     	const int mybatch = 1;
  299.  
  300.  
  301. 	//float* mat1[n * n];
  302. 	float mat1_size = sizeof(float) * n * n;
  303.     	float* mat1 = (float*) malloc(mat1_size);
  304.  
  305. 	fill(mat1, n);
  306.  
  307.     	float *result_flat = (float *)malloc(mybatch*n*n*sizeof(float));
  308.     	float **results = (float **)malloc(mybatch*sizeof(float *));
  309.  
  310.     	for (int i = 0; i < mybatch; i++)
  311.       		results[i] = result_flat + (i*n*n);
  312.  
  313.     	float **inputs = (float **)malloc(mybatch*sizeof(float *));
  314.     	
  315. 	//inputs[0]  = zero_pivot;
  316.  
  317. 	inputs[0]  = mat1;
  318.  
  319.  
  320. 	invert(inputs, results, n, mybatch);
  321.  
  322. 	if(n<11)
  323. 	{
  324.  
  325. 		for (int qq = 0; qq < mybatch; qq++)
  326. 		{
  327. 	      		if(mybatch==1)
  328. 				fprintf(stdout, "Input Matrix, M :\n\n");
  329. 			else
  330. 				fprintf(stdout, "Input Matrix %d:\n\n", qq);
  331. 			
  332. 	      		for(int i=0; i<n; i++)
  333. 	      		{
  334. 	        		for(int j=0; j<n; j++)
  335. 				{	
  336. 					if(inputs[qq][i*n+j])
  337. 		            			fprintf(stdout,"%12f ",inputs[qq][i*n+j]);
  338. 					else
  339. 						fprintf(stdout,"%11f ",inputs[qq][i*n+j]);
  340. 				}
  341. 	        			fprintf(stdout,"\n");
  342. 	      		}
  343. 	    	}
  344. 	    	fprintf(stdout,"\n\n");
  345.  
  346.  
  347.  
  348.  
  349. 	    	for (int qq = 0; qq < mybatch; qq++)
  350. 		{
  351.  
  352. 			if(mybatch==1)
  353. 				fprintf(stdout, "Inverse of the Input Matrix, Inv(M):\n\n");
  354. 			else
  355. 				fprintf(stdout, "Inverse Matrix %d:\n\n", qq);
  356. 	      		for(int i=0; i<n; i++)
  357. 	      		{
  358. 	        		for(int j=0; j<n; j++)
  359. 				{
  360. 					if(results[qq][i*n+j])
  361. 		            			fprintf(stdout,"%10f ",results[qq][i*n+j]);
  362. 					else
  363. 		            			fprintf(stdout,"%9f ",results[qq][i*n+j]);
  364. 	        		
  365. 				}
  366. 				fprintf(stdout,"\n");
  367. 	      		}
  368. 	    	}
  369. 	}
  370.  
  371.  
  372. 	else // order of the matrix is more than 10 x 10 then output the results in the file
  373. 	{
  374. 		printf("\nThe order of matrix is too large to display in terminal\n, Please open the file : Mat_Inv_out.txt located in the current folder. To see the output.\n\n");
  375. 		
  376. 		FILE *fp;
  377.  
  378.  
  379.  		fp = fopen("Mat_Inv_out", "w");
  380.  
  381.  		if (!fp) 
  382. 		{
  383.     			fprintf(stderr, "Failed to open Mat_Inv_out.\n");
  384. 		    	exit(1);
  385.   		}
  386.  
  387.  
  388.  
  389. 		for (int qq = 0; qq < mybatch; qq++)
  390. 		{
  391.  
  392. 			if(mybatch==1)
  393. 				fprintf(fp, "Input Matrix , M:\n\n");
  394. 			else
  395. 				fprintf(fp, "Input Matrix %d:\n\n", qq);
  396.  
  397.  
  398. 	      	
  399. 			
  400. 	      		for(int i=0; i<n; i++)
  401.       			{
  402.         			for(int j=0; j<n; j++)
  403. 				{
  404. 					if(inputs[qq][i*n+j])
  405. 		            			fprintf(fp,"%12f ",inputs[qq][i*n+j]);
  406. 					else
  407. 						fprintf(fp,"%11f ",inputs[qq][i*n+j]);
  408. 				}
  409. 		            		
  410.         			fprintf(fp,"\n");
  411. 	      		}
  412.     		}
  413. 	    	fprintf(fp,"\n\n");
  414.  
  415. 		for (int qq = 0; qq < mybatch; qq++)
  416. 		{
  417. 			if(mybatch==1)
  418. 				fprintf(fp, "Inverse of the Input Matrix, Inv(M):\n\n");
  419. 	      		else
  420. 				fprintf(fp, "Inverse %d:\n\n", qq);
  421. 	      		for(int i=0; i<n; i++)
  422. 	      		{
  423. 	        		for(int j=0; j<n; j++)	
  424. 				{
  425. 					if(results[qq][i*n+j])
  426. 		            			fprintf(fp,"%10f ",results[qq][i*n+j]);
  427. 					else
  428. 		            			fprintf(fp,"%9f ",results[qq][i*n+j]);
  429. 	        		
  430. 				}
  431.  
  432. 	        		fprintf(fp,"\n");
  433. 	      		}
  434. 	    	}
  435.  
  436. 		fclose(fp);
  437. 			
  438. 	}// end of if else condition for output
  439.  
  440. 	float *A, *B;
  441.  
  442. 	A=inputs[0];
  443. 	B=results[0];
  444. 	mul(A, B, n );
  445.  
  446. 	//mul(inputs[0][], results[0][], n );
  447.  
  448. }
  449.  
  450.  
  451. //
  452. //
  453.  
  454.  
  455. int main(int argc, char *argv[])
  456. {
  457. 	if(argc!=2)
  458. 	{
  459. 		printf("Usage: %s <matrix_width>\n", argv[0]);
  460. 		return 0;
  461. 	}
  462.  
  463. 	int w;
  464. 	w = atoi( argv[1] );
  465. 	
  466.    	test_invert(w);
  467.     	return 0;
  468. }
  469.  
  470. /*
  471. $ nvcc -arch=sm_20 -o t540 t540.cu -lcublas
  472. $ ./t540
  473. Input 0:
  474. 0.000000        3.000000        4.000000
  475. 1.000000        3.000000        10.000000
  476. 4.000000        9.000000        16.000000
  477. Input 1:
  478. 0.500000        3.000000        4.000000
  479. 1.000000        3.000000        10.000000
  480. 4.000000        9.000000        16.000000
  481. Input 2:
  482. 0.000000        3.000000        4.000000
  483. 1.000000        5.000000        6.000000
  484. 9.000000        8.000000        2.000000
  485. Input 3:
  486. 22.000000       3.000000        4.000000
  487. 1.000000        5.000000        6.000000
  488. 9.000000        8.000000        2.000000
  489. Inverse 0:
  490. -0.700000       -0.200000       0.300000
  491. 0.400000        -0.266667       0.066667
  492. -0.050000       0.200000        -0.050000
  493. Inverse 1:
  494. -1.076923       -0.307692       0.461538
  495. 0.615385        -0.205128       -0.025641
  496. -0.076923       0.192308        -0.038462
  497. Inverse 2:
  498. -4.750000       3.250000        -0.250000
  499. 6.500000        -4.500000       0.500000
  500. -4.625000       3.375000        -0.375000
  501. Inverse 3:
  502. 0.045894        -0.031401       0.002415
  503. -0.062802       -0.009662       0.154589
  504. 0.044686        0.179952        -0.129227
  505. $
  506. $ nvcc -arch=sm_20 -o t540 t540.cu -lcublas
  507. $ ./t540 
  508. Enter the order of the aquare matrix :4
  509. Input 0:
  510. -0.100222 -2.553872 -69.072723 0.016120 
  511. -2.752346 -1.230871 1.997387 0.606710 
  512. -0.029929 -0.583444 2.733107 0.254404 
  513. -1.844285 -0.070541 1.906255 10.758991 
  514. Inverse 0:
  515. 0.017501	-0.374555	0.713068	0.004234	
  516. -0.056876	-0.005437	-1.457745	0.034861	
  517. -0.012399	0.000729	0.052888	-0.001273	
  518. 0.004824	-0.064370	0.103305	0.094125	
  519. */

  1. //mat_inv_userInput.cu
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <cublas_v2.h>
  5. #include<math.h>
  6.  
  7.  
  8. #define cudacall(call)                                                                                                          \
  9.     do                                                                                                                          \
  10.     {                                                                                                                           \
  11.         cudaError_t err = (call);                                                                                               \
  12.         if(cudaSuccess != err)                                                                                                  \
  13.         {                                                                                                                       \
  14.             fprintf(stderr,"CUDA Error:\nFile = %s\nLine = %d\nReason = %s\n", __FILE__, __LINE__, cudaGetErrorString(err));    \
  15.             cudaDeviceReset();                                                                                                  \
  16.             exit(EXIT_FAILURE);                                                                                                 \
  17.         }                                                                                                                       \
  18.     }                                                                                                                           \
  19.     while (0)
  20.  
  21. #define cublascall(call)                                                                                        \
  22.     do                                                                                                          \
  23.     {                                                                                                           \
  24.         cublasStatus_t status = (call);                                                                         \
  25.         if(CUBLAS_STATUS_SUCCESS != status)                                                                     \
  26.         {                                                                                                       \
  27.             fprintf(stderr,"CUBLAS Error:\nFile = %s\nLine = %d\nCode = %d\n", __FILE__, __LINE__, status);     \
  28.             cudaDeviceReset();                                                                                  \
  29.             exit(EXIT_FAILURE);                                                                                 \
  30.         }                                                                                                       \
  31.                                                                                                                 \
  32.     }                                                                                                           \
  33.     while(0)
  34.  
  35.  
  36. //
  37. //
  38.  
  39. void invert(float** src, float** dst, int n, int batchSize)
  40. {
  41.     	cublasHandle_t handle;
  42.     	cublascall(cublasCreate_v2(&handle));
  43.  
  44.     	int *P, *INFO;
  45.  
  46.     	cudacall(cudaMalloc(&P, n * batchSize * sizeof(int)));
  47.     	cudacall(cudaMalloc(&INFO,  batchSize * sizeof(int)));
  48.  
  49.     	int lda = n;
  50.  
  51.     	float **A = (float **)malloc(batchSize*sizeof(float *));
  52.     	float **A_d, *A_dflat;
  53.  
  54.     	cudacall(cudaMalloc(&A_d,batchSize*sizeof(float *)));
  55.     	cudacall(cudaMalloc(&A_dflat, n*n*batchSize*sizeof(float)));
  56.  
  57. 	A[0] = A_dflat;
  58.     	for (int i = 1; i < batchSize; i++)
  59.       		A[i] = A[i-1]+(n*n);
  60.  
  61.     	cudacall(cudaMemcpy(A_d,A,batchSize*sizeof(float *),cudaMemcpyHostToDevice));
  62.    
  63.  	for (int i = 0; i < batchSize; i++)
  64.       		cudacall(cudaMemcpy(A_dflat+(i*n*n), src[i], n*n*sizeof(float), cudaMemcpyHostToDevice));
  65.  
  66.  
  67.     	cublascall(cublasSgetrfBatched(handle,n,A_d,lda,P,INFO,batchSize));
  68.  
  69.  
  70.     	int INFOh[batchSize];
  71.     	cudacall(cudaMemcpy(INFOh,INFO,batchSize*sizeof(int),cudaMemcpyDeviceToHost));
  72.  
  73.     	for (int i = 0; i < batchSize; i++)
  74.       		if(INFOh[i]  != 0)
  75.       		{
  76.         		fprintf(stderr, "Factorization of matrix %d Failed: Matrix may be singular\n", i);
  77.         		cudaDeviceReset();
  78.         		exit(EXIT_FAILURE);
  79.       		}
  80.  
  81.     	float **C = (float **)malloc(batchSize*sizeof(float *));
  82.     	float **C_d, *C_dflat;
  83.  
  84.     	cudacall(cudaMalloc(&C_d,batchSize*sizeof(float *)));
  85.     	cudacall(cudaMalloc(&C_dflat, n*n*batchSize*sizeof(float)));
  86.     	C[0] = C_dflat;
  87.     	for (int i = 1; i < batchSize; i++)
  88.       		C[i] = C[i-1] + (n*n);
  89.     	cudacall(cudaMemcpy(C_d,C,batchSize*sizeof(float *),cudaMemcpyHostToDevice));
  90.     	cublascall(cublasSgetriBatched(handle,n,(const float **)A_d,lda,P,C_d,lda,INFO,batchSize));
  91.  
  92.     	cudacall(cudaMemcpy(INFOh,INFO,batchSize*sizeof(int),cudaMemcpyDeviceToHost));
  93.  
  94.     	for (int i = 0; i < batchSize; i++)
  95. 	      	if(INFOh[i] != 0)
  96.       		{
  97.         		fprintf(stderr, "Inversion of matrix %d Failed: Matrix may be singular\n", i);
  98.         		cudaDeviceReset();
  99.       			exit(EXIT_FAILURE);
  100.       		}
  101.     	for (int i = 0; i < batchSize; i++)
  102.       	cudacall(cudaMemcpy(dst[i], C_dflat + (i*n*n), n*n*sizeof(float), cudaMemcpyDeviceToHost));
  103.     	
  104. 	cudaFree(A_d); cudaFree(A_dflat); free(A);
  105. 	cudaFree(C_d); cudaFree(C_dflat); free(C);
  106.     	cudaFree(P); cudaFree(INFO); cublasDestroy_v2(handle);
  107. }
  108.  
  109.  
  110.  
  111. //
  112. //
  113.  
  114.  
  115.  
  116. __global__ void MatrixMulKernel(float *Md, float *Nd, float *Pd, int Width) 
  117. {
  118.     //2D Thread ID
  119.     int col = blockIdx.x*blockDim.x+threadIdx.x;
  120.     int row = blockIdx.y*blockDim.y+threadIdx.y;
  121.  
  122.     //Pvalue stores the Pd element that is computed by the thread
  123.     float Pvalue = 0;
  124. 	if(col<Width && row < Width)
  125. 	{
  126. 		    for(int k = 0; k < Width ; ++k) 
  127. 		    {
  128. 		        float Mdelement = Md[row*Width + k];
  129. 		        float Ndelement = Nd[k*Width + col];
  130. 		        Pvalue += (Mdelement*Ndelement);
  131. 		
  132. 		    }
  133. 		    Pd[row*Width + col] = Pvalue;
  134. 	}
  135. }
  136.  
  137.  
  138.  
  139. void mul(float* M,float* N,int Width)
  140. {
  141. 		
  142. 	float * P = (float *) malloc(Width*Width*sizeof(float));
  143. 	float *Md, *Nd, *Pd;
  144.  
  145.  
  146.  
  147. 	unsigned long int size = Width*Width*sizeof(float);
  148.   
  149.  
  150.     //Transfer M and N to device memory
  151.     	cudaMalloc((void**)&Md, size);
  152.     	cudaMemcpy(Md,M,size,cudaMemcpyHostToDevice);
  153.  
  154.     	cudaMalloc((void**)&Nd, size);
  155.     	cudaMemcpy(Nd,N,size,cudaMemcpyHostToDevice);
  156.  
  157.     	//Allocate P on the device
  158.     	cudaMalloc((void**)&Pd,size);
  159.  
  160.     	//Setup the execution configuration
  161.     	dim3 dimBlock(Width,Width);
  162.     	dim3 dimGrid(1,1);
  163.  
  164.  
  165. 	if (Width*Width > 1024)
  166. 	{
  167. 		//printf("\n\n enter inside if condi\n\n");
  168. 		
  169. 		dimGrid.x = (Width-1)/32+1;
  170.         	dimGrid.y = (Width-1)/32+1;
  171. 	
  172. 		dimBlock.x = 32;
  173. 	        dimBlock.y = 32;
  174.  
  175.  
  176.  
  177. 	}
  178.  
  179.      
  180.     //Launch the device computation threads!
  181. 	MatrixMulKernel<<<dimGrid,dimBlock>>>(Md,Nd,Pd,Width);
  183.     //Transfer P from device to host
  184.     	cudaMemcpy(P,Pd,size,cudaMemcpyDeviceToHost);
  185.  
  186.     //Free device matrices
  187.     	cudaFree(Md);
  188.     	cudaFree(Nd);
  189.     	cudaFree(Pd);
  190.  
  191. 	int i;
  192.  
  193. 	fprintf(stdout,"\n\n");
  194.  
  195. 	if(Width<11)
  196. 	{
  197.  
  198.  
  199. 		fprintf(stdout,"\n\nMatrix Multiplication, M x Inv(M) :\n\n");
  200. 		for(i = 0; i < Width*Width; i++)
  201.    		{
  202. 			if(P[i])
  203. 				fprintf(stdout,"%10f ",P[i]) ;
  204. 			else
  205. 				fprintf(stdout,"%9f ",P[i]) ;
  206. 	        	
  207.  
  208.  
  209.  
  210. 			if((i+1)%Width==0)
  211. 				fprintf(stdout,"\n");
  212. 		}
  213.   
  214.  
  215. 	}
  216. 	else
  217. 	{
  218. 		FILE *fp;	
  219. 	
  220. 		fp = fopen("Mat_Inv_out", "a");
  221.  
  222.  		if (!fp) 
  223. 		{
  224. 	    		fprintf(stderr, "Failed to open matAdata.\n");
  225. 	    		exit(1);
  226. 	  	}
  227. 		fprintf(fp,"\n\nMatrix Multiplication, M x Inv(M) :\n\n");
  228. 	 	for(i = 0; i < Width*Width; i++)
  229.    		{	if(P[i])
  230. 				fprintf(fp,"%10f ",P[i]) ;
  231. 			else
  232. 				fprintf(fp,"%9f ",P[i]) ;
  233. 			
  234. 		        if((i+1)%Width==0)
  235. 				fprintf(fp,"\n");
  236. 		}
  237.     		fclose(fp);
  238. 	}
  239.  
  240. 	
  241. 	//printf("\n Matrix multiplication completed !!\n\n"); 
  242. 	free(M);
  243. 	free(N);
  244. 	free(P);
  245.  
  246. }
  247.  
  248.  
  249.  
  250. //
  251. //
  252.  
  253.  
  254.  
  255. void fill(float* h, int w)
  256. { 
  257. 	unsigned int i;
  258. 	for(i=0; i< w*w; i++)
  259. 	{
  260. 		scanf("%f",&h[i]);
  261. 	}
  262. 	
  263.  
  264.  
  265. } 
  266.  
  267.  
  268. //
  269. //
  270.  
  271.  
  272. void test_invert( )
  273. {
  274. 	int n;
  275.     	
  276. 	printf("Enter the order of the square matrix :");
  277. 	scanf("%d",&n);
  278.     	const int mybatch = 1;
  279.  
  280.  
  281. 	float mat1_size = sizeof(float) * n * n;
  282.     	float* mat1 = (float*) malloc(mat1_size);
  283.  
  284. 	fill(mat1, n);
  285.  
  286.     	float *result_flat = (float *)malloc(mybatch*n*n*sizeof(float));
  287.     	float **results = (float **)malloc(mybatch*sizeof(float *));
  288.  
  289.     	for (int i = 0; i < mybatch; i++)
  290.       		results[i] = result_flat + (i*n*n);
  291.  
  292.     	float **inputs = (float **)malloc(mybatch*sizeof(float *));
  293.     	
  294. 	//inputs[0]  = zero_pivot;
  295.  
  296. 	inputs[0]  = mat1;
  297.  
  298.  
  299. 	invert(inputs, results, n, mybatch);
  300.  
  301. 	if(n<11)
  302. 	{
  303.  
  304. 		for (int qq = 0; qq < mybatch; qq++)
  305. 		{
  306. 	      		if(mybatch==1)
  307. 				fprintf(stdout, "Input Matrix, M :\n\n");
  308. 			else
  309. 				fprintf(stdout, "Input Matrix %d:\n\n", qq);
  310. 			
  311. 	      		for(int i=0; i<n; i++)
  312. 	      		{
  313. 	        		for(int j=0; j<n; j++)
  314. 				{	
  315. 					if(inputs[qq][i*n+j])
  316. 		            			fprintf(stdout,"%12f ",inputs[qq][i*n+j]);
  317. 					else
  318. 						fprintf(stdout,"%11f ",inputs[qq][i*n+j]);
  319. 				}
  320. 	        			fprintf(stdout,"\n");
  321. 	      		}
  322. 	    	}
  323. 	    	fprintf(stdout,"\n\n");
  324.  
  325.  
  326.  
  327.  
  328. 	    	for (int qq = 0; qq < mybatch; qq++)
  329. 		{
  330.  
  331. 			if(mybatch==1)
  332. 				fprintf(stdout, "Inverse of the Input Matrix, Inv(M):\n\n");
  333. 			else
  334. 				fprintf(stdout, "Inverse Matrix %d:\n\n", qq);
  335. 	      		for(int i=0; i<n; i++)
  336. 	      		{
  337. 	        		for(int j=0; j<n; j++)
  338. 				{
  339. 					if(results[qq][i*n+j])
  340. 		            			fprintf(stdout,"%10f ",results[qq][i*n+j]);
  341. 					else
  342. 		            			fprintf(stdout,"%9f ",results[qq][i*n+j]);
  343. 	        		
  344. 				}
  345. 				fprintf(stdout,"\n");
  346. 	      		}
  347. 	    	}
  348. 	}
  349.  
  350.  
  351. 	else // order of the matrix is more than 10 x 10 then output the results in the file
  352. 	{
  353. 		printf("\nThe order of matrix is too large to display in terminal\n, Please open the file : Mat_Inv_out.txt located in the current folder. To see the output.\n\n");
  354. 		
  355. 		FILE *fp;
  356.  
  357.  
  358.  		fp = fopen("Mat_Inv_out", "w");
  359.  
  360.  		if (!fp) 
  361. 		{
  362.     			fprintf(stderr, "Failed to open Mat_Inv_out.\n");
  363. 		    	exit(1);
  364.   		}
  365.  
  366.  
  367.  
  368. 		for (int qq = 0; qq < mybatch; qq++)
  369. 		{
  370.  
  371. 			if(mybatch==1)
  372. 				fprintf(fp, "Input Matrix , M:\n\n");
  373. 			else
  374. 				fprintf(fp, "Input Matrix %d:\n\n", qq);
  375.  
  376.  
  377. 	      	
  378. 			
  379. 	      		for(int i=0; i<n; i++)
  380.       			{
  381.         			for(int j=0; j<n; j++)
  382. 				{
  383. 					if(inputs[qq][i*n+j])
  384. 		            			fprintf(fp,"%12f ",inputs[qq][i*n+j]);
  385. 					else
  386. 						fprintf(fp,"%11f ",inputs[qq][i*n+j]);
  387. 				}
  388. 		            		
  389.         			fprintf(fp,"\n");
  390. 	      		}
  391.     		}
  392. 	    	fprintf(fp,"\n\n");
  393.  
  394. 		for (int qq = 0; qq < mybatch; qq++)
  395. 		{
  396. 			if(mybatch==1)
  397. 				fprintf(fp, "Inverse of the Input Matrix, Inv(M):\n\n");
  398. 	      		else
  399. 				fprintf(fp, "Inverse %d:\n\n", qq);
  400. 	      		for(int i=0; i<n; i++)
  401. 	      		{
  402. 	        		for(int j=0; j<n; j++)	
  403. 				{
  404. 					if(results[qq][i*n+j])
  405. 		            			fprintf(fp,"%10f ",results[qq][i*n+j]);
  406. 					else
  407. 		            			fprintf(fp,"%9f ",results[qq][i*n+j]);
  408. 	        		
  409. 				}
  410.  
  411. 	        		fprintf(fp,"\n");
  412. 	      		}
  413. 	    	}
  414.  
  415. 		fclose(fp);
  416. 			
  417. 	}// end of if else condition for output
  418.  
  419. float *A, *B;
  421. A=inputs[0];
  422. B=results[0];
  423. mul(A, B, n );
  424.  
  425. //mul(inputs[0][], results[0][], n );
  426.  
  427. }
  428.  
  429. //
  430. //
  431.  
  432.  
  433. int main()
  434. {
  435. 	
  436. 	
  437.    	test_invert();
  438.     	return 0;
  439. }
  440.  
  441. /*
  442. $ nvcc -arch=sm_20 -o t540 t540.cu -lcublas
  443. $ ./t540
  444. Input 0:
  445. 0.000000        3.000000        4.000000
  446. 1.000000        3.000000        10.000000
  447. 4.000000        9.000000        16.000000
  448. Input 1:
  449. 0.500000        3.000000        4.000000
  450. 1.000000        3.000000        10.000000
  451. 4.000000        9.000000        16.000000
  452. Input 2:
  453. 0.000000        3.000000        4.000000
  454. 1.000000        5.000000        6.000000
  455. 9.000000        8.000000        2.000000
  456. Input 3:
  457. 22.000000       3.000000        4.000000
  458. 1.000000        5.000000        6.000000
  459. 9.000000        8.000000        2.000000
  460. Inverse 0:
  461. -0.700000       -0.200000       0.300000
  462. 0.400000        -0.266667       0.066667
  463. -0.050000       0.200000        -0.050000
  464. Inverse 1:
  465. -1.076923       -0.307692       0.461538
  466. 0.615385        -0.205128       -0.025641
  467. -0.076923       0.192308        -0.038462
  468. Inverse 2:
  469. -4.750000       3.250000        -0.250000
  470. 6.500000        -4.500000       0.500000
  471. -4.625000       3.375000        -0.375000
  472. Inverse 3:
  473. 0.045894        -0.031401       0.002415
  474. -0.062802       -0.009662       0.154589
  475. 0.044686        0.179952        -0.129227
  476. $
  477. $ nvcc -arch=sm_20 -o t540 t540.cu -lcublas
  478. $ ./t540 
  479. Enter the order of the aquare matrix :4
  480. Input 0:
  481. -0.100222 -2.553872 -69.072723 0.016120 
  482. -2.752346 -1.230871 1.997387 0.606710 
  483. -0.029929 -0.583444 2.733107 0.254404 
  484. -1.844285 -0.070541 1.906255 10.758991 
  485. Inverse 0:
  486. 0.017501	-0.374555	0.713068	0.004234	
  487. -0.056876	-0.005437	-1.457745	0.034861	
  488. -0.012399	0.000729	0.052888	-0.001273	
  489. 0.004824	-0.064370	0.103305	0.094125	
  490. */

  1. #Makefile
  2. default: all
  3. .PHONY: default
  4. 	
  5. all:
  6. 	nvcc -arch=sm_70 -o mat_inv mat_inv.cu -lcublas -lm
  7. .PHONY: all
  8.  
  9. user: 
  10. 	nvcc -arch=sm_20 -o user_mat_inv mat_inv_userInput.cu -lcublas -lm
  11.  
  12. clean:
  13. 	rm mat_inv user_mat_inv
  14. .PHONY: all

$make

$./mat_inv 3

std out:

  1. Input Matrix, M :
  2.  
  3.    -3.662523    -1.507356     1.094427
  4.    -0.853989     1.768155    45.548477
  5.    -5.994538     6.792342  -323.404816
  6.  
  7.  
  8. Inverse of the Input Matrix, Inv(M):
  9.  
  10.  -0.217016  -0.118222  -0.017385
  11.  -0.135258   0.293317   0.040853
  12.   0.001182   0.008352  -0.001912
  13.  
  14.  
  15.  
  16.  
  17. Matrix Multiplication, M x Inv(M) :
  18.  
  19.   1.000000   0.000000  -0.000000
  20.   0.000000   1.000000  -0.000000
  21.   0.000000  -0.000000   1.000000

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/我家小花儿/article/detail/245510
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号