rename test cholesky with 2 ctxs

examples/cholesky_2ctxs/cholesky/cholesky.h

@@ -0,0 +1,154 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010, 2011  Université de Bordeaux 1
+ * Copyright (C) 2010  Centre National de la Recherche Scientifique
+ *
+ * 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.
+ */
+
+#ifndef __DW_CHOLESKY_H__
+#define __DW_CHOLESKY_H__
+
+#include <limits.h>
+#include <string.h>
+#include <math.h>
+#include <sys/time.h>
+#ifdef STARPU_USE_CUDA
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <cublas.h>
+#endif
+
+#include <common/blas.h>
+#include <starpu.h>
+
+#define NMAXBLOCKS	32
+
+#define TAG11(k)	((starpu_tag_t)( (1ULL<<60) | (unsigned long long)(k)))
+#define TAG21(k,j)	((starpu_tag_t)(((3ULL<<60) | (((unsigned long long)(k))<<32)	\
+					| (unsigned long long)(j))))
+#define TAG22(k,i,j)	((starpu_tag_t)(((4ULL<<60) | ((unsigned long long)(k)<<32) 	\
+					| ((unsigned long long)(i)<<16)	\
+					| (unsigned long long)(j))))
+
+#define TAG11_AUX(k, prefix)	((starpu_tag_t)( (((unsigned long long)(prefix))<<60)  |  (1ULL<<56) | (unsigned long long)(k)))
+#define TAG21_AUX(k,j, prefix)	((starpu_tag_t)( (((unsigned long long)(prefix))<<60)  			\
+					|  ((3ULL<<56) | (((unsigned long long)(k))<<32)	\
+					| (unsigned long long)(j))))
+#define TAG22_AUX(k,i,j, prefix)    ((starpu_tag_t)(  (((unsigned long long)(prefix))<<60)	\
+					|  ((4ULL<<56) | ((unsigned long long)(k)<<32)  	\
+					| ((unsigned long long)(i)<<16) 			\
+					| (unsigned long long)(j))))
+
+#define BLOCKSIZE	(size/nblocks)
+
+#define BLAS3_FLOP(n1,n2,n3)    \
+        (2*((uint64_t)n1)*((uint64_t)n2)*((uint64_t)n3))
+
+//static unsigned size = 4*1024;
+//static unsigned nblocks = 16;
+static unsigned nbigblocks = 8;
+static unsigned pinned = 0;
+static unsigned noprio = 0;
+static unsigned check = 0;
+
+void chol_cpu_codelet_update_u11(void **, void *);
+void chol_cpu_codelet_update_u21(void **, void *);
+void chol_cpu_codelet_update_u22(void **, void *);
+
+#ifdef STARPU_USE_CUDA
+void chol_cublas_codelet_update_u11(void *descr[], void *_args);
+void chol_cublas_codelet_update_u21(void *descr[], void *_args);
+void chol_cublas_codelet_update_u22(void *descr[], void *_args);
+#endif
+
+double run_cholesky_implicit(unsigned sched_ctx, int start, int argc, char **argv, double *timing, pthread_barrier_t *barrier);
+
+extern struct starpu_perfmodel_t chol_model_11;
+extern struct starpu_perfmodel_t chol_model_21;
+extern struct starpu_perfmodel_t chol_model_22;
+
+static void __attribute__((unused)) parse_args(int argc, char **argv, unsigned *size, unsigned *nblocks)
+{
+	int i;
+	for (i = 1; i < argc; i++) {
+		if (strcmp(argv[i], "-size") == 0) {
+		        char *argptr;
+			(*size) = strtol(argv[++i], &argptr, 10);
+		}
+
+		if (strcmp(argv[i], "-nblocks") == 0) {
+		        char *argptr;
+			(*nblocks) = strtol(argv[++i], &argptr, 10);
+		}
+
+		if (strcmp(argv[i], "-nbigblocks") == 0) {
+		        char *argptr;
+			nbigblocks = strtol(argv[++i], &argptr, 10);
+		}
+
+		if (strcmp(argv[i], "-pin") == 0) {
+			pinned = 1;
+		}
+
+		if (strcmp(argv[i], "-no-prio") == 0) {
+			noprio = 1;
+		}
+
+		if (strcmp(argv[i], "-check") == 0) {
+			check = 1;
+		}
+
+		if (strcmp(argv[i], "-h") == 0) {
+			printf("usage : %s [-pin] [-size size] [-nblocks nblocks] [-check]\n", argv[0]);
+		}
+	}
+}
+
+static void __attribute__((unused)) parse_args_ctx(int start, int argc, char **argv, unsigned *size, unsigned *nblocks)
+{
+	int i;
+	for (i = start; i < argc; i++) {
+		if (strcmp(argv[i], "-size") == 0) {
+		        char *argptr;
+			(*size) = strtol(argv[++i], &argptr, 10);
+		}
+
+		if (strcmp(argv[i], "-nblocks") == 0) {
+		        char *argptr;
+			(*nblocks) = strtol(argv[++i], &argptr, 10);
+		}
+
+		if (strcmp(argv[i], "-nbigblocks") == 0) {
+		        char *argptr;
+			nbigblocks = strtol(argv[++i], &argptr, 10);
+		}
+
+		if (strcmp(argv[i], "-pin") == 0) {
+			pinned = 1;
+		}
+
+		if (strcmp(argv[i], "-no-prio") == 0) {
+			noprio = 1;
+		}
+
+		if (strcmp(argv[i], "-check") == 0) {
+			check = 1;
+		}
+
+		if (strcmp(argv[i], "-h") == 0) {
+			printf("usage : %s [-pin] [-size size] [-nblocks nblocks] [-check]\n", argv[0]);
+		}
+	}
+}
+
+#endif // __DW_CHOLESKY_H__

examples/cholesky_2ctxs/cholesky/cholesky_grain_tag.c

@@ -0,0 +1,382 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010, 2011  Université de Bordeaux 1
+ * Copyright (C) 2010  Mehdi Juhoor <mjuhoor@gmail.com>
+ * Copyright (C) 2010  Centre National de la Recherche Scientifique
+ *
+ * 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.
+ */
+
+#include "cholesky.h"
+
+/*
+ *	Some useful functions
+ */
+
+static struct starpu_task *create_task(starpu_tag_t id)
+{
+	struct starpu_task *task = starpu_task_create();
+		task->cl_arg = NULL;
+		task->use_tag = 1;
+		task->tag_id = id;
+
+	return task;
+}
+
+/*
+ *	Create the codelets
+ */
+
+static starpu_codelet cl11 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = chol_cpu_codelet_update_u11,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u11,
+#endif
+	.nbuffers = 1,
+	.model = &chol_model_11
+};
+
+static struct starpu_task * create_task_11(starpu_data_handle dataA, unsigned k, unsigned reclevel)
+{
+//	printf("task 11 k = %d TAG = %llx\n", k, (TAG11(k)));
+
+	struct starpu_task *task = create_task(TAG11_AUX(k, reclevel));
+	
+	task->cl = &cl11;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = starpu_data_get_sub_data(dataA, 2, k, k);
+	task->buffers[0].mode = STARPU_RW;
+
+	/* this is an important task */
+	task->priority = STARPU_MAX_PRIO;
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG11_AUX(k, reclevel), 1, TAG22_AUX(k-1, k, k, reclevel));
+	}
+
+	return task;
+}
+
+static starpu_codelet cl21 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = chol_cpu_codelet_update_u21,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u21,
+#endif
+	.nbuffers = 2,
+	.model = &chol_model_21
+};
+
+static void create_task_21(starpu_data_handle dataA, unsigned k, unsigned j, unsigned reclevel, struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = create_task(TAG21_AUX(k, j, reclevel));
+
+	task->cl = &cl21;	
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = starpu_data_get_sub_data(dataA, 2, k, k); 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = starpu_data_get_sub_data(dataA, 2, k, j); 
+	task->buffers[1].mode = STARPU_RW;
+
+	if (j == k+1) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG21_AUX(k, j, reclevel), 2, TAG11_AUX(k, reclevel), TAG22_AUX(k-1, k, j, reclevel));
+	}
+	else {
+		starpu_tag_declare_deps(TAG21_AUX(k, j, reclevel), 1, TAG11_AUX(k, reclevel));
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static starpu_codelet cl22 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = chol_cpu_codelet_update_u22,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u22,
+#endif
+	.nbuffers = 3,
+	.model = &chol_model_22
+};
+
+static void create_task_22(starpu_data_handle dataA, unsigned k, unsigned i, unsigned j, unsigned reclevel, struct starpu_sched_ctx *sched_ctx)
+{
+//	printf("task 22 k,i,j = %d,%d,%d TAG = %llx\n", k,i,j, TAG22_AUX(k,i,j));
+
+	struct starpu_task *task = create_task(TAG22_AUX(k, i, j, reclevel));
+
+	task->cl = &cl22;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = starpu_data_get_sub_data(dataA, 2, k, i); 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = starpu_data_get_sub_data(dataA, 2, k, j); 
+	task->buffers[1].mode = STARPU_R;
+	task->buffers[2].handle = starpu_data_get_sub_data(dataA, 2, i, j); 
+	task->buffers[2].mode = STARPU_RW;
+
+	if ( (i == k + 1) && (j == k +1) ) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG22_AUX(k, i, j, reclevel), 3, TAG22_AUX(k-1, i, j, reclevel), TAG21_AUX(k, i, reclevel), TAG21_AUX(k, j, reclevel));
+	}
+	else {
+		starpu_tag_declare_deps(TAG22_AUX(k, i, j, reclevel), 2, TAG21_AUX(k, i, reclevel), TAG21_AUX(k, j, reclevel));
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+
+
+/*
+ *	code to bootstrap the factorization 
+ *	and construct the DAG
+ */
+
+static void cholesky_grain_rec(float *matA, unsigned size, unsigned ld, unsigned nblocks, unsigned nbigblocks, unsigned reclevel, struct starpu_sched_ctx *sched_ctx)
+{
+	/* create a new codelet */
+	struct starpu_task *entry_task = NULL;
+
+	/* create all the DAG nodes */
+	unsigned i,j,k;
+
+	starpu_data_handle dataA;
+
+	/* monitor and partition the A matrix into blocks :
+	 * one block is now determined by 2 unsigned (i,j) */
+	starpu_matrix_data_register(&dataA, 0, (uintptr_t)matA, ld, size, size, sizeof(float));
+
+	starpu_data_set_sequential_consistency_flag(dataA, 0);
+
+	struct starpu_data_filter f;
+		f.filter_func = starpu_vertical_block_filter_func;
+		f.nchildren = nblocks;
+		f.get_nchildren = NULL;
+		f.get_child_ops = NULL;
+
+	struct starpu_data_filter f2;
+		f2.filter_func = starpu_block_filter_func;
+		f2.nchildren = nblocks;
+		f2.get_nchildren = NULL;
+		f2.get_child_ops = NULL;
+
+	starpu_data_map_filters(dataA, 2, &f, &f2);
+
+	for (k = 0; k < nbigblocks; k++)
+	{
+		struct starpu_task *task = create_task_11(dataA, k, reclevel);
+		/* we defer the launch of the first task */
+		if (k == 0) {
+			entry_task = task;
+		}
+		else {
+		  starpu_task_submit_to_ctx(task, sched_ctx);
+		}
+		
+		for (j = k+1; j<nblocks; j++)
+		{
+		  create_task_21(dataA, k, j, reclevel, sched_ctx);
+
+			for (i = k+1; i<nblocks; i++)
+			{
+				if (i <= j)
+				  create_task_22(dataA, k, i, j, reclevel, sched_ctx);
+			}
+		}
+	}
+
+	/* schedule the codelet */
+	int ret = starpu_task_submit_to_ctx(entry_task, sched_ctx);
+	if (STARPU_UNLIKELY(ret == -ENODEV))
+	{
+		fprintf(stderr, "No worker may execute this task\n");
+		exit(-1);
+	}
+
+	if (nblocks == nbigblocks)
+	{
+		/* stall the application until the end of computations */
+		starpu_tag_wait(TAG11_AUX(nblocks-1, reclevel));
+		starpu_data_unpartition(dataA, 0);
+		return;
+	}
+	else {
+		STARPU_ASSERT(reclevel == 0);
+		unsigned ndeps_tags = (nblocks - nbigblocks)*(nblocks - nbigblocks);
+
+		starpu_tag_t *tag_array = malloc(ndeps_tags*sizeof(starpu_tag_t));
+		STARPU_ASSERT(tag_array);
+
+		unsigned ind = 0;
+		for (i = nbigblocks; i < nblocks; i++)
+		for (j = nbigblocks; j < nblocks; j++)
+		{
+			if (i <= j)
+				tag_array[ind++] = TAG22_AUX(nbigblocks - 1, i, j, reclevel);
+		}
+
+		starpu_tag_wait_array(ind, tag_array);
+
+		free(tag_array);
+
+		starpu_data_unpartition(dataA, 0);
+		starpu_data_unregister(dataA);
+
+		float *newmatA = &matA[nbigblocks*(size/nblocks)*(ld+1)];
+
+		cholesky_grain_rec(newmatA, size/nblocks*(nblocks - nbigblocks), ld, (nblocks - nbigblocks)*2, (nblocks - nbigblocks)*2, reclevel+1, sched_ctx);
+	}
+}
+
+static void initialize_system(float **A, unsigned dim, unsigned pinned)
+{
+  //	starpu_init(NULL);
+
+	starpu_helper_cublas_init();
+
+	if (pinned)
+	{
+		starpu_data_malloc_pinned_if_possible((void **)A, dim*dim*sizeof(float));
+	} 
+	else {
+		*A = malloc(dim*dim*sizeof(float));
+	}
+}
+
+void cholesky_grain(float *matA, unsigned size, unsigned ld, unsigned nblocks, unsigned nbigblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	struct timeval start;
+	struct timeval end;
+
+	gettimeofday(&start, NULL);
+
+	cholesky_grain_rec(matA, size, ld, nblocks, nbigblocks, 0, sched_ctx);
+
+	gettimeofday(&end, NULL);
+
+	double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+	fprintf(stderr, "Computation took (in ms)\n");
+	printf("%2.2f\n", timing/1000);
+
+	double flop = (1.0f*size*size*size)/3.0f;
+	fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
+
+	starpu_helper_cublas_shutdown();
+
+	//	starpu_shutdown();
+}
+
+int run_cholesky_grain_tag(struct starpu_sched_ctx *sched_ctx, int argc, char **argv)
+{
+	/* create a simple definite positive symetric matrix example
+	 *
+	 *	Hilbert matrix : h(i,j) = 1/(i+j+1)
+	 * */
+
+	parse_args(argc, argv);
+
+	float *mat;
+
+	mat = malloc(size*size*sizeof(float));
+	initialize_system(&mat, size, pinned);
+
+	unsigned i,j;
+	for (i = 0; i < size; i++)
+	{
+		for (j = 0; j < size; j++)
+		{
+			mat[j +i*size] = (1.0f/(1.0f+i+j)) + ((i == j)?1.0f*size:0.0f);
+			//mat[j +i*size] = ((i == j)?1.0f*size:0.0f);
+		}
+	}
+
+
+#ifdef CHECK_OUTPUT
+	printf("Input :\n");
+
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+			}
+		}
+		printf("\n");
+	}
+#endif
+
+
+	cholesky_grain(mat, size, size, nblocks, nbigblocks, sched_ctx);
+
+#ifdef CHECK_OUTPUT
+	printf("Results :\n");
+
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+				mat[j+i*size] = 0.0f; // debug
+			}
+		}
+		printf("\n");
+	}
+
+	fprintf(stderr, "compute explicit LLt ...\n");
+	float *test_mat = malloc(size*size*sizeof(float));
+	STARPU_ASSERT(test_mat);
+
+	SSYRK("L", "N", size, size, 1.0f, 
+				mat, size, 0.0f, test_mat, size);
+
+	fprintf(stderr, "comparing results ...\n");
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", test_mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+			}
+		}
+		printf("\n");
+	}
+#endif
+
+	return 0;
+}

examples/cholesky_2ctxs/cholesky/cholesky_implicit.c

@@ -0,0 +1,293 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010, 2011  Université de Bordeaux 1
+ * Copyright (C) 2010  Mehdi Juhoor <mjuhoor@gmail.com>
+ * Copyright (C) 2010, 2011  Centre National de la Recherche Scientifique
+ *
+ * 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.
+ */
+
+#include "cholesky.h"
+
+/*
+ *	Create the codelets
+ */
+
+static starpu_codelet cl11 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.type = STARPU_SEQ,
+	.cpu_func = chol_cpu_codelet_update_u11,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u11,
+#endif
+	.nbuffers = 1,
+	.model = &chol_model_11
+};
+
+static starpu_codelet cl21 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.type = STARPU_SEQ,
+	.cpu_func = chol_cpu_codelet_update_u21,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u21,
+#endif
+	.nbuffers = 2,
+	.model = &chol_model_21
+};
+
+static starpu_codelet cl22 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.type = STARPU_SEQ,
+	.max_parallelism = INT_MAX,
+	.cpu_func = chol_cpu_codelet_update_u22,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u22,
+#endif
+	.nbuffers = 3,
+	.model = &chol_model_22
+};
+
+/*
+ *	code to bootstrap the factorization
+ *	and construct the DAG
+ */
+
+static void callback_turn_spmd_on(void *arg __attribute__ ((unused)))
+{
+	cl22.type = STARPU_SPMD;
+}
+
+static double _cholesky(starpu_data_handle dataA, unsigned nblocks, unsigned sched_ctx, double *timing)
+{
+	struct timeval start;
+	struct timeval end;
+
+	unsigned i,j,k;
+
+	int prio_level = noprio?STARPU_DEFAULT_PRIO:STARPU_MAX_PRIO;
+
+	gettimeofday(&start, NULL);
+
+	/* create all the DAG nodes */
+	for (k = 0; k < nblocks; k++)
+	{
+                starpu_data_handle sdatakk = starpu_data_get_sub_data(dataA, 2, k, k);
+
+		starpu_insert_task(&cl11,
+				   STARPU_PRIORITY, prio_level,
+				   STARPU_RW, sdatakk,
+				   STARPU_CALLBACK, (k == 3*nblocks/4)?callback_turn_spmd_on:NULL,
+				   STARPU_CTX, sched_ctx,
+				   0);
+
+		for (j = k+1; j<nblocks; j++)
+		{
+                        starpu_data_handle sdatakj = starpu_data_get_sub_data(dataA, 2, k, j);
+			starpu_insert_task(&cl21,
+					   STARPU_PRIORITY, (j == k+1)?prio_level:STARPU_DEFAULT_PRIO,
+					   STARPU_R, sdatakk,
+					   STARPU_RW, sdatakj,
+					   STARPU_CTX, sched_ctx,
+					   0);
+
+			for (i = k+1; i<nblocks; i++)
+			{
+				if (i <= j)
+                                {
+					starpu_data_handle sdataki = starpu_data_get_sub_data(dataA, 2, k, i);
+					starpu_data_handle sdataij = starpu_data_get_sub_data(dataA, 2, i, j);
+					
+					starpu_insert_task(&cl22,
+							   STARPU_PRIORITY, ((i == k+1) && (j == k+1))?prio_level:STARPU_DEFAULT_PRIO,
+							   STARPU_R, sdataki,
+							   STARPU_R, sdatakj,
+							   STARPU_RW, sdataij,
+							   STARPU_CTX, sched_ctx,
+							   0);
+                                }
+			}
+		}
+	}
+
+	if(sched_ctx != 0)
+		starpu_wait_for_all_tasks_of_sched_ctx(sched_ctx);
+	else
+		starpu_task_wait_for_all();
+
+	starpu_data_unpartition(dataA, 0);
+
+	gettimeofday(&end, NULL);
+
+	(*timing) = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+
+	unsigned long n = starpu_matrix_get_nx(dataA);
+
+	double flop = (1.0f*n*n*n)/3.0f;
+
+	double gflops = (flop/(*timing)/1000.0f);
+	(*timing) /= 1000000.0f; //sec
+	//	(*timing) /= 60.0f; //min
+	return gflops;
+}
+
+static double cholesky(float *matA, unsigned size, unsigned ld, unsigned nblocks, unsigned sched_ctx, double *timing)
+{
+	starpu_data_handle dataA;
+
+	/* monitor and partition the A matrix into blocks :
+	 * one block is now determined by 2 unsigned (i,j) */
+	starpu_matrix_data_register(&dataA, 0, (uintptr_t)matA, ld, size, size, sizeof(float));
+
+	struct starpu_data_filter f = {
+		.filter_func = starpu_vertical_block_filter_func,
+		.nchildren = nblocks
+	};
+
+	struct starpu_data_filter f2 = {
+		.filter_func = starpu_block_filter_func,
+		.nchildren = nblocks
+	};
+
+	starpu_data_map_filters(dataA, 2, &f, &f2);
+	double gflops = _cholesky(dataA, nblocks, sched_ctx, timing);
+	starpu_data_unregister(dataA);
+	return gflops;
+}
+
+double run_cholesky_implicit(unsigned sched_ctx, int start, int argc, char **argv, double *timing, pthread_barrier_t *barrier)
+{
+	/* create a simple definite positive symetric matrix example
+	 *
+	 *	Hilbert matrix : h(i,j) = 1/(i+j+1)
+	 * */
+
+	unsigned size = 4 * 1024;
+	unsigned nblocks = 16;
+	parse_args_ctx(start, argc, argv, &size, &nblocks);
+
+	//	starpu_init(NULL);
+
+	//	starpu_helper_cublas_init();
+
+	float *mat;
+
+	starpu_malloc((void **)&mat, (size_t)size*size*sizeof(float));
+
+	unsigned i,j;
+	for (i = 0; i < size; i++)
+	{
+		for (j = 0; j < size; j++)
+		{
+			mat[j +i*size] = (1.0f/(1.0f+i+j)) + ((i == j)?1.0f*size:0.0f);
+			//mat[j +i*size] = ((i == j)?1.0f*size:0.0f);
+		}
+	}
+
+//#define PRINT_OUTPUT
+#ifdef PRINT_OUTPUT
+	printf("Input :\n");
+
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+			}
+		}
+		printf("\n");
+	}
+#endif
+	//	if(barrier != NULL)
+	//	  pthread_barrier_wait(barrier);
+	double gflops = cholesky(mat, size, size, nblocks, sched_ctx, timing);
+
+#ifdef PRINT_OUTPUT
+	printf("Results :\n");
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+				mat[j+i*size] = 0.0f; // debug
+			}
+		}
+		printf("\n");
+	}
+#endif
+
+	if (check)
+	{
+		fprintf(stderr, "compute explicit LLt ...\n");
+		for (j = 0; j < size; j++)
+		{
+			for (i = 0; i < size; i++)
+			{
+				if (i > j) {
+					mat[j+i*size] = 0.0f; // debug
+				}
+			}
+		}
+		float *test_mat = malloc(size*size*sizeof(float));
+		STARPU_ASSERT(test_mat);
+	
+		SSYRK("L", "N", size, size, 1.0f,
+					mat, size, 0.0f, test_mat, size);
+	
+		fprintf(stderr, "comparing results ...\n");
+#ifdef PRINT_OUTPUT
+		for (j = 0; j < size; j++)
+		{
+			for (i = 0; i < size; i++)
+			{
+				if (i <= j) {
+					printf("%2.2f\t", test_mat[j +i*size]);
+				}
+				else {
+					printf(".\t");
+				}
+			}
+			printf("\n");
+		}
+#endif
+	
+		for (j = 0; j < size; j++)
+		{
+			for (i = 0; i < size; i++)
+			{
+				if (i <= j) {
+	                                float orig = (1.0f/(1.0f+i+j)) + ((i == j)?1.0f*size:0.0f);
+	                                float err = abs(test_mat[j +i*size] - orig);
+	                                if (err > 0.00001) {
+	                                        fprintf(stderr, "Error[%d, %d] --> %2.2f != %2.2f (err %2.2f)\n", i, j, test_mat[j +i*size], orig, err);
+	                                        assert(0);
+	                                }
+	                        }
+			}
+	        }
+	}
+	starpu_free((void *)mat);
+	//	starpu_helper_cublas_shutdown();
+	//	starpu_shutdown();
+
+	return gflops;
+}

examples/cholesky_2ctxs/cholesky/cholesky_implicit_all_machine.c

@@ -0,0 +1,280 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010, 2011  Université de Bordeaux 1
+ * Copyright (C) 2010  Mehdi Juhoor <mjuhoor@gmail.com>
+ * Copyright (C) 2010, 2011  Centre National de la Recherche Scientifique
+ *
+ * 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.
+ */
+
+#include "cholesky.h"
+
+/*
+ *	Create the codelets
+ */
+
+static starpu_codelet cl11 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.type = STARPU_SEQ,
+	.cpu_func = chol_cpu_codelet_update_u11,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u11,
+#endif
+	.nbuffers = 1,
+	.model = &chol_model_11
+};
+
+static starpu_codelet cl21 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.type = STARPU_SEQ,
+	.cpu_func = chol_cpu_codelet_update_u21,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u21,
+#endif
+	.nbuffers = 2,
+	.model = &chol_model_21
+};
+
+static starpu_codelet cl22 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.type = STARPU_SEQ,
+	.max_parallelism = INT_MAX,
+	.cpu_func = chol_cpu_codelet_update_u22,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u22,
+#endif
+	.nbuffers = 3,
+	.model = &chol_model_22
+};
+
+/*
+ *	code to bootstrap the factorization
+ *	and construct the DAG
+ */
+
+static void callback_turn_spmd_on(void *arg __attribute__ ((unused)))
+{
+	cl22.type = STARPU_SPMD;
+}
+
+static double _cholesky(starpu_data_handle dataA, unsigned nblocks)
+{
+	struct timeval start;
+	struct timeval end;
+
+	unsigned i,j,k;
+
+	int prio_level = noprio?STARPU_DEFAULT_PRIO:STARPU_MAX_PRIO;
+
+	gettimeofday(&start, NULL);
+
+	/* create all the DAG nodes */
+	for (k = 0; k < nblocks; k++)
+	{
+                starpu_data_handle sdatakk = starpu_data_get_sub_data(dataA, 2, k, k);
+
+                starpu_insert_task(&cl11,
+                                   STARPU_PRIORITY, prio_level,
+                                   STARPU_RW, sdatakk,
+				   STARPU_CALLBACK, (k == 3*nblocks/4)?callback_turn_spmd_on:NULL,
+                                   0);
+
+		for (j = k+1; j<nblocks; j++)
+		{
+                        starpu_data_handle sdatakj = starpu_data_get_sub_data(dataA, 2, k, j);
+
+                        starpu_insert_task(&cl21,
+                                           STARPU_PRIORITY, (j == k+1)?prio_level:STARPU_DEFAULT_PRIO,
+                                           STARPU_R, sdatakk,
+                                           STARPU_RW, sdatakj,
+                                           0);
+
+			for (i = k+1; i<nblocks; i++)
+			{
+				if (i <= j)
+                                {
+					starpu_data_handle sdataki = starpu_data_get_sub_data(dataA, 2, k, i);
+					starpu_data_handle sdataij = starpu_data_get_sub_data(dataA, 2, i, j);
+					
+					starpu_insert_task(&cl22,
+                                                           STARPU_PRIORITY, ((i == k+1) && (j == k+1))?prio_level:STARPU_DEFAULT_PRIO,
+                                                           STARPU_R, sdataki,
+                                                           STARPU_R, sdatakj,
+                                                           STARPU_RW, sdataij,
+                                                           0);
+                                }
+			}
+		}
+	}
+
+	starpu_task_wait_for_all();
+
+	starpu_data_unpartition(dataA, 0);
+
+	gettimeofday(&end, NULL);
+
+	double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+	unsigned long n = starpu_matrix_get_nx(dataA);
+
+	double flop = (1.0f*n*n*n)/3.0f;
+	return (flop/timing/1000.0f);
+}
+
+static double cholesky(float *matA, unsigned size, unsigned ld, unsigned nblocks)
+{
+	starpu_data_handle dataA;
+
+	/* monitor and partition the A matrix into blocks :
+	 * one block is now determined by 2 unsigned (i,j) */
+	starpu_matrix_data_register(&dataA, 0, (uintptr_t)matA, ld, size, size, sizeof(float));
+
+	struct starpu_data_filter f;
+		f.filter_func = starpu_vertical_block_filter_func;
+		f.nchildren = nblocks;
+		f.get_nchildren = NULL;
+		f.get_child_ops = NULL;
+
+	struct starpu_data_filter f2;
+		f2.filter_func = starpu_block_filter_func;
+		f2.nchildren = nblocks;
+		f2.get_nchildren = NULL;
+		f2.get_child_ops = NULL;
+
+	starpu_data_map_filters(dataA, 2, &f, &f2);
+
+	return _cholesky(dataA, nblocks);
+}
+
+double run_cholesky_implicit_all_machine(int argc, char **argv)
+{
+	/* create a simple definite positive symetric matrix example
+	 *
+	 *	Hilbert matrix : h(i,j) = 1/(i+j+1)
+	 * */
+
+	parse_args(argc, argv);
+
+	//	starpu_init(NULL);
+
+	//	starpu_helper_cublas_init();
+
+	float *mat;
+	starpu_data_malloc_pinned_if_possible((void **)&mat, (size_t)size*size*sizeof(float));
+
+	unsigned i,j;
+	for (i = 0; i < size; i++)
+	{
+		for (j = 0; j < size; j++)
+		{
+			mat[j +i*size] = (1.0f/(1.0f+i+j)) + ((i == j)?1.0f*size:0.0f);
+			//mat[j +i*size] = ((i == j)?1.0f*size:0.0f);
+		}
+	}
+
+//#define PRINT_OUTPUT
+#ifdef PRINT_OUTPUT
+	printf("Input :\n");
+
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+			}
+		}
+		printf("\n");
+	}
+#endif
+
+	double gflops = cholesky(mat, size, size, nblocks);
+
+#ifdef PRINT_OUTPUT
+	printf("Results :\n");
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+				mat[j+i*size] = 0.0f; // debug
+			}
+		}
+		printf("\n");
+	}
+#endif
+
+	if (check)
+	{
+		fprintf(stderr, "compute explicit LLt ...\n");
+		for (j = 0; j < size; j++)
+		{
+			for (i = 0; i < size; i++)
+			{
+				if (i > j) {
+					mat[j+i*size] = 0.0f; // debug
+				}
+			}
+		}
+		float *test_mat = malloc(size*size*sizeof(float));
+		STARPU_ASSERT(test_mat);
+	
+		SSYRK("L", "N", size, size, 1.0f,
+					mat, size, 0.0f, test_mat, size);
+	
+		fprintf(stderr, "comparing results ...\n");
+#ifdef PRINT_OUTPUT
+		for (j = 0; j < size; j++)
+		{
+			for (i = 0; i < size; i++)
+			{
+				if (i <= j) {
+					printf("%2.2f\t", test_mat[j +i*size]);
+				}
+				else {
+					printf(".\t");
+				}
+			}
+			printf("\n");
+		}
+#endif
+	
+		for (j = 0; j < size; j++)
+		{
+			for (i = 0; i < size; i++)
+			{
+				if (i <= j) {
+	                                float orig = (1.0f/(1.0f+i+j)) + ((i == j)?1.0f*size:0.0f);
+	                                float err = abs(test_mat[j +i*size] - orig);
+	                                if (err > 0.00001) {
+	                                        fprintf(stderr, "Error[%d, %d] --> %2.2f != %2.2f (err %2.2f)\n", i, j, test_mat[j +i*size], orig, err);
+	                                        assert(0);
+	                                }
+	                        }
+			}
+	        }
+	}
+
+	//	starpu_helper_cublas_shutdown();
+	//	starpu_shutdown();
+
+	return gflops;
+}

examples/cholesky_2ctxs/cholesky/cholesky_kernels.c

@@ -0,0 +1,230 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010, 2011  Université de Bordeaux 1
+ * Copyright (C) 2010  Centre National de la Recherche Scientifique
+ *
+ * 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.
+ */
+
+#include <starpu_config.h>
+#include "cholesky.h"
+#include "../../common/blas.h"
+#ifdef STARPU_USE_CUDA
+#include <starpu_cuda.h>
+#endif
+
+/*
+ *   U22 
+ */
+
+static inline void chol_common_cpu_codelet_update_u22(void *descr[], int s, __attribute__((unused)) void *_args)
+{
+	//printf("22\n");
+	float *left 	= (float *)STARPU_MATRIX_GET_PTR(descr[0]);
+	float *right 	= (float *)STARPU_MATRIX_GET_PTR(descr[1]);
+	float *center 	= (float *)STARPU_MATRIX_GET_PTR(descr[2]);
+
+	unsigned dx = STARPU_MATRIX_GET_NY(descr[2]);
+	unsigned dy = STARPU_MATRIX_GET_NX(descr[2]);
+	unsigned dz = STARPU_MATRIX_GET_NY(descr[0]);
+
+	unsigned ld21 = STARPU_MATRIX_GET_LD(descr[0]);
+	unsigned ld12 = STARPU_MATRIX_GET_LD(descr[1]);
+	unsigned ld22 = STARPU_MATRIX_GET_LD(descr[2]);
+
+	if (s == 0)
+	{
+		int worker_size = starpu_combined_worker_get_size();
+
+		if (worker_size == 1)
+		{
+			/* Sequential CPU kernel */
+			SGEMM("N", "T", dy, dx, dz, -1.0f, left, ld21, 
+				right, ld12, 1.0f, center, ld22);
+		}
+		else {
+			/* Parallel CPU kernel */
+			int rank = starpu_combined_worker_get_rank();
+
+			int block_size = (dx + worker_size - 1)/worker_size;
+			int new_dx = STARPU_MIN(dx, block_size*(rank+1)) - block_size*rank;
+			
+			float *new_left = &left[block_size*rank];
+			float *new_center = &center[block_size*rank];
+
+			SGEMM("N", "T", dy, new_dx, dz, -1.0f, new_left, ld21, 
+				right, ld12, 1.0f, new_center, ld22);
+		}
+	}
+	else
+	{
+		/* CUDA kernel */
+#ifdef STARPU_USE_CUDA
+		cublasSgemm('n', 't', dy, dx, dz, 
+				-1.0f, left, ld21, right, ld12, 
+				 1.0f, center, ld22);
+		cudaStreamSynchronize(starpu_cuda_get_local_stream());
+#endif
+
+	}
+}
+
+void chol_cpu_codelet_update_u22(void *descr[], void *_args)
+{
+	chol_common_cpu_codelet_update_u22(descr, 0, _args);
+}
+
+#ifdef STARPU_USE_CUDA
+void chol_cublas_codelet_update_u22(void *descr[], void *_args)
+{
+	chol_common_cpu_codelet_update_u22(descr, 1, _args);
+}
+#endif// STARPU_USE_CUDA
+
+/* 
+ * U21
+ */
+
+static inline void chol_common_codelet_update_u21(void *descr[], int s, __attribute__((unused)) void *_args)
+{
+//	printf("21\n");
+	float *sub11;
+	float *sub21;
+
+	sub11 = (float *)STARPU_MATRIX_GET_PTR(descr[0]);
+	sub21 = (float *)STARPU_MATRIX_GET_PTR(descr[1]);
+
+	unsigned ld11 = STARPU_MATRIX_GET_LD(descr[0]);
+	unsigned ld21 = STARPU_MATRIX_GET_LD(descr[1]);
+
+	unsigned nx21 = STARPU_MATRIX_GET_NY(descr[1]);
+	unsigned ny21 = STARPU_MATRIX_GET_NX(descr[1]);
+
+	switch (s) {
+		case 0:
+			STRSM("R", "L", "T", "N", nx21, ny21, 1.0f, sub11, ld11, sub21, ld21);
+			break;
+#ifdef STARPU_USE_CUDA
+		case 1:
+			cublasStrsm('R', 'L', 'T', 'N', nx21, ny21, 1.0f, sub11, ld11, sub21, ld21);
+			cudaStreamSynchronize(starpu_cuda_get_local_stream());
+			break;
+#endif
+		default:
+			STARPU_ABORT();
+			break;
+	}
+}
+
+void chol_cpu_codelet_update_u21(void *descr[], void *_args)
+{
+	 chol_common_codelet_update_u21(descr, 0, _args);
+}
+
+#ifdef STARPU_USE_CUDA
+void chol_cublas_codelet_update_u21(void *descr[], void *_args)
+{
+	chol_common_codelet_update_u21(descr, 1, _args);
+}
+#endif 
+
+/*
+ *	U11
+ */
+
+static inline void chol_common_codelet_update_u11(void *descr[], int s, __attribute__((unused)) void *_args) 
+{
+//	printf("11\n");
+	float *sub11;
+
+	sub11 = (float *)STARPU_MATRIX_GET_PTR(descr[0]); 
+
+	unsigned nx = STARPU_MATRIX_GET_NY(descr[0]);
+	unsigned ld = STARPU_MATRIX_GET_LD(descr[0]);
+
+	unsigned z;
+
+	switch (s) {
+		case 0:
+
+			/*
+			 *	- alpha 11 <- lambda 11 = sqrt(alpha11)
+			 *	- alpha 21 <- l 21	= alpha 21 / lambda 11
+			 *	- A22 <- A22 - l21 trans(l21)
+			 */
+
+			for (z = 0; z < nx; z++)
+			{
+				float lambda11;
+				lambda11 = sqrt(sub11[z+z*ld]);
+				sub11[z+z*ld] = lambda11;
+
+				STARPU_ASSERT(lambda11 != 0.0f);
+		
+				SSCAL(nx - z - 1, 1.0f/lambda11, &sub11[(z+1)+z*ld], 1);
+		
+				SSYR("L", nx - z - 1, -1.0f, 
+							&sub11[(z+1)+z*ld], 1,
+							&sub11[(z+1)+(z+1)*ld], ld);
+			}
+			break;
+#ifdef STARPU_USE_CUDA
+		case 1:
+			{
+			float *lambda11;
+			cudaHostAlloc((void **)&lambda11, sizeof(float), 0);
+
+			for (z = 0; z < nx; z++)
+			{
+
+				cudaMemcpyAsync(lambda11, &sub11[z+z*ld], sizeof(float), cudaMemcpyDeviceToHost, starpu_cuda_get_local_stream());
+				cudaStreamSynchronize(starpu_cuda_get_local_stream());
+
+				STARPU_ASSERT(*lambda11 != 0.0f);
+				
+				*lambda11 = sqrt(*lambda11);
+
+//				cublasSetVector(1, sizeof(float), lambda11, sizeof(float), &sub11[z+z*ld], sizeof(float));
+				cudaMemcpyAsync(&sub11[z+z*ld], lambda11, sizeof(float), cudaMemcpyHostToDevice, starpu_cuda_get_local_stream());
+
+				cublasSscal(nx - z - 1, 1.0f/(*lambda11), &sub11[(z+1)+z*ld], 1);
+
+				cublasSsyr('U', nx - z - 1, -1.0f,
+							&sub11[(z+1)+z*ld], 1,
+							&sub11[(z+1)+(z+1)*ld], ld);
+			}
+
+			cudaStreamSynchronize(starpu_cuda_get_local_stream());
+			cudaFreeHost(lambda11);
+			}
+		
+
+			break;
+#endif
+		default:
+			STARPU_ABORT();
+			break;
+	}
+}
+
+
+void chol_cpu_codelet_update_u11(void *descr[], void *_args)
+{
+	chol_common_codelet_update_u11(descr, 0, _args);
+}
+
+#ifdef STARPU_USE_CUDA
+void chol_cublas_codelet_update_u11(void *descr[], void *_args)
+{
+	chol_common_codelet_update_u11(descr, 1, _args);
+}
+#endif// STARPU_USE_CUDA

examples/cholesky_2ctxs/cholesky/cholesky_models.c

@@ -0,0 +1,152 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010  Université de Bordeaux 1
+ * Copyright (C) 2010  Centre National de la Recherche Scientifique
+ *
+ * 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.
+ */
+
+/*
+ * As a convention, in that file, descr[0] is represented by A,
+ * 				  descr[1] is B ...
+ */
+
+/*
+ *	Number of flops of Gemm 
+ */
+
+#include <starpu.h>
+
+//#define USE_PERTURBATION	1
+
+#ifdef USE_PERTURBATION
+#define PERTURBATE(a)	((starpu_drand48()*2.0f*(AMPL) + 1.0f - (AMPL))*(a))
+#else
+#define PERTURBATE(a)	(a)
+#endif
+
+static double cpu_chol_task_11_cost(starpu_buffer_descr *descr)
+{
+	uint32_t n;
+
+	n = starpu_matrix_get_nx(descr[0].handle);
+
+	double cost = (((double)(n)*n*n)/1000.0f*0.894/0.79176);
+
+#ifdef STARPU_MODEL_DEBUG
+	printf("cpu_chol_task_11_cost n %d cost %e\n", n, cost);
+#endif
+
+	return PERTURBATE(cost);
+}
+
+static double cuda_chol_task_11_cost(starpu_buffer_descr *descr)
+{
+	uint32_t n;
+
+	n = starpu_matrix_get_nx(descr[0].handle);
+
+	double cost = (((double)(n)*n*n)/50.0f/10.75/5.088633/0.9883);
+
+#ifdef STARPU_MODEL_DEBUG
+	printf("cuda_chol_task_11_cost n %d cost %e\n", n, cost);
+#endif
+
+	return PERTURBATE(cost);
+}
+
+static double cpu_chol_task_21_cost(starpu_buffer_descr *descr)
+{
+	uint32_t n;
+
+	n = starpu_matrix_get_nx(descr[0].handle);
+
+	double cost = (((double)(n)*n*n)/7706.674/0.95/0.9965);
+
+#ifdef STARPU_MODEL_DEBUG
+	printf("cpu_chol_task_21_cost n %d cost %e\n", n, cost);
+#endif
+
+	return PERTURBATE(cost);
+}
+
+static double cuda_chol_task_21_cost(starpu_buffer_descr *descr)
+{
+	uint32_t n;
+
+	n = starpu_matrix_get_nx(descr[0].handle);
+
+	double cost = (((double)(n)*n*n)/50.0f/10.75/87.29520);
+
+#ifdef STARPU_MODEL_DEBUG
+	printf("cuda_chol_task_21_cost n %d cost %e\n", n, cost);
+#endif
+
+	return PERTURBATE(cost);
+}
+
+static double cpu_chol_task_22_cost(starpu_buffer_descr *descr)
+{
+	uint32_t n;
+
+	n = starpu_matrix_get_nx(descr[0].handle);
+
+	double cost = (((double)(n)*n*n)/50.0f/10.75/8.0760);
+
+#ifdef STARPU_MODEL_DEBUG
+	printf("cpu_chol_task_22_cost n %d cost %e\n", n, cost);
+#endif
+
+	return PERTURBATE(cost);
+}
+
+static double cuda_chol_task_22_cost(starpu_buffer_descr *descr)
+{
+	uint32_t n;
+
+	n = starpu_matrix_get_nx(descr[0].handle);
+
+	double cost = (((double)(n)*n*n)/50.0f/10.75/76.30666);
+
+#ifdef STARPU_MODEL_DEBUG
+	printf("cuda_chol_task_22_cost n %d cost %e\n", n, cost);
+#endif
+
+	return PERTURBATE(cost);
+}
+
+struct starpu_perfmodel_t chol_model_11 = {
+	.per_arch = { 
+		[STARPU_CPU_DEFAULT] = { .cost_model = cpu_chol_task_11_cost },
+		[STARPU_CUDA_DEFAULT] = { .cost_model = cuda_chol_task_11_cost }
+	},
+	.type = STARPU_HISTORY_BASED,
+	.symbol = "chol_model_11"
+};
+
+struct starpu_perfmodel_t chol_model_21 = {
+	.per_arch = { 
+		[STARPU_CPU_DEFAULT] = { .cost_model = cpu_chol_task_21_cost },
+		[STARPU_CUDA_DEFAULT] = { .cost_model = cuda_chol_task_21_cost }
+	},
+	.type = STARPU_HISTORY_BASED,
+	.symbol = "chol_model_21"
+};
+
+struct starpu_perfmodel_t chol_model_22 = {
+	.per_arch = { 
+		[STARPU_CPU_DEFAULT] = { .cost_model = cpu_chol_task_22_cost },
+		[STARPU_CUDA_DEFAULT] = { .cost_model = cuda_chol_task_22_cost }
+	},
+	.type = STARPU_HISTORY_BASED,
+	.symbol = "chol_model_22"
+};

examples/cholesky_2ctxs/cholesky/cholesky_tag.c

@@ -0,0 +1,370 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010, 2011  Université de Bordeaux 1
+ * Copyright (C) 2010  Mehdi Juhoor <mjuhoor@gmail.com>
+ * Copyright (C) 2010  Centre National de la Recherche Scientifique
+ *
+ * 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.
+ */
+
+#include "cholesky.h"
+
+/*
+ *	Some useful functions
+ */
+
+static struct starpu_task *create_task(starpu_tag_t id)
+{
+	struct starpu_task *task = starpu_task_create();
+		task->cl_arg = NULL;
+		task->use_tag = 1;
+		task->tag_id = id;
+
+	return task;
+}
+
+/*
+ *	Create the codelets
+ */
+
+static starpu_codelet cl11 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = chol_cpu_codelet_update_u11,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u11,
+#endif
+	.nbuffers = 1,
+	.model = &chol_model_11
+};
+
+static struct starpu_task * create_task_11(starpu_data_handle dataA, unsigned k)
+{
+//	printf("task 11 k = %d TAG = %llx\n", k, (TAG11(k)));
+
+	struct starpu_task *task = create_task(TAG11(k));
+	
+	task->cl = &cl11;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = starpu_data_get_sub_data(dataA, 2, k, k);
+	task->buffers[0].mode = STARPU_RW;
+
+	/* this is an important task */
+	if (!noprio)
+		task->priority = STARPU_MAX_PRIO;
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG11(k), 1, TAG22(k-1, k, k));
+	}
+
+	return task;
+}
+
+static starpu_codelet cl21 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = chol_cpu_codelet_update_u21,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u21,
+#endif
+	.nbuffers = 2,
+	.model = &chol_model_21
+};
+
+static void create_task_21(starpu_data_handle dataA, unsigned k, unsigned j, struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = create_task(TAG21(k, j));
+
+	task->cl = &cl21;	
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = starpu_data_get_sub_data(dataA, 2, k, k); 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = starpu_data_get_sub_data(dataA, 2, k, j); 
+	task->buffers[1].mode = STARPU_RW;
+
+	if (!noprio && (j == k+1)) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG21(k, j), 2, TAG11(k), TAG22(k-1, k, j));
+	}
+	else {
+		starpu_tag_declare_deps(TAG21(k, j), 1, TAG11(k));
+	}
+
+	int ret = starpu_task_submit_to_ctx(task, sched_ctx);
+        if (STARPU_UNLIKELY(ret == -ENODEV)) {
+                fprintf(stderr, "No worker may execute this task\n");
+                exit(0);
+        }
+
+}
+
+static starpu_codelet cl22 =
+{
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = chol_cpu_codelet_update_u22,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u22,
+#endif
+	.nbuffers = 3,
+	.model = &chol_model_22
+};
+
+static void create_task_22(starpu_data_handle dataA, unsigned k, unsigned i, unsigned j, struct starpu_sched_ctx *sched_ctx)
+{
+//	printf("task 22 k,i,j = %d,%d,%d TAG = %llx\n", k,i,j, TAG22(k,i,j));
+
+	struct starpu_task *task = create_task(TAG22(k, i, j));
+
+	task->cl = &cl22;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = starpu_data_get_sub_data(dataA, 2, k, i); 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = starpu_data_get_sub_data(dataA, 2, k, j); 
+	task->buffers[1].mode = STARPU_R;
+	task->buffers[2].handle = starpu_data_get_sub_data(dataA, 2, i, j); 
+	task->buffers[2].mode = STARPU_RW;
+
+	if (!noprio && (i == k + 1) && (j == k +1) ) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG22(k, i, j), 3, TAG22(k-1, i, j), TAG21(k, i), TAG21(k, j));
+	}
+	else {
+		starpu_tag_declare_deps(TAG22(k, i, j), 2, TAG21(k, i), TAG21(k, j));
+	}
+
+	int ret = starpu_task_submit_to_ctx(task, sched_ctx);
+        if (STARPU_UNLIKELY(ret == -ENODEV)) {
+                fprintf(stderr, "No worker may execute this task\n");
+                exit(0);
+        }
+}
+
+
+
+/*
+ *	code to bootstrap the factorization 
+ *	and construct the DAG
+ */
+
+static void _cholesky(starpu_data_handle dataA, unsigned nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	struct timeval start;
+	struct timeval end;
+
+	struct starpu_task *entry_task = NULL;
+
+	/* create all the DAG nodes */
+	unsigned i,j,k;
+
+	gettimeofday(&start, NULL);
+
+	for (k = 0; k < nblocks; k++)
+	{
+		struct starpu_task *task = create_task_11(dataA, k);
+		/* we defer the launch of the first task */
+		if (k == 0) {
+			entry_task = task;
+		}
+		else {
+		  int ret = starpu_task_submit_to_ctx(task, sched_ctx);
+                        if (STARPU_UNLIKELY(ret == -ENODEV)) {
+                                fprintf(stderr, "No worker may execute this task\n");
+                                exit(0);
+                        }
+
+		}
+		
+		for (j = k+1; j<nblocks; j++)
+		{
+		  create_task_21(dataA, k, j, sched_ctx);
+
+			for (i = k+1; i<nblocks; i++)
+			{
+				if (i <= j)
+				  create_task_22(dataA, k, i, j, sched_ctx);
+			}
+		}
+	}
+
+	/* schedule the codelet */
+	int ret = starpu_task_submit_to_ctx(entry_task, sched_ctx);
+        if (STARPU_UNLIKELY(ret == -ENODEV)) {
+                fprintf(stderr, "No worker may execute this task\n");
+                exit(0);
+        }
+
+
+	/* stall the application until the end of computations */
+	starpu_tag_wait(TAG11(nblocks-1));
+
+	starpu_data_unpartition(dataA, 0);
+
+	gettimeofday(&end, NULL);
+
+
+	double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+	fprintf(stderr, "Computation took (in ms)\n");
+	printf("%2.2f\n", timing/1000);
+
+	unsigned n = starpu_matrix_get_nx(dataA);
+
+	double flop = (1.0f*n*n*n)/3.0f;
+	fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
+}
+
+static void initialize_system(float **A, unsigned dim, unsigned pinned)
+{
+  //	starpu_init(NULL);
+	
+	starpu_helper_cublas_init();
+
+	if (pinned)
+	{
+		starpu_data_malloc_pinned_if_possible((void **)A, (size_t)dim*dim*sizeof(float));
+	} 
+	else {
+		*A = malloc(dim*dim*sizeof(float));
+	}
+}
+
+static void cholesky(float *matA, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	starpu_data_handle dataA;
+
+	/* monitor and partition the A matrix into blocks :
+	 * one block is now determined by 2 unsigned (i,j) */
+	starpu_matrix_data_register(&dataA, 0, (uintptr_t)matA, ld, size, size, sizeof(float));
+
+	starpu_data_set_sequential_consistency_flag(dataA, 0);
+
+	struct starpu_data_filter f;
+		f.filter_func = starpu_vertical_block_filter_func;
+		f.nchildren = nblocks;
+		f.get_nchildren = NULL;
+		f.get_child_ops = NULL;
+
+	struct starpu_data_filter f2;
+		f2.filter_func = starpu_block_filter_func;
+		f2.nchildren = nblocks;
+		f2.get_nchildren = NULL;
+		f2.get_child_ops = NULL;
+
+	starpu_data_map_filters(dataA, 2, &f, &f2);
+
+	_cholesky(dataA, nblocks, sched_ctx);
+
+	starpu_helper_cublas_shutdown();
+
+	//	starpu_shutdown();
+}
+
+int run_cholesky_tag(struct starpu_sched_ctx *sched_ctx, int argc, char **argv)
+{
+	/* create a simple definite positive symetric matrix example
+	 *
+	 *	Hilbert matrix : h(i,j) = 1/(i+j+1)
+	 * */
+
+	parse_args(argc, argv);
+
+	float *mat;
+
+	mat = malloc(size*size*sizeof(float));
+	initialize_system(&mat, size, pinned);
+
+	unsigned i,j;
+	for (i = 0; i < size; i++)
+	{
+		for (j = 0; j < size; j++)
+		{
+			mat[j +i*size] = (1.0f/(1.0f+i+j)) + ((i == j)?1.0f*size:0.0f);
+			//mat[j +i*size] = ((i == j)?1.0f*size:0.0f);
+		}
+	}
+
+
+#ifdef CHECK_OUTPUT
+	printf("Input :\n");
+
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+			}
+		}
+		printf("\n");
+	}
+#endif
+
+
+	cholesky(mat, size, size, nblocks, sched_ctx);
+
+#ifdef CHECK_OUTPUT
+	printf("Results :\n");
+
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+				mat[j+i*size] = 0.0f; // debug
+			}
+		}
+		printf("\n");
+	}
+
+	fprintf(stderr, "compute explicit LLt ...\n");
+	float *test_mat = malloc(size*size*sizeof(float));
+	STARPU_ASSERT(test_mat);
+
+	SSYRK("L", "N", size, size, 1.0f, 
+				mat, size, 0.0f, test_mat, size);
+
+	fprintf(stderr, "comparing results ...\n");
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < size; i++)
+		{
+			if (i <= j) {
+				printf("%2.2f\t", test_mat[j +i*size]);
+			}
+			else {
+				printf(".\t");
+			}
+		}
+		printf("\n");
+	}
+#endif
+
+	return 0;
+}

examples/cholesky_2ctxs/cholesky/cholesky_tile_tag.c

@@ -0,0 +1,307 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010, 2011  Université de Bordeaux 1
+ * Copyright (C) 2010  Centre National de la Recherche Scientifique
+ *
+ * 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.
+ */
+
+#include "cholesky.h"
+
+/*
+ *	Some useful functions
+ */
+
+/* A [ y ] [ x ] */
+float *ch_A[NMAXBLOCKS][NMAXBLOCKS];
+starpu_data_handle ch_A_state[NMAXBLOCKS][NMAXBLOCKS];
+
+static struct starpu_task *create_task(starpu_tag_t id)
+{
+	struct starpu_task *task = starpu_task_create();
+		task->cl_arg = NULL;
+		task->use_tag = 1;
+		task->tag_id = id;
+
+	return task;
+}
+
+/*
+ *	Create the codelets
+ */
+
+static starpu_codelet cl11 =
+{
+	.where = STARPU_CPU|STARPU_CUDA|STARPU_GORDON,
+	.cpu_func = chol_cpu_codelet_update_u11,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u11,
+#endif
+#ifdef STARPU_USE_GORDON
+#ifdef SPU_FUNC_POTRF
+	.gordon_func = SPU_FUNC_POTRF,
+#else
+#warning SPU_FUNC_POTRF is not available
+#endif
+#endif
+	.nbuffers = 1,
+	.model = &chol_model_11
+};
+
+static struct starpu_task * create_task_11(unsigned k, unsigned nblocks)
+{
+//	printf("task 11 k = %d TAG = %llx\n", k, (TAG11(k)));
+
+	struct starpu_task *task = create_task(TAG11(k));
+	
+	task->cl = &cl11;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = ch_A_state[k][k];
+	task->buffers[0].mode = STARPU_RW;
+
+	/* this is an important task */
+	task->priority = STARPU_MAX_PRIO;
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG11(k), 1, TAG22(k-1, k, k));
+	}
+
+	return task;
+}
+
+static starpu_codelet cl21 =
+{
+	.where = STARPU_CPU|STARPU_CUDA|STARPU_GORDON,
+	.cpu_func = chol_cpu_codelet_update_u21,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u21,
+#endif
+#ifdef STARPU_USE_GORDON
+#ifdef SPU_FUNC_STRSM
+	.gordon_func = SPU_FUNC_STRSM,
+#else
+#warning SPU_FUNC_STRSM is not available
+#endif
+#endif
+	.nbuffers = 2,
+	.model = &chol_model_21
+};
+
+static void create_task_21(unsigned k, unsigned j, struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = create_task(TAG21(k, j));
+
+	task->cl = &cl21;	
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = ch_A_state[k][k]; 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = ch_A_state[j][k]; 
+	task->buffers[1].mode = STARPU_RW;
+
+	if (j == k+1) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG21(k, j), 2, TAG11(k), TAG22(k-1, k, j));
+	}
+	else {
+		starpu_tag_declare_deps(TAG21(k, j), 1, TAG11(k));
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static starpu_codelet cl22 =
+{
+	.where = STARPU_CPU|STARPU_CUDA|STARPU_GORDON,
+	.cpu_func = chol_cpu_codelet_update_u22,
+#ifdef STARPU_USE_CUDA
+	.cuda_func = chol_cublas_codelet_update_u22,
+#endif
+#ifdef STARPU_USE_GORDON
+#ifdef SPU_FUNC_SGEMM
+	.gordon_func = SPU_FUNC_SGEMM,
+#else
+#warning SPU_FUNC_SGEMM is not available
+#endif
+#endif
+	.nbuffers = 3,
+	.model = &chol_model_22
+};
+
+static void create_task_22(unsigned k, unsigned i, unsigned j, struct starpu_sched_ctx *sched_ctx)
+{
+//	printf("task 22 k,i,j = %d,%d,%d TAG = %llx\n", k,i,j, TAG22(k,i,j));
+
+	struct starpu_task *task = create_task(TAG22(k, i, j));
+
+	task->cl = &cl22;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = ch_A_state[i][k]; 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = ch_A_state[j][k]; 
+	task->buffers[1].mode = STARPU_R;
+	task->buffers[2].handle = ch_A_state[j][i]; 
+	task->buffers[2].mode = STARPU_RW;
+
+	if ( (i == k + 1) && (j == k +1) ) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG22(k, i, j), 3, TAG22(k-1, i, j), TAG21(k, i), TAG21(k, j));
+	}
+	else {
+		starpu_tag_declare_deps(TAG22(k, i, j), 2, TAG21(k, i), TAG21(k, j));
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+
+
+/*
+ *	code to bootstrap the factorization 
+ *	and construct the DAG
+ */
+
+static double cholesky_no_stride(struct starpu_sched_ctx *sched_ctx)
+{
+	struct timeval start;
+	struct timeval end;
+
+	struct starpu_task *entry_task = NULL;
+
+	/* create all the DAG nodes */
+	unsigned i,j,k;
+
+	for (k = 0; k < nblocks; k++)
+	{
+	  struct starpu_task *task = create_task_11(k, nblocks);
+		/* we defer the launch of the first task */
+		if (k == 0) {
+			entry_task = task;
+		}
+		else {
+		  starpu_task_submit_to_ctx(task, sched_ctx);
+		}
+		
+		for (j = k+1; j<nblocks; j++)
+		{
+		  create_task_21(k, j, sched_ctx);
+
+			for (i = k+1; i<nblocks; i++)
+			{
+				if (i <= j)
+				  create_task_22(k, i, j, sched_ctx);
+			}
+		}
+	}
+
+	/* schedule the codelet */
+	gettimeofday(&start, NULL);
+		
+	starpu_task_submit_to_ctx(entry_task, sched_ctx);
+
+	/* stall the application until the end of computations */
+	starpu_tag_wait(TAG11(nblocks-1));
+	printf("cholesky finish wait for %d blocks \n", nblocks - 1);
+
+        gettimeofday(&end, NULL);
+
+        double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+        double flop = (1.0f*size*size*size)/3.0f;
+	return flop/timing/1000.0f;
+}
+
+double run_cholesky_tile_tag(struct starpu_sched_ctx *sched_ctx, int argc, char **argv)
+{
+	unsigned x, y;
+	unsigned i, j;
+
+	parse_args(argc, argv);
+	assert(nblocks <= NMAXBLOCKS);
+
+	//	fprintf(stderr, "BLOCK SIZE = %d\n", size / nblocks);
+
+	//	starpu_init(NULL);
+
+	/* Disable sequential consistency */
+	starpu_data_set_default_sequential_consistency_flag(0);
+
+	//	starpu_helper_cublas_init();
+
+	for (y = 0; y < nblocks; y++)
+	for (x = 0; x < nblocks; x++)
+	{
+		if (x <= y) {
+			ch_A[y][x] = malloc(BLOCKSIZE*BLOCKSIZE*sizeof(float));
+			assert(ch_A[y][x]);
+		}
+	}
+
+
+	for (y = 0; y < nblocks; y++)
+	for (x = 0; x < nblocks; x++)
+	{
+		if (x <= y) {
+#ifdef STARPU_HAVE_POSIX_MEMALIGN
+			posix_memalign((void **)&ch_A[y][x], 128, BLOCKSIZE*BLOCKSIZE*sizeof(float));
+#else
+			ch_A[y][x] = malloc(BLOCKSIZE*BLOCKSIZE*sizeof(float));
+#endif
+			assert(ch_A[y][x]);
+		}
+	}
+
+	/* create a simple definite positive symetric matrix example
+	 *
+	 *	Hilbert matrix : h(i,j) = 1/(i+j+1) ( + n In to make is stable ) 
+	 * */
+	for (y = 0; y < nblocks; y++)
+	for (x = 0; x < nblocks; x++)
+	if (x <= y) {
+		for (i = 0; i < BLOCKSIZE; i++)
+		for (j = 0; j < BLOCKSIZE; j++)
+		{
+			ch_A[y][x][i*BLOCKSIZE + j] =
+				(float)(1.0f/((float) (1.0+(x*BLOCKSIZE+i)+(y*BLOCKSIZE+j))));
+
+			/* make it a little more numerically stable ... ;) */
+			if ((x == y) && (i == j))
+				ch_A[y][x][i*BLOCKSIZE + j] += (float)(2*size);
+		}
+	}
+
+
+
+	for (y = 0; y < nblocks; y++)
+	for (x = 0; x < nblocks; x++)
+	{
+		if (x <= y) {
+			starpu_matrix_data_register(&ch_A_state[y][x], 0, (uintptr_t)ch_A[y][x], 
+				BLOCKSIZE, BLOCKSIZE, BLOCKSIZE, sizeof(float));
+		}
+	}
+
+	return cholesky_no_stride(sched_ctx);
+
+	//	starpu_shutdown();
+	//	return 0;
+}

examples/cholesky_2ctxs/cholesky_2ctxs.c

@@ -0,0 +1,211 @@
+#include "cholesky/cholesky.h"
+#include <pthread.h>
+
+typedef struct {
+  int start;
+  int argc;
+  char **argv;
+  unsigned ctx;
+  int the_other_ctx;
+  int *procs;
+  int ncpus;
+} params;
+
+typedef struct {
+  double flops;
+  double avg_timing;
+} retvals;
+
+#define NSAMPLES 3
+int first = 1;
+pthread_mutex_t mut;
+
+pthread_barrier_t barrier;
+
+void* func_cholesky(void *val){
+  params *p = (params*)val;
+  unsigned sched_ctx = p->ctx;
+  int the_other_ctx = p->the_other_ctx;
+
+  int i;
+  retvals *rv  = (retvals*)malloc(sizeof(retvals));
+  rv->flops = 0;
+  rv->avg_timing = 0;
+  double timing = 0;
+  for(i = 0; i < NSAMPLES; i++)
+    {
+      rv->flops += run_cholesky_implicit(sched_ctx, p->start, p->argc, p->argv, &timing, &barrier);
+      rv->avg_timing += timing;
+    }
+
+
+  pthread_mutex_lock(&mut);
+  if(first){
+      starpu_delete_sched_ctx(p->ctx, the_other_ctx);
+  }
+
+  first = 0;
+  pthread_mutex_unlock(&mut);
+ 
+
+  rv->flops /= NSAMPLES;
+  rv->avg_timing /= NSAMPLES;
+  return (void*)rv;
+}
+
+void cholesky_vs_cholesky(params *p1, params *p2, params *p3, 
+			  unsigned cpu1, unsigned cpu2,
+			  unsigned gpu, unsigned gpu1, unsigned gpu2){
+
+  int ncpus1 = cpu1 + gpu + gpu1;
+  int ncpus2 = cpu2 + gpu + gpu2;
+
+  /* 2 cholesky in different ctxs */
+  starpu_init(NULL);
+  starpu_helper_cublas_init();
+
+  int procs[ncpus1];
+  int i;
+  int k = 0;
+
+  for(i = 0; i < gpu; i++)
+    {
+      procs[k++] = i;
+      //      printf("%d ", i);
+    }
+
+  for(i = gpu; i < gpu + gpu1; i++)
+    {
+      procs[k++] = i;
+      //printf("%d ", i);
+    }
+
+  for(i = 3; i < 3 + cpu1; i++)
+    {
+      procs[k++] = i;
+      //printf("%d ", i);
+    }
+  //printf("\n");
+
+
+  p1->ctx = starpu_create_sched_ctx("heft", procs, ncpus1, "cholesky1");
+  p2->the_other_ctx = (int)p1->ctx;
+  p1->procs = procs;
+  p1->ncpus = ncpus1;
+  int procs2[ncpus2];
+
+  k = 0;
+
+  for(i = 0; i < gpu; i++){
+    procs2[k++] = i;
+    //    printf("%d ", i);
+  }
+
+  for(i = gpu + gpu1; i < gpu + gpu1 + gpu2; i++){
+    procs2[k++] = i;
+    //    printf("%d ", i);
+  }
+
+  for(i = 3  + cpu1; i < 3 + cpu1 + cpu2; i++){
+    procs2[k++] = i;
+    //    printf("%d ", i);
+  }
+
+  //  printf("\n");
+
+  p2->ctx = starpu_create_sched_ctx("prio", procs2, ncpus2, "cholesky2");
+  p1->the_other_ctx = (int)p2->ctx;
+  p2->procs = procs2;
+  p2->ncpus = ncpus2;
+
+  pthread_t tid[2];
+  pthread_barrier_init(&barrier, NULL, 2);
+  pthread_mutex_init(&mut, NULL);
+
+  struct timeval start;
+  struct timeval end;
+
+  gettimeofday(&start, NULL);
+
+
+  pthread_create(&tid[0], NULL, (void*)func_cholesky, (void*)p1);
+  pthread_create(&tid[1], NULL, (void*)func_cholesky, (void*)p2);
+
+  void *gflops_cholesky1;
+  void *gflops_cholesky2;
+ 
+  pthread_join(tid[0], &gflops_cholesky1);
+  pthread_join(tid[1], &gflops_cholesky2);
+
+  gettimeofday(&end, NULL);
+
+  pthread_mutex_destroy(&mut);
+  starpu_helper_cublas_shutdown();
+  starpu_shutdown();
+  
+  double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+  timing /= 1000000;
+  //  timing /= 60;
+
+  printf("%2.2f %2.2f ", ((retvals*)gflops_cholesky1)->flops, ((retvals*)gflops_cholesky2)->flops);
+  printf("%2.2f %2.2f %2.2f\n", ((retvals*)gflops_cholesky1)->avg_timing, ((retvals*)gflops_cholesky2)->avg_timing, timing);
+  /* printf("%2.2f %2.2f ", ((retvals*)gflops_cholesky1)->flops, 0.0 );     */
+  /*  printf("%2.2f %2.2f %2.2f\n", ((retvals*)gflops_cholesky1)->avg_timing, 0.0, timing); */
+
+}
+
+int main(int argc, char **argv)
+{
+  unsigned cpu1 = 0, cpu2 = 0;
+
+  unsigned gpu = 0, gpu1 = 0, gpu2 = 0;
+  int i;
+  
+  for (i = 9; i < argc; i++) {
+
+    if (strcmp(argv[i], "-cpu1") == 0) {
+      char *argptr;
+      cpu1 = strtol(argv[++i], &argptr, 10);
+    }    
+
+    if (strcmp(argv[i], "-cpu2") == 0) {
+      char *argptr;
+      cpu2 = strtol(argv[++i], &argptr, 10);
+    }    
+
+    if (strcmp(argv[i], "-gpu") == 0) {
+      char *argptr;
+      gpu = strtol(argv[++i], &argptr, 10);
+    }    
+
+    if (strcmp(argv[i], "-gpu1") == 0) {
+      char *argptr;
+      gpu1 = strtol(argv[++i], &argptr, 10);
+    }    
+
+    if (strcmp(argv[i], "-gpu2") == 0) {
+      char *argptr;
+      gpu2 = strtol(argv[++i], &argptr, 10);
+    }    
+
+
+  }
+
+  params p1;
+  p1.start = 1;
+  p1.argc = 5;
+  p1.argv = argv;
+
+  params p2;
+  p2.start = 5;
+  p2.argc = 9;
+  p2.argv = argv;
+
+  params p3;
+  p3.argc = argc;
+  p3.argv = argv;
+  p3.ctx = 0;
+  cholesky_vs_cholesky(&p1, &p2,&p3, cpu1, cpu2, gpu, gpu1, gpu2);
+
+  return 0;
+}