starpu-max/julia/tst/mult.c

/* StarPU --- Runtime system for heterogeneous multicore architectures.
 *
 * Copyright (C) 2018                                     Alexis Juven
 * Copyright (C) 2012,2013                                Inria
 * Copyright (C) 2009-2011,2013-2015                      Université de Bordeaux
 * Copyright (C) 2010                                     Mehdi Juhoor
 * Copyright (C) 2010-2013,2015,2017                      CNRS
 *
 * StarPU is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or (at
 * your option) any later version.
 *
 * StarPU is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *
 * See the GNU Lesser General Public License in COPYING.LGPL for more details.
 */

/*
 * This example shows a simple implementation of a blocked matrix
 * multiplication. Note that this is NOT intended to be an efficient
 * implementation of sgemm! In this example, we show:
 *  - how to declare dense matrices (starpu_matrix_data_register)
 *  - how to manipulate matrices within codelets (eg. descr[0].blas.ld)
 *  - how to use filters to partition the matrices into blocks
 *    (starpu_data_partition and starpu_data_map_filters)
 *  - how to unpartition data (starpu_data_unpartition) and how to stop
 *    monitoring data (starpu_data_unregister)
 *  - how to manipulate subsets of data (starpu_data_get_sub_data)
 *  - how to construct an autocalibrated performance model (starpu_perfmodel)
 *  - how to submit asynchronous tasks
 */

#include <string.h>
#include <math.h>
#include <sys/types.h>
#include <signal.h>

#include <starpu.h>



/*
 * That program should compute C = A * B
 *
 *   A of size (z,y)
 *   B of size (x,z)
 *   C of size (x,y)

              |---------------|
            z |       B       |
              |---------------|
       z              x
     |----|   |---------------|
     |    |   |               |
     |    |   |               |
     | A  | y |       C       |
     |    |   |               |
     |    |   |               |
     |----|   |---------------|

 */





void gpu_mult(void **, void *);
void cpu_mult(void **, void *);





static struct starpu_perfmodel model =
{
		.type = STARPU_HISTORY_BASED,
		.symbol = "history_perf"
};

static struct starpu_codelet cl =
{
		.cpu_funcs = {cpu_mult},
		.cpu_funcs_name = {"cpu_mult"},
		.cuda_funcs = {gpu_mult},
		.nbuffers = 3,
		.modes = {STARPU_R, STARPU_R, STARPU_W},
		.model = &model
};






void multiply_with_starpu(float *A, float *B, float *C,  unsigned xdim,  unsigned ydim,  unsigned zdim, unsigned nslicesx, unsigned nslicesy)
{
	starpu_data_handle_t A_handle, B_handle, C_handle;


	starpu_matrix_data_register(&A_handle, STARPU_MAIN_RAM, (uintptr_t)A,
			ydim, ydim, zdim, sizeof(float));
	starpu_matrix_data_register(&B_handle, STARPU_MAIN_RAM, (uintptr_t)B,
			zdim, zdim, xdim, sizeof(float));
	starpu_matrix_data_register(&C_handle, STARPU_MAIN_RAM, (uintptr_t)C,
			ydim, ydim, xdim, sizeof(float));


	struct starpu_data_filter vert =
	{
			.filter_func = starpu_matrix_filter_vertical_block,
			.nchildren = nslicesx
	};

	struct starpu_data_filter horiz =
	{
			.filter_func = starpu_matrix_filter_block,
			.nchildren = nslicesy
	};


	starpu_data_partition(B_handle, &vert);
	starpu_data_partition(A_handle, &horiz);
	starpu_data_map_filters(C_handle, 2, &vert, &horiz);

	unsigned taskx, tasky;

	for (taskx = 0; taskx < nslicesx; taskx++){
		for (tasky = 0; tasky < nslicesy; tasky++){

			struct starpu_task *task = starpu_task_create();

			task->cl = &cl;
			task->handles[0] = starpu_data_get_sub_data(A_handle, 1, tasky);
			task->handles[1] = starpu_data_get_sub_data(B_handle, 1, taskx);
			task->handles[2] = starpu_data_get_sub_data(C_handle, 2, taskx, tasky);

			starpu_task_submit(task);

		}
	}

	starpu_task_wait_for_all();


	starpu_data_unpartition(A_handle, STARPU_MAIN_RAM);
	starpu_data_unpartition(B_handle, STARPU_MAIN_RAM);
	starpu_data_unpartition(C_handle, STARPU_MAIN_RAM);

	starpu_data_unregister(A_handle);
	starpu_data_unregister(B_handle);
	starpu_data_unregister(C_handle);

}



void init_rand(float * m, unsigned width, unsigned height)
{
	unsigned i,j;

	for (j = 0 ; j < height ; j++){
		for (i = 0 ; i < width ; i++){
			m[j+i*height] = (float)(starpu_drand48());
		}
	}
}


void init_zero(float * m, unsigned width, unsigned height)
{
	memset(m, 0, sizeof(float) * width * height);
}



void sort(unsigned int size, double t[])
{
	unsigned int j;

	int is_sort = 0;

	while(!is_sort){

		is_sort = 1;

		for (j = 0 ; j < size - 1 ; j++){

			if (t[j] > t[j+1]){
				double tmp = t[j];
				t[j] = t[j+1];
				t[j+1] = tmp;
				is_sort = 0;
			}
		}
	}


}


double median_time(unsigned nb_test, unsigned xdim, unsigned ydim, unsigned zdim, unsigned nsclicesx, unsigned nsclicesy)
{
	unsigned i;

	float * A = (float *) malloc(zdim*ydim*sizeof(float));
	float * B = (float *) malloc(xdim*zdim*sizeof(float));
	float * C = (float *) malloc(xdim*ydim*sizeof(float));

	double exec_times[nb_test];

	for (i = 0 ; i < nb_test ; i++){

		double start, stop, exec_t;

		init_rand(A, zdim, ydim);
		init_rand(B, xdim, zdim);
		init_zero(C, xdim, ydim);

		start = starpu_timing_now();
		multiply_with_starpu(A, B, C, xdim, ydim, zdim, nsclicesx, nsclicesy);
		stop = starpu_timing_now();

		exec_t = (stop - start)/1.e6;
		exec_times[i] = exec_t;
	}

	sort(nb_test, exec_times);

	free(A);
	free(B);
	free(C);

	return exec_times[nb_test/2];
}


void display_times(unsigned start_dim, unsigned step_dim, unsigned stop_dim, unsigned nb_tests, unsigned nsclicesx, unsigned nsclicesy)
{
	unsigned dim;

	for (dim = start_dim ; dim <= stop_dim ; dim += step_dim){
		double t = median_time(nb_tests, dim, dim, dim, nsclicesx, nsclicesy);
		printf("%u ; %f\n", dim, t);
	}

}


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

	if (argc != 7){
		printf("Usage : %s start_dim step_dim stop_dim nb_tests nsclicesx nsclicesy\n", argv[0]);
		return 1;
	}


	if (starpu_init(NULL) != EXIT_SUCCESS){
		fprintf(stderr, "ERROR\n");
		return 77;
	}

	unsigned start_dim = (unsigned) atoi(argv[1]);
	unsigned step_dim = (unsigned) atoi(argv[2]);
	unsigned stop_dim = (unsigned) atoi(argv[3]);
	unsigned nb_tests = (unsigned) atoi(argv[4]);
	unsigned nsclicesx = (unsigned) atoi(argv[5]);
	unsigned nsclicesy = (unsigned) atoi(argv[6]);

	display_times(start_dim, step_dim, stop_dim, nb_tests, nsclicesx, nsclicesy);

	starpu_shutdown();

	return 0;
}