Demo entry 1068178

LaTeX

   

Submitted by anonymous on Jan 13, 2015 at 21:33
Language: CUDA. Code size: 10.5 kB.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <limits.h>
#include <cuda.h>
#include <curand.h>
#include <curand_kernel.h>

#define INDEX(row,col,m) ((row)*(m) + col - (row)*(row+1)/2)
#define DIM_TRIANG(m) ((m)*(m+1)/2)

const int BlockSize = 16;
const int sizef = sizeof(float);

texture<float> row_k_T;

//legge una matrice m*m da file, verifica che sia simmetrica e ne memorizza solo la sua triangolare superiore
void read_mat_sym(FILE *f, int m, float *T){

	for( int i = 0; i < m; i++ ){
		
		for( int j = 0; j < i; j++ ){
			
			float useless;
			fscanf(f, "%f", &useless);
			if( useless != T[INDEX(j,i,m)] ){
				printf("La matrice data non è simmetrica.\n");
				exit(1);
			}	
		}
		
		for( int j = i; j < m; j++ ){
			
			fscanf(f, "%f", &T[INDEX(i,j,m)] );
		}
	}
}


//stampa una matrice m*m triang. sup. a schermo
void print_mat_tsup(int m, float *T){

	for( int i = 0; i < m; i++ ){
		for( int j = 0; j < i; j++ ){
			printf("0.000000 ");
		}
		for( int j = i; j < m; j++ ){
			printf("%f ",T[INDEX(i,j,m)]);
		}
		printf("\n");
	}
	printf("\n");
}


//stampa una matrice m*m triang. sup. su file
void fprint_mat_tsup(FILE *f, int m, float *T){

	for( int i = 0; i < m; i++ ){
		for( int j = 0; j < i; j++ ){
			fprintf(f,"0.000000 ");
		}
		for( int j = i; j < m; j++ ){
			fprintf(f,"%f ",T[INDEX(i,j,m)]);
		}
		fprintf(f,"\n");
	}
}


//calcola il prodotto della trasposta di una matrice m*m triang. sup. per la matrice stessa (At * A)
__global__ void mat_prod_tsup(int m, float* dev_A, float* dev_C) {
   int Row = blockIdx.y * blockDim.y + threadIdx.y;
   int Col = blockIdx.x * blockDim.x + threadIdx.x;

   float sum = 0.0;
   int numBlocks = gridDim.x;
   for (int n = 0; n < numBlocks; n++) {

	  __shared__ float At_s[BlockSize][BlockSize];
	  __shared__ float A_s[BlockSize][BlockSize];

	  int c = n * BlockSize + threadIdx.x;
	  int r = n * BlockSize + threadIdx.y;

	  At_s[threadIdx.y][threadIdx.x] = (c <= Row) ? (dev_A[INDEX(c,Row,m)]) : 0;
	  A_s[threadIdx.y][threadIdx.x] = (r <= Col) ? (dev_A[INDEX(r,Col,m)]) : 0;
	  __syncthreads();

	  int k_end = (n == numBlocks - 1 ? m - n * BlockSize : BlockSize);
	  for (int k = 0; k < k_end; k++) {
		 sum += At_s[threadIdx.y][k] * A_s[k][threadIdx.x];
	  }
	  __syncthreads();
   }

   if ( Row <= Col && Col < m )
	  dev_C[INDEX(Row,Col,m)] = sum;
}


// calcola la differenza di due matrici m*m triangolari superiori
__global__ void diff_tsup(int m, float *dev_A, float *dev_B, float *dev_C){

	int Row = blockIdx.y * blockDim.y + threadIdx.y;
	int Col = blockIdx.x * blockDim.x + threadIdx.x;

	if( Row <= Col && Col < m )
		dev_C[INDEX(Row,Col,m)] = dev_A[INDEX(Row,Col,m)] - dev_B[INDEX(Row,Col,m)];
}


//calcola la norma infinito di un vettore
float norm_infty(int s, float *v){

	float M = 0, temp = 0;

	for( int i=0; i<s; i++){
		temp = fabsf(v[i]);
		if ( M < temp )
			M = temp;
	}

	return M;
}


// crea matrici m*m sim. def. pos. la cui fatt. di Cholesky è una triang. sup. con tutti "coef" e ne memorizza solo la parte triangolare superiore
__global__ void special_mat(int m, float *dev_A, float coef){

	int Col = blockDim.x * blockIdx.x + threadIdx.x;
	int Row = blockDim.y * blockIdx.y + threadIdx.y;

	if ( Row <= Col && Col < m ){
		dev_A[INDEX(Row,Col,m)] = (Row + 1) * coef * coef;
	}

}


//calcola la norma infinito di A - Rt*R
float err_chol(int m, float *A, float *R){

	float err = 0.0, *delta;
	float *dev_A, *dev_R, *dev_Ac;
	
	int mem_size = DIM_TRIANG(m);
	
	dim3 dimBlock( BlockSize, BlockSize );
	dim3 dimGrid( (m + BlockSize -1)/BlockSize, (m + BlockSize -1)/BlockSize );
	
	delta = (float *)malloc(mem_size * sizef);
	cudaMalloc((void**) &dev_A, mem_size * sizef);
	cudaMalloc((void**) &dev_Ac, mem_size * sizef);
	cudaMalloc((void**) &dev_R, mem_size * sizef);

	cudaMemcpy(dev_A, A, mem_size * sizef, cudaMemcpyHostToDevice);
	cudaMemcpy(dev_R, R, mem_size * sizef, cudaMemcpyHostToDevice);
	
	mat_prod_tsup<<< dimGrid, dimBlock >>>(m, dev_R, dev_Ac);
	diff_tsup<<< dimGrid, dimBlock >>>(m,dev_A,dev_Ac,dev_R);
	
	cudaMemcpy(delta, dev_R, mem_size * sizef, cudaMemcpyDeviceToHost);

	err = norm_infty(mem_size,delta);

	cudaFree(dev_A);
	cudaFree(dev_R);
	cudaFree(dev_Ac);
	
	return err;
}





/***********************************************************/
/****************ALGORITMO CHOLESKY C***********************/
/***********************************************************/

// calcola la fattorizzazione di Cholesky in C partendo dalla triangolare superiore
void chol_c(int m, float *A, float *R) {
	
	int mem_size = DIM_TRIANG(m); 
	memcpy(R,A,mem_size * sizef);

	for(int k=0; k<m; k++){
		for(int j=k+1; j<m; j++){
			for(int i=j; i<m; i++){
				R[INDEX(j,i,m)] = R[INDEX(j,i,m)] -R[INDEX(k,i,m)]*R[INDEX(k,j,m)]/R[INDEX(k,k,m)];
			}
		}

		float temp = R[INDEX(k,k,m)];
		if ( temp <= 0.0 ){
			printf("La matrice non è definita positiva o si è verificato un problema di aritmetica floating point.\n");
			exit(1);
		}

		for(int i=k; i<m; i++){
			R[INDEX(k,i,m)] = R[INDEX(k,i,m)]/sqrt(temp);
		}
	}

}

/************************FINE*******************************/
/****************ALGORITMO CHOLESKY C***********************/
/***********************************************************/





/***********************************************************/
/****************ALGORITMO CHOLESKY CUDA********************/
/***********************************************************/

//compie un'iterazione dell'algoritmo con memoria texture
__global__ void chol_iteration(int m, float *dev_R, int k){

	int Col = blockDim.x * blockIdx.x + threadIdx.x + k;
	int Row = blockDim.y * blockIdx.y + threadIdx.y + k;
	
	if ( Row <= Col && Col < m ){
		if ( Row > k )
			dev_R[INDEX(Row,Col,m)] = dev_R[INDEX(Row,Col,m)] - tex1Dfetch(row_k_T, Col - k) * tex1Dfetch(row_k_T, Row - k) / tex1Dfetch(row_k_T, 0);
		else
			dev_R[INDEX(Row,Col,m)] = dev_R[INDEX(Row,Col,m)] / sqrt(tex1Dfetch(row_k_T, 0));
	}
}

void chol_cuda(int m, float *A, float *R){

	dim3 dimBlock( BlockSize, BlockSize );

	float *dev_R, *row_k, temp = A[0];
	int mem_size = DIM_TRIANG(m);
	
	cudaMalloc((void**) &dev_R, mem_size * sizef);
	cudaMemcpy(dev_R, A, mem_size * sizef, cudaMemcpyHostToDevice);
	
	for (int k=0; k<m; k++){
		if ( temp <= 0.00 ){
			printf("La matrice non è definita positiva o si è verificato un problema di aritmetica floating point.\n");
			exit(1);
		}
	
		cudaMalloc((void**) &row_k, (m-k) * sizef);
		cudaBindTexture(NULL, row_k_T, row_k, (m-k) * sizef);
		
		cudaMemcpy(row_k, &dev_R[INDEX(k,k,m)], (m-k)*sizef, cudaMemcpyDeviceToDevice);
		
		dim3 dG( ((m-k) + BlockSize -1)/BlockSize, ((m-k) + BlockSize -1)/BlockSize );
		chol_iteration<<< dG, dimBlock >>>(m,dev_R,k);
		
		cudaMemcpy(&temp, &dev_R[INDEX(k+1,k+1,m)], sizef, cudaMemcpyDeviceToHost);
		
		cudaUnbindTexture(row_k_T);
	}

	cudaMemcpy(R, dev_R, mem_size * sizef, cudaMemcpyDeviceToHost);
	cudaFree(dev_R);
}

/************************FINE*******************************/
/****************ALGORITMO CHOLESKY CUDA********************/
/***********************************************************/





int main(){

	char file[20];
	char *t="test";

	printf("Per confrontare gli algoritmi scrivi 'test' altrimenti scrivi il nome del file contenente la matrice da fattorizzare\n");
	scanf("%19s",file);
	
	
// Confronto tra algoritmo in c e in cuda su matrici m*m, m=2^n con n massimo richiesto da riga di comando
	if ( strncmp(file, t, 4) == 0 ){
		int n;
		float coeff;

		printf("Fino a che potenza di 2 vuoi arrivare come dimensione delle matrici test?\n");
		scanf("%d",&n);
		printf("Inserire il coefficiente di generazione delle matrici test:\n");
		scanf("%f",&coeff);

		float *dev_A, *A, *R_c, *R_cuda;
		float time_c[n], time_cuda[n];
		float err_c[n], err_cuda[n];

		for( int i = 1; i <=n; i++ ){
			int m = pow(2,i);
			int mem_size = DIM_TRIANG(m);

			cudaMalloc((void**) &dev_A, mem_size * sizef);
			A = (float *)malloc(mem_size*sizef);
			R_c = (float *)malloc(mem_size*sizef);
			R_cuda= (float *)malloc(mem_size*sizef);
			
			dim3 dimBlock( BlockSize, BlockSize );
			dim3 dimGrid( (m + BlockSize -1)/BlockSize, (m + BlockSize -1)/BlockSize );	
			special_mat<<< dimGrid,dimBlock >>>(m,dev_A,coeff);
			cudaMemcpy(A,dev_A,mem_size*sizef,cudaMemcpyDeviceToHost);
			cudaFree(dev_A);

			clock_t start, stop;

			start = clock();
			chol_c(m,A,R_c);
			stop = clock();
			time_c[i-1] = (stop-start)/(float) CLOCKS_PER_SEC;

			start = clock();
			chol_cuda(m,A,R_cuda);
			stop = clock();
			time_cuda[i-1] = (stop-start)/(float) CLOCKS_PER_SEC;
	 
			printf("\nCon matrice %d x %d c impiega %f secondi mentre cuda impiega %f secondi\n", m, m, time_c[i-1], time_cuda[i-1]);

			//controllo della precisione
			err_c[i-1] = err_chol(m,A,R_c);
			err_cuda[i-1] = err_chol(m,A,R_cuda);
			
			printf("L'errore commesso da c è %f, quello commesso da cuda %f\n", err_c[i-1], err_cuda[i-1]);
		}
	}
	
// calcola la fattorizzazione di Cholesky di una matrice letta da file (nella prima riga del file dev'esserci la dimensione della matrice).
// Richiede da riga di comando se si vuole utilizzare l'algoritmo cuda o c
	else{
	
		FILE *f, *g;
		int m; // dimensione matrice
		
		f=fopen(file,"r");
		
		if( f == NULL ){
			printf("Non esiste un file con questo nome.\n");
			exit(1);
		}
		fscanf(f,"%d", &m);

		int mem_size = DIM_TRIANG(m);
		float *A, *R, err = 0.0;
		float time = 0.0;
		clock_t start, stop;

		A = (float *)calloc(mem_size, sizef);
		R = (float *)calloc(mem_size, sizef);

		read_mat_sym(f, m, A);

		char opt[6];
		char *type = "cuda";

		printf("Per usare l'algoritmo cuda scrivi 'cuda' altrimenti verrà utilizzato l'algoritmo c.\n");
		scanf("%5s",opt);

		if ( strncmp(opt, type, 4) == 0 ){
			start = clock();
			chol_cuda(m,A,R);
			stop = clock();
			time = (stop-start)/(float) CLOCKS_PER_SEC;
			err = err_chol(m,A,R);
		}
		else{
			start = clock();
			chol_c(m,A,R);
			stop = clock();
			time = (stop-start)/(float) CLOCKS_PER_SEC;
			err = err_chol(m,A,R);
		}

		printf("L'algoritmo scelto ha impiegato %f secondi con un errore di %f\n", time, err);

		char where[6];
		char *F = "file";

		printf("Per l'output sul file chol.dat scrivi 'file' altrimenti la matrice verra' stampata a schermo.\n");
		scanf("%5s",where);
		
		if ( strncmp(F, where, 4) == 0 ){
			g=fopen("chol.dat","w");
			fprint_mat_tsup(g,m,R);
		}
		else{
			print_mat_tsup(m,R);
		}
	}

	return 0;

}

This snippet took 0.03 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).