codice

/**************************************************************************************
 *
 * CdL Magistrale in Ingegneria Informatica
 * Corso di Architetture e Programmazione dei Sistemi di Elaborazione - a.a. 2017/18
 *
 * Progetto dell'algoritmo di PageRank
 * in linguaggio assembly x86-32 + SSE
 *
 * Fabrizio Angiulli, 23 novembre 2017
 *
 **************************************************************************************/

/*

 Software necessario per l'esecuzione:

     NASM (www.nasm.us)
     GCC (gcc.gnu.org)

 entrambi sono disponibili come pacchetti software
 installabili mediante il packaging tool del sistema
 operativo; per esempio, su Ubuntu, mediante i comandi:

     sudo apt-get install nasm
     sudo apt-get install gcc

 potrebbe essere necessario installare le seguenti librerie:

     sudo apt-get install lib32gcc-4.8-dev (o altra versione)
     sudo apt-get install libc6-dev-i386

 Per generare il file eseguibile:

 nasm -f elf32 pagerank32.nasm && gcc -O0 -m32 -msse pagerank32.o pagerank32c.c -o pagerank32c -lm && ./pagerank32c

 oppure

 ./runpagerank32

*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <xmmintrin.h>


#define	MATRIX		double*
#define	VECTOR		double*
#define	GRAPH		double*


typedef struct {
	int x;
	int y;
} location;


typedef struct {
	char* file_name;
	MATRIX P; // codifica dense
	GRAPH G; // codifica full
	int N; // numero di nodi
	int M; // numero di archi
	double c; // default=0.85
	double eps; // default=1e-5
	int format; // 0=sparse, 1=full
	int prec; // 0=single, 1=double
	int opt; // 0=nopt, 1=opt
	double* pagerank;
	int silent;
	int display;
} params;


/*
 *
 *	Le funzioni sono state scritte assumento che le matrici siano memorizzate
 * 	mediante un array (float*), in modo da occupare un unico blocco
 * 	di memoria, ma a scelta del candidato possono essere
 * 	memorizzate mediante array di array (float**).
 *
 * 	L'assunzione corrente è che le matrici siano in row-major order.
 *
 */


void* get_block(int size, int elements) {
	return _mm_malloc(elements*size,16);
}


void free_block(void* p) {
	_mm_free(p);
}


MATRIX alloc_matrix(int rows, int cols) {
	return (MATRIX) get_block(sizeof(double),rows*cols);
}


void dealloc_matrix(MATRIX mat) {
	free_block(mat);
}


/*
 *
 * 	load_dense
 * 	===========
 *
 *	Legge da file una matrice di N righe
 * 	e M colonne e la memorizza in un array lineare in row-major order
 *
 * 	Codifica del file:
 * 	primi 4 byte: numero di righe (N) --> numero intero a 32 bit
 * 	successivi 4 byte: numero di colonne (M) --> numero intero a 32 bit
 * 	successivi N*M*4 byte: matrix data in row-major order --> numeri floating-point a precisione doppia
 *
 */
MATRIX load_dense(char* filename, int *n, int *m) {
	FILE* fp;
	int rows, cols, status, i;
	char fpath[256];

	sprintf(fpath, "%s.matrix", filename);
	fp = fopen(fpath, "rb");

	if (fp == NULL) {
		printf("'%s' : bad matrix file name!\n", fpath);
		exit(0);
	}

	status = fread(&rows, sizeof(int), 1, fp);
	status = fread(&cols, sizeof(int), 1, fp);

	MATRIX data = alloc_matrix(rows,cols);
	status = fread(data, sizeof(double), rows*cols, fp);
	fclose(fp);

	*n = rows;
	*m = rows*cols;

	return data;
}

/*
 *
 * 	load_sparse
 * 	===========
 *
 *	Legge da file il grafo
 *
 * 	Codifica del file:
 * 	primi 4 byte: numero di nodi (N) --> numero intero
 * 	successivi 4 byte: numero di archi (M) --> numero intero
 * 	successivi M*2*4 byte: M archi rappresentati come coppie (i,j) di interi a 32 bit
 *
 */
GRAPH load_sparse(char* filename, int *n, int *m) {
	FILE* fp;
	int nodes, arcs, status, i, j, k;
	char fpath[256];
	int sorg, dest;

	sprintf(fpath, "%s.graph", filename);
	fp = fopen(fpath, "rb");

	if (fp == NULL) {
		printf("'%s' : bad graph file name!\n", fpath);
		exit(0);
	}

	status = fread(&nodes, sizeof(int), 1, fp);
	status = fread(&arcs, sizeof(int), 1, fp);

	//printf("nodi: %d", nodes);
	//printf("archi: %d",arcs);

	GRAPH g = alloc_matrix(nodes,nodes);

	for ( j=0;j<nodes;j++)
		for(k=0;k<nodes;k++)
			g[i*nodes+j]=0.0;

	for (i = 0; i < arcs; i++) {
		status = fread(&sorg, sizeof(int), 1, fp);
		status = fread(&dest, sizeof(int), 1, fp);
		g[sorg*nodes+dest]=1.0;
	}
	fclose(fp);

	//for ( j=0;j<nodes;j++)
			//printf("colonna: %d \n elemento:%f",j, g[j]);

	*n = nodes;
	*m = arcs;

	return g;

}


void save_pageranks(char* filename, int n, VECTOR pagerank) {
	FILE* fp;
	int i;
	char fpath[256];

	sprintf(fpath, "%s_pageranks.txt", filename);
	fp = fopen(fpath, "w");
	for (i = 0; i < n; i++)
		fprintf(fp, "%.14g\n", pagerank[i]);
	fclose(fp);
}

double* outdegree(GRAPH m, int n){
	//serve solo per il caso sparso ecco perchè Graph
	int i,j;
	double outdegree;
	double* v= (double*) get_block(sizeof(double),n);
	if(v==NULL){
		printf("Memoria non disponibile!");
		exit(-1);
	}
	double nodi=0;
	for(i=0;i<n;i++){
		outdegree=0;
		for(j=0;j<n;j++){
			//printf("riga: %d \n colonna:%d \n grado:%f",i,j,outdegree);
			outdegree=outdegree+m[i*n+j];
		}
		v[i]=outdegree;
		nodi+=outdegree;
		//printf("%f",outdegree);
	}
	//printf("%f",nodi);
	return v;
}

void matriceProb(MATRIX p, int n, double* v){
	int i,j;
	//Calcolo D
	MATRIX d=alloc_matrix(n, n);
	if(d==NULL){
		printf("Memoria non disponibile!");
		dealloc_matrix(p);
		exit(-1);
	}
	double k= 1.0/n;
	for( i=0;i<n;i++){
		for( j=0;j<n;j++){
			d[i*n+j]=v[i]*k;
			//printf("%g - ",d[i*n+j]);
		}
	}

	for( i=0;i<n;i++){
		for( j=0;j<n;j++){
			p[i*n+j]= p[i*n+j]+d[i*n+j];
			//printf("%g - ",p[i*n+j]);
		}
	}

	dealloc_matrix(d);

}

MATRIX matriceTransizione(GRAPH m, int n, double c){
	int i,j,k;
	double* d=outdegree(m,n);
	MATRIX p=alloc_matrix(n, n);
	if(p==NULL){
		printf("Memoria non disponibile!");
		exit(-1);
	}

	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			//printf("- %g- ",m[i*n+j]);
			if(d[i]!=0.0) p[i*n+j]=m[i*n+j]/d[i];
			else p[i*n+j]=0.0;
			//printf("%g - ",p[i*n+j]);
		}
	}

	//Calcolo il vettore sigma
	double* v= (double*) get_block(sizeof(double),n);
	if(v==NULL){
		printf("Memoria non disponibile!");
		dealloc_matrix(p);
		exit(-1);
	}
	for (k=0;k<n;k++){
		if(d[k]== 0.0)
			v[k]=1.0;
		else v[k]=0.0;
	}

	free_block(d);
	matriceProb(p,n,v);
	free_block(v);

	double cost=(1.0-c)/n;
	printf("%g - ",cost);
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			p[i*n+j]= p[i*n+j]*c+cost;
			//printf("%g - ",p[i*n+j]);
		}
	}
	return p;
}

VECTOR autovettore_single(VECTOR pi,MATRIX p,int n){
	int i,j;
	float s;
	VECTOR r=(VECTOR) get_block(sizeof(float),n);
	if(r == NULL){
		printf("Memoria non disponibile!");
		free(pi);
		exit(-1);
	}

	for(i=0;i<n;i++){
		s=0;
		for(j=0;j<n;j++){
		//Si accede a p come se fosse CMO dato che serve la 		//trasposta
			s+=p[i+j*n]*pi[j];
		}
		r[i]=s;
	}
	return r;
}

VECTOR autovettore_double(VECTOR pi,MATRIX p,int n){

	int i,j;
	double s;
	VECTOR r=(VECTOR) get_block(sizeof(double),n);
	if(r == NULL){
		printf("Memoria non disponibile!");
		free(pi);
		exit(-1);
	}

	for(i=0;i<n;i++){
		s=0;
		for(j=0;j<n;j++){
		//Si accede a p come se fosse CMO dato che serve la 		//trasposta
			s+=p[i+j*n]*pi[j];
		}
		r[i]=s;
	}
	return r;
}

double terminationError(VECTOR i,VECTOR j, int n){
	int k;
	double norma=0.0;
	for(k=0;k<n;k++){
		double diff=i[k]-j[k];
		norma+= diff*diff;
	}
	return sqrt(norma);
}

void powerIteration(VECTOR pi,MATRIX p,int n,double eps,int prec){
	VECTOR v;
	if(prec==0) v=autovettore_single(pi,p,n);
	else autovettore_double(pi,p,n);
	double delta=terminationError(v,pi,n);
	void* temp;

	while(delta>eps){
		temp=v;
		if(prec==0) v=autovettore_single(pi,p,n);
		else autovettore_double(pi,p,n);
		pi=temp;
		delta=terminationError(v,pi,n);
	}

	pi=v;
}

extern void pagerank32(params* input); //funzione assembly

void pagerank(params* input) {

	if (input->format == 0) input->P = matriceTransizione(input->G, input->N , input->c);

	if(input->opt==1)pagerank32(input); // Chiamata alla funzione assembly
	else{
		int i;
		VECTOR v;
		//VECTOR è un double*...è corretto così?
		if (input->prec== 0)
			v=(VECTOR) get_block(sizeof(float),input->N);
		else 	v=(VECTOR) get_block(sizeof(double),input->N);

		if(v==NULL){
			printf("Memoria non disponibile!");
			free(input); //Il main lo rialloca??
			exit(-1);
		}
		//inizializzazione di v
		for(i=0; i<input->N; i++)
			v[i]=1/input->N;
		powerIteration(v,input->P,input->N,input->eps,input->prec);
		input->pagerank=v;
	}
}


#define	SPARSE	0
#define	DENSE	1
#define	SINGLE	0
#define	DOUBLE	1
#define	NOPT	0
#define	OPT		1



int main(int argc, char** argv) {

	params* input = malloc(sizeof(params));

	input->file_name = NULL;
	input->P = NULL; // dense format
	input->G = NULL; // sparse format
	input->N = 0; // number of nodes
	input->M = 0; // number of arcs
	input->c = 0.85;
	input->eps = 1e-5;
	input->format = SPARSE; // 0=sparse, 1=dense
	input->prec = SINGLE; // 0=single, 1=double
	input->opt = NOPT; // 0=nopt, 1=opt
	input->silent = 0; // 0=false,<>0=true
	input->display = 0; // 0=false <>0=true

	int i, j;

	int par = 1;
	while (par < argc) {
		if (par == 1) {
			input->file_name = argv[par];
			par++;
		} else if (strcmp(argv[par],"-s") == 0) {
			input->silent = 1;
			par++;
		} else if (strcmp(argv[par],"-d") == 0) {
			input->display = 1;
			par++;
		} else if (strcmp(argv[par],"-c") == 0) {
			par++;
			if (par >= argc) {
				printf("Missing c value!\n");
				exit(1);
			}
			input->c = atof(argv[par]);
			par++;
		} else if (strcmp(argv[par],"-eps") == 0) {
			par++;
			if (par >= argc) {
				printf("Missing eps value!\n");
				exit(1);
			}
			input->eps = atof(argv[par]);
			par++;
		} else if (strcmp(argv[par],"-sparse") == 0) {
			input->format = SPARSE;
			par++;
		} else if (strcmp(argv[par],"-dense") == 0) {
			input->format = DENSE;
			par++;
		} else if (strcmp(argv[par],"-single") == 0) {
			input->prec = SINGLE;
			par++;
		} else if (strcmp(argv[par],"-double") == 0) {
			input->prec = DOUBLE;
			par++;
		} else if (strcmp(argv[par],"-nopt") == 0) {
			input->opt = NOPT;
			par++;
		} else if (strcmp(argv[par],"-opt") == 0) {
			input->opt = OPT;
			par++;
		} else
			par++;
	}

	if (!input->silent) {
		printf("Usage: %s <input_file_name> [-d][-s][-sparse|-dense][-single|-double][-nopt|-opt][-c <value>][-eps <value>]\n", argv[0]);
		printf("\nParameters:\n");
		printf("\t-d : display input and output\n");
		printf("\t-s : silent\n");
		printf("\t-sparse/-full: input format (sparse=list of arcs,full=matrix)\n");
		printf("\t-single/-double: floating-point precision (only sparse format)\n");
		printf("\t-nopt/-opt: disable/enable optimizations\n");
		printf("\t-c <value> : 1-teleportation_probability (default 0.85)\n");
		printf("\t-eps <value> : termination error (default 1e-5)\n");
		printf("\n");
	}

	if (input->file_name == NULL || strlen(input->file_name) == 0) {
		printf("Missing input file name!\n");
		exit(1);
	}

	if (input->format == 0)
		input->G = load_sparse(input->file_name, &input->N, &input->M);
	else
		input->P = load_dense(input->file_name, &input->N, &input->M);

	if (!input->silent) {
		printf("Input file name: '%s'\n", input->file_name);
		printf("Number of nodes: %d\n", input->N);
		printf("Number of arcs: %d\n", input->M);
		printf("Parameter c: %f\n", input->c);
		printf("Parameter eps: %f\n", input->eps);
	}

	clock_t t = clock();
	pagerank(input);
	t = clock() - t;

	if (!input->silent)
		printf("\nExecution time = %.3f seconds\n", ((float)t)/CLOCKS_PER_SEC);
	else
		printf("%.3f\n", ((float)t)/CLOCKS_PER_SEC);

	if (input->pagerank != NULL)
	{
		if (!input->silent && input->display) {
			printf("\nPageRanks:\n");
			for (i = 0; i < input->N; i++) {
				printf("%d %.14g\n", i+1, input->pagerank[i]);
			}
		}
		save_pageranks(input->file_name, input->N, input->pagerank);
	}

	return 0;
}