Example with cholesky and lu running in two different contexts

examples/Makefile.am

@@ -14,6 +14,8 @@
 #
 # See the GNU Lesser General Public License in COPYING.LGPL for more details.
 
+AUTOMAKE_OPTIONS = subdir-objects
+
 AM_CFLAGS = $(HWLOC_CFLAGS) -Wall
 LIBS = $(top_builddir)/src/libstarpu.la $(HWLOC_LIBS) @LIBS@
 AM_CPPFLAGS = -I$(top_srcdir)/include/ -I$(top_srcdir)/examples/ -I$(top_builddir)/include
@@ -47,6 +49,7 @@ EXTRA_DIST = 					\
 	lu/xlu_implicit_pivot.c			\
 	lu/xlu_kernels.c			\
 	lu/lu_example.c				\
+	cholesky_and_lu/cholesky_and_lu.c				\
 	incrementer/incrementer_kernels_opencl_kernel.cl 	\
 	basic_examples/variable_kernels_opencl_kernel.cl	\
 	matvecmult/matvecmult_kernel.cl				\
@@ -423,6 +426,24 @@ lu_lu_implicit_example_double_SOURCES =		\
 
 endif
 
+###########################
+# Cholesky and Lu example #
+###########################
+if !NO_BLAS_LIB
+
+examplebin_PROGRAMS += 				\
+	cholesky_and_lu/cholesky_and_lu
+
+cholesky_and_lu_cholesky_and_lu_SOURCES =			\
+	cholesky_and_lu/cholesky_and_lu.c			\
+	cholesky_and_lu/cholesky/cholesky_tile_tag.c		\
+	cholesky_and_lu/cholesky/cholesky_models.c		\
+	cholesky_and_lu/cholesky/cholesky_kernels.c		\
+	cholesky_and_lu/lu/slu.c				\
+	cholesky_and_lu/lu/slu_pivot.c				\
+	cholesky_and_lu/lu/slu_kernels.c			\
+	common/blas.c
+endif
 
 ################
 # Heat example #

examples/cholesky_and_lu/cholesky/cholesky.h

@@ -0,0 +1,121 @@
+/* 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
+
+int run_cholesky_grain_tag(struct starpu_sched_ctx *sched_ctx, int argc, char **argv);
+int run_cholesky_implicit(struct starpu_sched_ctx *sched_ctx, int argc, char **argv);
+int run_cholesky_tag(struct starpu_sched_ctx *sched_ctx, int argc, char **argv);
+int run_cholesky_tile_tag(struct starpu_sched_ctx *sched_ctx, int argc, char **argv);
+int finish_cholesky_tile_tag();
+
+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)
+{
+	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]);
+		}
+	}
+}
+
+#endif // __DW_CHOLESKY_H__

examples/cholesky_and_lu/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_and_lu/cholesky/cholesky_implicit.c

@@ -0,0 +1,283 @@
+/* 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 void _cholesky(starpu_data_handle dataA, unsigned nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	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_to_ctx(sched_ctx, &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_to_ctx(sched_ctx, &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_to_ctx(sched_ctx, &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));
+	fprintf(stderr, "Computation took (in ms)\n");
+	printf("%2.2f\n", timing/1000);
+
+	unsigned long 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 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));
+
+	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);
+}
+
+int run_cholesky_implicit(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);
+
+	//	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
+
+	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 0;
+}

examples/cholesky_and_lu/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_and_lu/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_and_lu/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_and_lu/cholesky/cholesky_tile_tag.c

@@ -0,0 +1,316 @@
+/* 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"
+
+/* A [ y ] [ x ] */
+float *A[NMAXBLOCKS][NMAXBLOCKS];
+starpu_data_handle A_state[NMAXBLOCKS][NMAXBLOCKS];
+struct timeval start;
+struct timeval end;
+
+
+/*
+ *	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|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 = 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 = A_state[k][k]; 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = 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 = A_state[i][k]; 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = A_state[j][k]; 
+	task->buffers[1].mode = STARPU_R;
+	task->buffers[2].handle = 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 void cholesky_no_stride(struct starpu_sched_ctx *sched_ctx)
+{
+	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);
+
+}
+
+int 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) {
+			A[y][x] = malloc(BLOCKSIZE*BLOCKSIZE*sizeof(float));
+			assert(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 **)&A[y][x], 128, BLOCKSIZE*BLOCKSIZE*sizeof(float));
+#else
+			A[y][x] = malloc(BLOCKSIZE*BLOCKSIZE*sizeof(float));
+#endif
+			assert(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++)
+		{
+			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))
+				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(&A_state[y][x], 0, (uintptr_t)A[y][x], 
+				BLOCKSIZE, BLOCKSIZE, BLOCKSIZE, sizeof(float));
+		}
+	}
+
+	cholesky_no_stride(sched_ctx);
+
+	//	starpu_shutdown();
+	return 0;
+}
+
+
+int finish_cholesky_tile_tag(){	
+  //	starpu_helper_cublas_shutdown();
+
+	/* stall the application until the end of computations */
+	starpu_tag_wait(TAG11(nblocks-1));
+
+	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));
+
+	return 0;
+}

examples/cholesky_and_lu/cholesky_and_lu.c

@@ -0,0 +1,27 @@
+#include "cholesky/cholesky.h"
+#include "lu/lu_example_float.c"
+
+int main(int argc, char **argv)
+{
+  starpu_init(NULL);
+
+  struct starpu_sched_ctx sched_ctx;
+  int procs[] = {1, 2, 3};
+  starpu_create_sched_ctx(&sched_ctx, "random", procs, 3);
+
+  run_cholesky_tile_tag(&sched_ctx, argc, argv);
+
+  struct starpu_sched_ctx sched_ctx2;
+  int procs2[] = {0, 4, 5, 6, 7};
+  starpu_create_sched_ctx(&sched_ctx2, "random", procs2, 5);
+
+  run_lu(&sched_ctx2, argc, argv);
+
+  finish_cholesky_tile_tag();
+  finish_lu();
+  //starpu_task_wait_for_all();
+
+  starpu_shutdown();
+  
+  return 0;
+}

examples/cholesky_and_lu/lu/dlu.c

@@ -0,0 +1,19 @@
+/* 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.
+ */
+
+#include "double.h"
+#include "xlu.c"

examples/cholesky_and_lu/lu/dlu_implicit.c

@@ -0,0 +1,19 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 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.
+ */
+
+#include "double.h"
+#include "xlu_implicit.c"

examples/cholesky_and_lu/lu/dlu_implicit_pivot.c

@@ -0,0 +1,19 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 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.
+ */
+
+#include "double.h"
+#include "xlu_implicit_pivot.c"

examples/cholesky_and_lu/lu/dlu_kernels.c

@@ -0,0 +1,19 @@
+/* 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.
+ */
+
+#include "double.h"
+#include "xlu_kernels.c"

examples/cholesky_and_lu/lu/dlu_pivot.c

@@ -0,0 +1,19 @@
+/* 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.
+ */
+
+#include "double.h"
+#include "xlu_pivot.c"

examples/cholesky_and_lu/lu/double.h

@@ -0,0 +1,46 @@
+/* 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.
+ */
+
+#define TYPE double
+
+#define STARPU_LU(name)       starpu_dlu_##name
+
+#ifdef STARPU_HAVE_MAGMA
+#include <magmablas.h>
+#define CUBLAS_GEMM	magmablas_dgemm
+#define CUBLAS_TRSM	magmablas_dtrsm
+#else
+#define CUBLAS_GEMM	cublasDgemm
+#define CUBLAS_TRSM	cublasDtrsm
+#endif
+#define CUBLAS_SCAL	cublasDscal
+#define CUBLAS_GER	cublasDger
+#define CUBLAS_SWAP	cublasDswap
+#define CUBLAS_IAMAX	cublasIdamax
+
+#define CPU_GEMM	DGEMM
+#define CPU_TRSM	DTRSM
+#define CPU_SCAL	DSCAL
+#define CPU_GER		DGER
+#define CPU_SWAP	DSWAP
+
+#define CPU_TRMM	DTRMM
+#define CPU_AXPY	DAXPY
+#define CPU_ASUM	DASUM
+#define CPU_IAMAX	IDAMAX
+
+#define PIVOT_THRESHHOLD	10e-10

examples/cholesky_and_lu/lu/float.h

@@ -0,0 +1,48 @@
+/* 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.
+ */
+
+
+#define TYPE float
+
+#define STARPU_LU(name)       starpu_slu_##name
+
+#ifdef STARPU_HAVE_MAGMA
+#include <magmablas.h>
+#define CUBLAS_GEMM	magmablas_sgemm
+#define CUBLAS_TRSM	magmablas_strsm
+#else
+#define CUBLAS_GEMM	cublasSgemm
+#define CUBLAS_TRSM	cublasStrsm
+#endif
+
+#define CUBLAS_SCAL	cublasSscal
+#define CUBLAS_GER	cublasSger
+#define CUBLAS_SWAP	cublasSswap
+#define CUBLAS_IAMAX	cublasIsamax
+
+#define CPU_GEMM	SGEMM
+#define CPU_TRSM	STRSM
+#define CPU_SCAL	SSCAL
+#define CPU_GER		SGER
+#define CPU_SWAP	SSWAP
+
+#define CPU_TRMM	STRMM
+#define CPU_AXPY	SAXPY
+#define CPU_ASUM	SASUM
+#define CPU_IAMAX	ISAMAX
+
+#define PIVOT_THRESHHOLD	10e-5

examples/cholesky_and_lu/lu/lu_example.c

@@ -0,0 +1,382 @@
+/* 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.
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+#include <starpu.h>
+#include <starpu_profiling.h>
+#include <starpu_bound.h>
+#include "xlu.h"
+#include "xlu_kernels.h"
+
+static unsigned long lu_size = 4096;
+static unsigned lu_nblocks = 16;
+static unsigned lu_check = 0;
+static unsigned pivot = 0;
+static unsigned no_stride = 0;
+static unsigned profile = 0;
+static unsigned bound = 0;
+static unsigned bounddeps = 0;
+static unsigned boundprio = 0;
+struct timeval lu_start;
+struct timeval lu_end;
+unsigned *ipiv;
+
+TYPE *A, *A_saved;
+
+/* in case we use non-strided blocks */
+TYPE **A_blocks;
+
+static void lu_parse_args(int argc, char **argv)
+{
+	int i;
+	for (i = 1; i < argc; i++) {
+		if (strcmp(argv[i], "-size") == 0) {
+			char *argptr;
+			lu_size = strtol(argv[++i], &argptr, 10);
+		}
+
+		if (strcmp(argv[i], "-nblocks") == 0) {
+			char *argptr;
+			lu_nblocks = strtol(argv[++i], &argptr, 10);
+		}
+
+		if (strcmp(argv[i], "-check") == 0) {
+			lu_check = 1;
+		}
+
+		if (strcmp(argv[i], "-piv") == 0) {
+			pivot = 1;
+		}
+
+		if (strcmp(argv[i], "-no-stride") == 0) {
+			no_stride = 1;
+		}
+
+		if (strcmp(argv[i], "-profile") == 0) {
+			profile = 1;
+		}
+
+		if (strcmp(argv[i], "-bound") == 0) {
+			bound = 1;
+		}
+		if (strcmp(argv[i], "-bounddeps") == 0) {
+			bound = 1;
+			bounddeps = 1;
+		}
+		if (strcmp(argv[i], "-bounddepsprio") == 0) {
+			bound = 1;
+			bounddeps = 1;
+			boundprio = 1;
+		}
+	}
+}
+
+static void display_matrix(TYPE *m, unsigned n, unsigned ld, char *str)
+{
+#if 0
+	fprintf(stderr, "***********\n");
+	fprintf(stderr, "Display matrix %s\n", str);
+	unsigned i,j;
+	for (j = 0; j < n; j++)
+	{
+		for (i = 0; i < n; i++)
+		{
+			fprintf(stderr, "%2.2f\t", m[i+j*ld]);
+		}
+		fprintf(stderr, "\n");
+	}
+	fprintf(stderr, "***********\n");
+#endif
+}
+
+void copy_blocks_into_matrix(void)
+{
+	unsigned blocklu_size = (lu_size/lu_nblocks);
+
+	unsigned i, j;
+	unsigned bi, bj;
+	for (bj = 0; bj < lu_nblocks; bj++)
+	for (bi = 0; bi < lu_nblocks; bi++)
+	{
+		for (j = 0; j < blocklu_size; j++)
+		for (i = 0; i < blocklu_size; i++)
+		{
+			A[(i+bi*blocklu_size) + (j + bj*blocklu_size)*lu_size] =
+				A_blocks[bi+lu_nblocks*bj][i + j * blocklu_size];
+		}
+
+		//free(A_blocks[bi+nblocks*bj]);
+	}
+}
+
+
+
+void copy_matrix_into_blocks(void)
+{
+	unsigned blocklu_size = (lu_size/lu_nblocks);
+
+	unsigned i, j;
+	unsigned bi, bj;
+	for (bj = 0; bj < lu_nblocks; bj++)
+	for (bi = 0; bi < lu_nblocks; bi++)
+	{
+		starpu_data_malloc_pinned_if_possible((void **)&A_blocks[bi+lu_nblocks*bj], (size_t)blocklu_size*blocklu_size*sizeof(TYPE));
+
+		for (j = 0; j < blocklu_size; j++)
+		for (i = 0; i < blocklu_size; i++)
+		{
+			A_blocks[bi+lu_nblocks*bj][i + j * blocklu_size] =
+			A[(i+bi*blocklu_size) + (j + bj*blocklu_size)*lu_size];
+		}
+	}
+}
+
+static void init_matrix(void)
+{
+	/* allocate matrix */
+	starpu_data_malloc_pinned_if_possible((void **)&A, (size_t)lu_size*lu_size*sizeof(TYPE));
+	STARPU_ASSERT(A);
+
+	starpu_srand48((long int)time(NULL));
+	//starpu_srand48(0);
+
+	/* initialize matrix content */
+	unsigned long i,j;
+	for (j = 0; j < lu_size; j++)
+	{
+		for (i = 0; i < lu_size; i++)
+		{
+			A[i + j*lu_size] = (TYPE)starpu_drand48();
+		}
+	}
+
+}
+
+static void save_matrix(void)
+{
+	A_saved = malloc((size_t)lu_size*lu_size*sizeof(TYPE));
+	STARPU_ASSERT(A_saved);
+
+	memcpy(A_saved, A, (size_t)lu_size*lu_size*sizeof(TYPE));
+}
+
+static double frobenius_norm(TYPE *v, unsigned n)
+{
+	double sum2 = 0.0;
+
+	/* compute sqrt(Sum(|x|^2)) */
+
+	unsigned i,j;
+	for (j = 0; j < n; j++)
+	for (i = 0; i < n; i++)
+	{
+		double a = fabsl((double)v[i+n*j]);
+		sum2 += a*a;
+	}
+
+	return sqrt(sum2);
+}
+
+static void pivot_saved_matrix(unsigned *ipiv)
+{
+	unsigned k;
+	for (k = 0; k < lu_size; k++)
+	{
+		if (k != ipiv[k])
+		{
+	//		fprintf(stderr, "SWAP %d and %d\n", k, ipiv[k]);
+			CPU_SWAP(lu_size, &A_saved[k*lu_size], 1, &A_saved[ipiv[k]*lu_size], 1);
+		}
+	}
+}
+
+static void check_result(void)
+{
+	unsigned i,j;
+	TYPE *L, *U;
+
+	L = malloc((size_t)lu_size*lu_size*sizeof(TYPE));
+	U = malloc((size_t)lu_size*lu_size*sizeof(TYPE));
+
+	memset(L, 0, lu_size*lu_size*sizeof(TYPE));
+	memset(U, 0, lu_size*lu_size*sizeof(TYPE));
+
+	/* only keep the lower part */
+	for (j = 0; j < lu_size; j++)
+	{
+		for (i = 0; i < j; i++)
+		{
+			L[j+i*lu_size] = A[j+i*lu_size];
+		}
+
+		/* diag i = j */
+		L[j+j*lu_size] = A[j+j*lu_size];
+		U[j+j*lu_size] = 1.0;
+
+		for (i = j+1; i < lu_size; i++)
+		{
+			U[j+i*lu_size] = A[j+i*lu_size];
+		}
+	}
+
+	display_matrix(L, lu_size, lu_size, "L");
+	display_matrix(U, lu_size, lu_size, "U");
+
+	/* now A_err = L, compute L*U */
+	CPU_TRMM("R", "U", "N", "U", lu_size, lu_size, 1.0f, U, lu_size, L, lu_size);
+
+	display_matrix(A_saved, lu_size, lu_size, "P A_saved");
+	display_matrix(L, lu_size, lu_size, "LU");
+
+	/* compute "LU - A" in L*/
+	CPU_AXPY(lu_size*lu_size, -1.0, A_saved, 1, L, 1);
+	display_matrix(L, lu_size, lu_size, "Residuals");
+	
+	TYPE err = CPU_ASUM(lu_size*lu_size, L, 1);
+	int max = CPU_IAMAX(lu_size*lu_size, L, 1);
+
+	fprintf(stderr, "Avg error : %e\n", err/(lu_size*lu_size));
+	fprintf(stderr, "Max error : %e\n", L[max]);
+
+	double residual = frobenius_norm(L, lu_size);
+	double matnorm = frobenius_norm(A_saved, lu_size);
+
+	fprintf(stderr, "||%sA-LU|| / (||A||*N) : %e\n", pivot?"P":"", residual/(matnorm*lu_size));
+
+	if (residual/(matnorm*lu_size) > 1e-5)
+		exit(-1);
+}
+
+int run_lu(struct starpu_sched_ctx *sched_ctx, int argc, char **argv)
+{
+	lu_parse_args(argc, argv);
+
+	//	starpu_init(NULL);
+
+	//	starpu_helper_cublas_init();
+
+	init_matrix();
+
+	if (lu_check)
+		save_matrix();
+
+	display_matrix(A, lu_size, lu_size, "A");
+
+	if (bound)
+		starpu_bound_start(bounddeps, boundprio);
+
+	if (profile)
+		starpu_profiling_status_set(STARPU_PROFILING_ENABLE);
+
+	/* Factorize the matrix (in place) */
+	if (pivot)
+	{
+ 		ipiv = malloc(lu_size*sizeof(unsigned));
+		if (no_stride)
+		{
+			/* in case the LU decomposition uses non-strided blocks, we _copy_ the matrix into smaller blocks */
+			A_blocks = malloc(lu_nblocks*lu_nblocks*sizeof(TYPE **));
+			copy_matrix_into_blocks();
+
+			STARPU_LU(lu_decomposition_pivot_no_stride)(A_blocks, ipiv, lu_size, lu_size, lu_nblocks, sched_ctx);
+
+			copy_blocks_into_matrix();
+			free(A_blocks);
+		}
+		else 
+		{
+			gettimeofday(&lu_start, NULL);
+
+			STARPU_LU(lu_decomposition_pivot)(A, ipiv, lu_size, lu_size, lu_nblocks, sched_ctx);
+		}
+	}
+	else
+	{
+	  STARPU_LU(lu_decomposition)(A, lu_size, lu_size, lu_nblocks, sched_ctx);
+	}
+
+
+	//	starpu_shutdown();
+
+	return 0;
+}
+
+int finish_lu()
+{
+
+	if (pivot)
+	{
+		if (no_stride)
+		{
+			finish_lu_decomposition_pivot_no_stride(lu_nblocks);
+		}
+		else 
+		{
+			finish_lu_decomposition_pivot(lu_nblocks);
+			gettimeofday(&lu_end, NULL);
+			
+			double timing = (double)((lu_end.tv_sec - lu_start.tv_sec)*1000000 + (lu_end.tv_usec - lu_start.tv_usec));
+			
+			unsigned n = lu_size;
+			double flop = (2.0f*n*n*n)/3.0f;
+			fprintf(stderr, "Synthetic GFlops (TOTAL) : \n");
+			fprintf(stdout, "%d	%6.2f\n", n, (flop/timing/1000.0f));
+
+		}
+	}
+	else
+	{
+	  finish_lu_decomposition(lu_nblocks);
+	}
+
+  
+  
+	if (profile)
+	{
+		starpu_profiling_status_set(STARPU_PROFILING_DISABLE);
+		starpu_bus_profiling_helper_display_summary();
+	}
+
+	if (bound) {
+		double min;
+		starpu_bound_stop();
+		if (bounddeps) {
+			FILE *f = fopen("lu.pl", "w");
+			starpu_bound_print_lp(f);
+			fprintf(stderr,"system printed to lu.pl\n");
+		} else {
+			starpu_bound_compute(&min, NULL, 0);
+			if (min != 0.)
+				fprintf(stderr, "theoretical min: %lf ms\n", min);
+		}
+	}
+
+	if (lu_check)
+	{
+		if (pivot)
+			pivot_saved_matrix(ipiv);
+
+		check_result();
+	}
+
+	starpu_helper_cublas_shutdown();
+	return 0;
+}

examples/cholesky_and_lu/lu/lu_example_double.c

@@ -0,0 +1,19 @@
+/* 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.
+ */
+
+#include "double.h"
+#include "lu_example.c"

examples/cholesky_and_lu/lu/lu_example_float.c

@@ -0,0 +1,19 @@
+/* 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.
+ */
+
+#include "float.h"
+#include "lu_example.c"

examples/cholesky_and_lu/lu/slu.c

@@ -0,0 +1,19 @@
+/* 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.
+ */
+
+#include "float.h"
+#include "xlu.c"

examples/cholesky_and_lu/lu/slu_implicit.c

@@ -0,0 +1,19 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 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.
+ */
+
+#include "float.h"
+#include "xlu_implicit.c"

examples/cholesky_and_lu/lu/slu_implicit_pivot.c

@@ -0,0 +1,19 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 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.
+ */
+
+#include "float.h"
+#include "xlu_implicit_pivot.c"

examples/cholesky_and_lu/lu/slu_kernels.c

@@ -0,0 +1,19 @@
+/* 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.
+ */
+
+#include "float.h"
+#include "xlu_kernels.c"

examples/cholesky_and_lu/lu/slu_pivot.c

@@ -0,0 +1,19 @@
+/* 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.
+ */
+
+#include "float.h"
+#include "xlu_pivot.c"

examples/cholesky_and_lu/lu/xlu.c

@@ -0,0 +1,256 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010  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 "xlu.h"
+#include "xlu_kernels.h"
+
+#define TAG11(k)	((starpu_tag_t)( (1ULL<<60) | (unsigned long long)(k)))
+#define TAG12(k,i)	((starpu_tag_t)(((2ULL<<60) | (((unsigned long long)(k))<<32)	\
+					| (unsigned long long)(i))))
+#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))))
+
+static unsigned no_prio = 0;
+struct timeval xlu_start;
+struct timeval xlu_end;
+starpu_data_handle xlu_dataA;
+
+
+
+/*
+ *	Construct the DAG
+ */
+
+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;
+}
+
+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 (!no_prio)
+		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 void create_task_12(starpu_data_handle dataA, unsigned k, unsigned j, struct starpu_sched_ctx *sched_ctx)
+{
+//	printf("task 12 k,i = %d,%d TAG = %llx\n", k,i, TAG12(k,i));
+
+	struct starpu_task *task = create_task(TAG12(k, j));
+	
+	task->cl = &cl12;
+
+	/* 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, j, k); 
+	task->buffers[1].mode = STARPU_RW;
+
+	if (!no_prio && (j == k+1)) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG12(k, j), 2, TAG11(k), TAG22(k-1, k, j));
+	}
+	else {
+		starpu_tag_declare_deps(TAG12(k, j), 1, TAG11(k));
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_21(starpu_data_handle dataA, unsigned k, unsigned i, struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = create_task(TAG21(k, i));
+
+	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, i); 
+	task->buffers[1].mode = STARPU_RW;
+
+	if (!no_prio && (i == k+1)) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	/* enforce dependencies ... */
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG21(k, i), 2, TAG11(k), TAG22(k-1, i, k));
+	}
+	else {
+		starpu_tag_declare_deps(TAG21(k, i), 1, TAG11(k));
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+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); /* produced by TAG21(k, i) */ 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = starpu_data_get_sub_data(dataA, 2, j, k); /* produced by TAG12(k, j) */
+	task->buffers[1].mode = STARPU_R;
+	task->buffers[2].handle = starpu_data_get_sub_data(dataA, 2, j, i); /* produced by TAG22(k-1, i, j) */
+	task->buffers[2].mode = STARPU_RW;
+
+	if (!no_prio &&  (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), TAG12(k, j), TAG21(k, i));
+	}
+	else {
+		starpu_tag_declare_deps(TAG22(k, i, j), 2, TAG12(k, j), TAG21(k, i));
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+/*
+ *	code to bootstrap the factorization 
+ */
+
+static void dw_codelet_facto_v3(starpu_data_handle dataA, unsigned nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	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(dataA, k);
+
+		/* we defer the launch of the first task */
+		if (k == 0) {
+			entry_task = task;
+		}
+		else {
+		  starpu_task_submit_to_ctx(task, sched_ctx);
+		}
+		
+		for (i = k+1; i<nblocks; i++)
+		{
+			create_task_12(dataA, k, i, sched_ctx);
+			create_task_21(dataA, k, i, sched_ctx);
+		}
+
+		for (i = k+1; i<nblocks; i++)
+		{
+			for (j = k+1; j<nblocks; j++)
+			{
+			  create_task_22(dataA, k, i, j, sched_ctx);
+			}
+		}
+	}
+
+	/* schedule the codelet */
+	gettimeofday(&xlu_start, NULL);
+	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);
+	}
+
+}
+
+void STARPU_LU(lu_decomposition)(TYPE *matA, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	/* monitor and partition the A matrix into blocks :
+	 * one block is now determined by 2 unsigned (i,j) */
+	starpu_matrix_data_register(&xlu_dataA, 0, (uintptr_t)matA, ld, size, size, sizeof(TYPE));
+
+	/* We already enforce deps by hand */
+	starpu_data_set_sequential_consistency_flag(xlu_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(xlu_dataA, 2, &f, &f2);
+
+	dw_codelet_facto_v3(xlu_dataA, nblocks, sched_ctx);
+}
+
+void finish_lu_decomposition(unsigned nblocks)
+{
+	/* stall the application until the end of computations */
+	starpu_tag_wait(TAG11(nblocks-1));
+
+	gettimeofday(&xlu_end, NULL);
+
+	double timing = (double)((xlu_end.tv_sec - xlu_start.tv_sec)*1000000 + (xlu_end.tv_usec - xlu_start.tv_usec));
+	fprintf(stderr, "Computation took (in ms)\n");
+	printf("%2.2f\n", timing/1000);
+
+	unsigned n = starpu_matrix_get_nx(xlu_dataA);
+	double flop = (2.0f*n*n*n)/3.0f;
+	fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
+
+	/* gather all the data */
+	starpu_data_unpartition(xlu_dataA, 0);
+}

examples/cholesky_and_lu/lu/xlu.h

@@ -0,0 +1,121 @@
+/* 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.
+ */
+
+#ifndef __XLU_H__
+#define __XLU_H__
+
+#include <sys/time.h>
+
+/* for STARPU_USE_CUDA */
+#include <starpu_config.h>
+#include <starpu.h>
+#include <starpu_cuda.h>
+
+#include <common/blas.h>
+
+#define BLAS3_FLOP(n1,n2,n3)    \
+        (2*((uint64_t)n1)*((uint64_t)n2)*((uint64_t)n3))
+
+#ifdef CHECK_RESULTS
+static void __attribute__ ((unused)) compare_A_LU(float *A, float *LU,
+				unsigned size, unsigned ld)
+{
+	unsigned i,j;
+	float *L;
+	float *U;
+
+	L = malloc(size*size*sizeof(float));
+	U = malloc(size*size*sizeof(float));
+
+	memset(L, 0, size*size*sizeof(float));
+	memset(U, 0, size*size*sizeof(float));
+
+	/* only keep the lower part */
+	for (j = 0; j < size; j++)
+	{
+		for (i = 0; i < j; i++)
+		{
+			L[j+i*size] = LU[j+i*ld];
+		}
+
+		/* diag i = j */
+		L[j+j*size] = LU[j+j*ld];
+		U[j+j*size] = 1.0f;
+
+		for (i = j+1; i < size; i++)
+		{
+			U[j+i*size] = LU[j+i*ld];
+		}
+	}
+
+        /* now A_err = L, compute L*U */
+	STRMM("R", "U", "N", "U", size, size, 1.0f, U, size, L, size);
+
+	float max_err = 0.0f;
+	for (i = 0; i < size ; i++)
+	{
+		for (j = 0; j < size; j++) 
+		{
+			max_err = STARPU_MAX(max_err, fabs(  L[j+i*size] - A[j+i*ld]  ));
+		}
+	}
+
+	printf("max error between A and L*U = %f \n", max_err);
+}
+#endif // CHECK_RESULTS
+
+void dw_cpu_codelet_update_u11(void **, void *);
+void dw_cpu_codelet_update_u12(void **, void *);
+void dw_cpu_codelet_update_u21(void **, void *);
+void dw_cpu_codelet_update_u22(void **, void *);
+
+#ifdef STARPU_USE_CUDA
+void dw_cublas_codelet_update_u11(void *descr[], void *_args);
+void dw_cublas_codelet_update_u12(void *descr[], void *_args);
+void dw_cublas_codelet_update_u21(void *descr[], void *_args);
+void dw_cublas_codelet_update_u22(void *descr[], void *_args);
+#endif
+
+void dw_callback_codelet_update_u11(void *);
+void dw_callback_codelet_update_u12_21(void *);
+void dw_callback_codelet_update_u22(void *);
+
+void dw_callback_v2_codelet_update_u11(void *);
+void dw_callback_v2_codelet_update_u12(void *);
+void dw_callback_v2_codelet_update_u21(void *);
+void dw_callback_v2_codelet_update_u22(void *);
+
+extern struct starpu_perfmodel_t model_11;
+extern struct starpu_perfmodel_t model_12;
+extern struct starpu_perfmodel_t model_21;
+extern struct starpu_perfmodel_t model_22;
+
+struct piv_s {
+	unsigned *piv; /* complete pivot array */
+	unsigned first; /* first element */
+	unsigned last; /* last element */
+};
+
+void STARPU_LU(lu_decomposition)(TYPE *matA, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx);
+void STARPU_LU(lu_decomposition_pivot_no_stride)(TYPE **matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx);
+void STARPU_LU(lu_decomposition_pivot)(TYPE *matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx);
+
+void finish_lu_decomposition(unsigned nblocks);
+void finish_lu_decomposition_pivot_no_stride(unsigned nblocks);
+void finish_lu_decomposition_pivot(unsigned nblocks);
+
+#endif // __XLU_H__

examples/cholesky_and_lu/lu/xlu_implicit.c

@@ -0,0 +1,164 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2010  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 "xlu.h"
+#include "xlu_kernels.h"
+
+static unsigned no_prio = 0;
+
+static void create_task_11(starpu_data_handle dataA, unsigned k, struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = starpu_task_create();
+	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 (!no_prio)
+		task->priority = STARPU_MAX_PRIO;
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_12(starpu_data_handle dataA, unsigned k, unsigned j, struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = starpu_task_create();
+	task->cl = &cl12;
+
+	/* 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, j, k); 
+	task->buffers[1].mode = STARPU_RW;
+
+	if (!no_prio && (j == k+1))
+		task->priority = STARPU_MAX_PRIO;
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_21(starpu_data_handle dataA, unsigned k, unsigned i, struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = starpu_task_create();
+
+	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, i); 
+	task->buffers[1].mode = STARPU_RW;
+
+	if (!no_prio && (i == k+1))
+		task->priority = STARPU_MAX_PRIO;
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_22(starpu_data_handle dataA, unsigned k, unsigned i, unsigned j, struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = starpu_task_create();
+
+	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, j, k);
+	task->buffers[1].mode = STARPU_R;
+	task->buffers[2].handle = starpu_data_get_sub_data(dataA, 2, j, i);
+	task->buffers[2].mode = STARPU_RW;
+
+	if (!no_prio &&  (i == k + 1) && (j == k +1) )
+		task->priority = STARPU_MAX_PRIO;
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+/*
+ *	code to bootstrap the factorization 
+ */
+
+static void dw_codelet_facto_v3(starpu_data_handle dataA, unsigned nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	struct timeval start;
+	struct timeval end;
+
+	/* create all the DAG nodes */
+	unsigned i,j,k;
+
+	gettimeofday(&start, NULL);
+
+	for (k = 0; k < nblocks; k++)
+	{
+	  create_task_11(dataA, k, sched_ctx);
+		
+		for (i = k+1; i<nblocks; i++)
+		{
+			create_task_12(dataA, k, i, sched_ctx);
+			create_task_21(dataA, k, i, sched_ctx);
+		}
+
+		for (i = k+1; i<nblocks; i++)
+			for (j = k+1; j<nblocks; j++)
+			  create_task_22(dataA, k, i, j, sched_ctx);
+	}
+
+	/* stall the application until the end of computations */
+	starpu_task_wait_for_all();
+
+	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 = (2.0f*n*n*n)/3.0f;
+	fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
+}
+
+void STARPU_LU(lu_decomposition)(TYPE *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(TYPE));
+	
+	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);
+
+	dw_codelet_facto_v3(dataA, nblocks, sched_ctx);
+
+	/* gather all the data */
+	starpu_data_unpartition(dataA, 0);
+}

examples/cholesky_and_lu/lu/xlu_implicit_pivot.c

@@ -0,0 +1,283 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2010  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 "xlu.h"
+#include "xlu_kernels.h"
+
+static unsigned no_prio = 0;
+
+/*
+ *	Construct the DAG
+ */
+
+static void create_task_pivot(starpu_data_handle *dataAp, unsigned nblocks,
+			      struct piv_s *piv_description,
+			      unsigned k, unsigned i,
+			      starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = starpu_task_create();
+
+	task->cl = &cl_pivot;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, i);
+	task->buffers[0].mode = STARPU_RW;
+
+	task->cl_arg = &piv_description[k];
+
+	/* this is an important task */
+	if (!no_prio && (i == k+1))
+		task->priority = STARPU_MAX_PRIO;
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_11_pivot(starpu_data_handle *dataAp, unsigned nblocks,
+				 unsigned k, struct piv_s *piv_description,
+				 starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = starpu_task_create();
+
+	task->cl = &cl11_pivot;
+
+	task->cl_arg = &piv_description[k];
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, k);
+	task->buffers[0].mode = STARPU_RW;
+
+	/* this is an important task */
+	if (!no_prio)
+		task->priority = STARPU_MAX_PRIO;
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_12(starpu_data_handle *dataAp, unsigned nblocks, unsigned k, unsigned j,
+		starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = starpu_task_create();
+	
+	task->cl = &cl12;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, k);
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = get_block(dataAp, nblocks, j, k);
+	task->buffers[1].mode = STARPU_RW;
+
+	if (!no_prio && (j == k+1))
+		task->priority = STARPU_MAX_PRIO;
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_21(starpu_data_handle *dataAp, unsigned nblocks, unsigned k, unsigned i,
+				starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = starpu_task_create();
+
+	task->cl = &cl21;
+	
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, k); 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = get_block(dataAp, nblocks, k, i); 
+	task->buffers[1].mode = STARPU_RW;
+
+	if (!no_prio && (i == k+1))
+		task->priority = STARPU_MAX_PRIO;
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_22(starpu_data_handle *dataAp, unsigned nblocks, unsigned k, unsigned i, unsigned j,
+				starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = starpu_task_create();
+
+	task->cl = &cl22;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, i);
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = get_block(dataAp, nblocks, j, k);
+	task->buffers[1].mode = STARPU_R;
+	task->buffers[2].handle = get_block(dataAp, nblocks, j, i);
+	task->buffers[2].mode = STARPU_RW;
+
+	if (!no_prio &&  (i == k + 1) && (j == k +1) )
+		task->priority = STARPU_MAX_PRIO;
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+/*
+ *	code to bootstrap the factorization 
+ */
+
+static double dw_codelet_facto_pivot(starpu_data_handle *dataAp,
+					struct piv_s *piv_description,
+					unsigned nblocks,
+					starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx)
+{
+	struct timeval start;
+	struct timeval end;
+
+	gettimeofday(&start, NULL);
+
+	/* create all the DAG nodes */
+	unsigned i,j,k;
+	for (k = 0; k < nblocks; k++)
+	{
+	  create_task_11_pivot(dataAp, nblocks, k, piv_description, get_block, sched_ctx);
+
+		for (i = 0; i < nblocks; i++)
+		{
+			if (i != k)
+			  create_task_pivot(dataAp, nblocks, piv_description, k, i, get_block, sched_ctx);
+		}
+	
+		for (i = k+1; i<nblocks; i++)
+		{
+		  create_task_12(dataAp, nblocks, k, i, get_block, sched_ctx);
+		  create_task_21(dataAp, nblocks, k, i, get_block, sched_ctx);
+		}
+
+		for (i = k+1; i<nblocks; i++)
+		for (j = k+1; j<nblocks; j++)
+		  create_task_22(dataAp, nblocks, k, i, j, get_block, sched_ctx);
+	}
+
+	/* stall the application until the end of computations */
+	starpu_task_wait_for_all();
+
+	gettimeofday(&end, NULL);
+
+	double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+	return timing;
+}
+
+starpu_data_handle get_block_with_striding(starpu_data_handle *dataAp,
+			unsigned nblocks __attribute__((unused)), unsigned j, unsigned i)
+{
+	/* we use filters */
+	return starpu_data_get_sub_data(*dataAp, 2, j, i);
+}
+
+
+void STARPU_LU(lu_decomposition_pivot)(TYPE *matA, unsigned *ipiv, 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(TYPE));
+
+	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);
+
+	unsigned i;
+	for (i = 0; i < size; i++)
+		ipiv[i] = i;
+
+	struct piv_s *piv_description = malloc(nblocks*sizeof(struct piv_s));
+	unsigned block;
+	for (block = 0; block < nblocks; block++)
+	{
+		piv_description[block].piv = ipiv;
+		piv_description[block].first = block * (size / nblocks);
+		piv_description[block].last = (block + 1) * (size / nblocks);
+	}
+
+	double timing;
+	timing = dw_codelet_facto_pivot(&dataA, piv_description, nblocks, get_block_with_striding, sched_ctx);
+
+	fprintf(stderr, "Computation took (in ms)\n");
+	fprintf(stderr, "%2.2f\n", timing/1000);
+
+	unsigned n = starpu_matrix_get_nx(dataA);
+	double flop = (2.0f*n*n*n)/3.0f;
+	fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
+
+	/* gather all the data */
+	starpu_data_unpartition(dataA, 0);
+}
+
+
+starpu_data_handle get_block_with_no_striding(starpu_data_handle *dataAp, unsigned nblocks, unsigned j, unsigned i)
+{
+	/* dataAp is an array of data handle */
+	return dataAp[i+j*nblocks];
+}
+
+void STARPU_LU(lu_decomposition_pivot_no_stride)(TYPE **matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned nblocks)
+{
+	starpu_data_handle *dataAp = malloc(nblocks*nblocks*sizeof(starpu_data_handle));
+
+	/* monitor and partition the A matrix into blocks :
+	 * one block is now determined by 2 unsigned (i,j) */
+	unsigned bi, bj;
+	for (bj = 0; bj < nblocks; bj++)
+	for (bi = 0; bi < nblocks; bi++)
+	{
+		starpu_matrix_data_register(&dataAp[bi+nblocks*bj], 0,
+			(uintptr_t)matA[bi+nblocks*bj], size/nblocks,
+			size/nblocks, size/nblocks, sizeof(TYPE));
+	}
+
+	unsigned i;
+	for (i = 0; i < size; i++)
+		ipiv[i] = i;
+
+	struct piv_s *piv_description = malloc(nblocks*sizeof(struct piv_s));
+	unsigned block;
+	for (block = 0; block < nblocks; block++)
+	{
+		piv_description[block].piv = ipiv;
+		piv_description[block].first = block * (size / nblocks);
+		piv_description[block].last = (block + 1) * (size / nblocks);
+	}
+
+	double timing;
+	timing = dw_codelet_facto_pivot(dataAp, piv_description, nblocks, get_block_with_no_striding, sched_ctx);
+
+	fprintf(stderr, "Computation took (in ms)\n");
+	fprintf(stderr, "%2.2f\n", timing/1000);
+
+	unsigned n = starpu_matrix_get_nx(dataAp[0])*nblocks;
+	double flop = (2.0f*n*n*n)/3.0f;
+	fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
+
+	for (bj = 0; bj < nblocks; bj++)
+	for (bi = 0; bi < nblocks; bi++)
+	{
+		starpu_data_unregister(dataAp[bi+nblocks*bj]);
+	}
+}

examples/cholesky_and_lu/lu/xlu_kernels.c

@@ -0,0 +1,584 @@
+/* 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.
+ */
+
+#include "xlu.h"
+#include <math.h>
+
+#define str(s) #s
+#define xstr(s)        str(s)
+#define STARPU_LU_STR(name)  xstr(STARPU_LU(name))
+
+/*
+ *   U22 
+ */
+
+static inline void STARPU_LU(common_u22)(void *descr[],
+				int s, __attribute__((unused)) void *_args)
+{
+	TYPE *right 	= (TYPE *)STARPU_MATRIX_GET_PTR(descr[0]);
+	TYPE *left 	= (TYPE *)STARPU_MATRIX_GET_PTR(descr[1]);
+	TYPE *center 	= (TYPE *)STARPU_MATRIX_GET_PTR(descr[2]);
+
+	unsigned dx = STARPU_MATRIX_GET_NX(descr[2]);
+	unsigned dy = STARPU_MATRIX_GET_NY(descr[2]);
+	unsigned dz = STARPU_MATRIX_GET_NY(descr[0]);
+
+	unsigned ld12 = STARPU_MATRIX_GET_LD(descr[0]);
+	unsigned ld21 = STARPU_MATRIX_GET_LD(descr[1]);
+	unsigned ld22 = STARPU_MATRIX_GET_LD(descr[2]);
+
+#ifdef STARPU_USE_CUDA
+	cublasStatus status;
+	cudaError_t cures;
+#endif
+
+	switch (s) {
+		case 0:
+			CPU_GEMM("N", "N", dy, dx, dz, 
+				(TYPE)-1.0, right, ld21, left, ld12,
+				(TYPE)1.0, center, ld22);
+			break;
+
+#ifdef STARPU_USE_CUDA
+		case 1:
+			CUBLAS_GEMM('n', 'n', dx, dy, dz,
+				(TYPE)-1.0, right, ld21, left, ld12,
+				(TYPE)1.0f, center, ld22);
+
+			status = cublasGetError();
+			if (STARPU_UNLIKELY(status != CUBLAS_STATUS_SUCCESS))
+				STARPU_ABORT();
+
+			if (STARPU_UNLIKELY((cures = cudaThreadSynchronize()) != cudaSuccess))
+				STARPU_CUDA_REPORT_ERROR(cures);
+
+			break;
+#endif
+		default:
+			STARPU_ABORT();
+			break;
+	}
+}
+
+void STARPU_LU(cpu_u22)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u22)(descr, 0, _args);
+}
+
+#ifdef STARPU_USE_CUDA
+void STARPU_LU(cublas_u22)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u22)(descr, 1, _args);
+}
+#endif// STARPU_USE_CUDA
+
+static struct starpu_perfmodel_t STARPU_LU(model_22) = {
+	.type = STARPU_HISTORY_BASED,
+#ifdef STARPU_ATLAS
+	.symbol = STARPU_LU_STR(lu_model_22_atlas)
+#elif defined(STARPU_GOTO)
+	.symbol = STARPU_LU_STR(lu_model_22_goto)
+#else
+	.symbol = STARPU_LU_STR(lu_model_22)
+#endif
+};
+
+starpu_codelet cl22 = {
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = STARPU_LU(cpu_u22),
+#ifdef STARPU_USE_CUDA
+	.cuda_func = STARPU_LU(cublas_u22),
+#endif
+	.nbuffers = 3,
+	.model = &STARPU_LU(model_22)
+};
+
+/*
+ * U12
+ */
+
+static inline void STARPU_LU(common_u12)(void *descr[],
+				int s, __attribute__((unused)) void *_args)
+{
+	TYPE *sub11;
+	TYPE *sub12;
+
+	sub11 = (TYPE *)STARPU_MATRIX_GET_PTR(descr[0]);	
+	sub12 = (TYPE *)STARPU_MATRIX_GET_PTR(descr[1]);
+
+	unsigned ld11 = STARPU_MATRIX_GET_LD(descr[0]);
+	unsigned ld12 = STARPU_MATRIX_GET_LD(descr[1]);
+
+	unsigned nx12 = STARPU_MATRIX_GET_NX(descr[1]);
+	unsigned ny12 = STARPU_MATRIX_GET_NY(descr[1]);
+
+#ifdef STARPU_USE_CUDA
+	cublasStatus status;
+	cudaError_t cures;
+#endif
+
+	/* solve L11 U12 = A12 (find U12) */
+	switch (s) {
+		case 0:
+			CPU_TRSM("L", "L", "N", "N", nx12, ny12,
+					(TYPE)1.0, sub11, ld11, sub12, ld12);
+			break;
+#ifdef STARPU_USE_CUDA
+		case 1:
+			CUBLAS_TRSM('L', 'L', 'N', 'N', ny12, nx12,
+					(TYPE)1.0, sub11, ld11, sub12, ld12);
+
+			status = cublasGetError();
+			if (STARPU_UNLIKELY(status != CUBLAS_STATUS_SUCCESS))
+				STARPU_ABORT();
+
+			if (STARPU_UNLIKELY((cures = cudaThreadSynchronize()) != cudaSuccess))
+				STARPU_CUDA_REPORT_ERROR(cures);
+
+			break;
+#endif
+		default:
+			STARPU_ABORT();
+			break;
+	}
+}
+
+void STARPU_LU(cpu_u12)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u12)(descr, 0, _args);
+}
+
+#ifdef STARPU_USE_CUDA
+void STARPU_LU(cublas_u12)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u12)(descr, 1, _args);
+}
+#endif // STARPU_USE_CUDA
+
+static struct starpu_perfmodel_t STARPU_LU(model_12) = {
+	.type = STARPU_HISTORY_BASED,
+#ifdef STARPU_ATLAS
+	.symbol = STARPU_LU_STR(lu_model_12_atlas)
+#elif defined(STARPU_GOTO)
+	.symbol = STARPU_LU_STR(lu_model_12_goto)
+#else
+	.symbol = STARPU_LU_STR(lu_model_12)
+#endif
+};
+
+starpu_codelet cl12 = {
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = STARPU_LU(cpu_u12),
+#ifdef STARPU_USE_CUDA
+	.cuda_func = STARPU_LU(cublas_u12),
+#endif
+	.nbuffers = 2,
+	.model = &STARPU_LU(model_12)
+};
+
+/* 
+ * U21
+ */
+
+static inline void STARPU_LU(common_u21)(void *descr[],
+				int s, __attribute__((unused)) void *_args)
+{
+	TYPE *sub11;
+	TYPE *sub21;
+
+	sub11 = (TYPE *)STARPU_MATRIX_GET_PTR(descr[0]);
+	sub21 = (TYPE *)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_NX(descr[1]);
+	unsigned ny21 = STARPU_MATRIX_GET_NY(descr[1]);
+	
+#ifdef STARPU_USE_CUDA
+	cublasStatus status;
+	cudaError_t cures;
+#endif
+
+	switch (s) {
+		case 0:
+			CPU_TRSM("R", "U", "N", "U", nx21, ny21,
+					(TYPE)1.0, sub11, ld11, sub21, ld21);
+			break;
+#ifdef STARPU_USE_CUDA
+		case 1:
+			CUBLAS_TRSM('R', 'U', 'N', 'U', ny21, nx21,
+					(TYPE)1.0, sub11, ld11, sub21, ld21);
+
+			status = cublasGetError();
+			if (status != CUBLAS_STATUS_SUCCESS)
+				STARPU_ABORT();
+
+			cudaThreadSynchronize();
+
+			break;
+#endif
+		default:
+			STARPU_ABORT();
+			break;
+	}
+}
+
+void STARPU_LU(cpu_u21)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u21)(descr, 0, _args);
+}
+
+#ifdef STARPU_USE_CUDA
+void STARPU_LU(cublas_u21)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u21)(descr, 1, _args);
+}
+#endif 
+
+static struct starpu_perfmodel_t STARPU_LU(model_21) = {
+	.type = STARPU_HISTORY_BASED,
+#ifdef STARPU_ATLAS
+	.symbol = STARPU_LU_STR(lu_model_21_atlas)
+#elif defined(STARPU_GOTO)
+	.symbol = STARPU_LU_STR(lu_model_21_goto)
+#else
+	.symbol = STARPU_LU_STR(lu_model_21)
+#endif
+};
+
+starpu_codelet cl21 = {
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = STARPU_LU(cpu_u21),
+#ifdef STARPU_USE_CUDA
+	.cuda_func = STARPU_LU(cublas_u21),
+#endif
+	.nbuffers = 2,
+	.model = &STARPU_LU(model_21)
+};
+
+/*
+ *	U11
+ */
+
+static inline void STARPU_LU(common_u11)(void *descr[],
+				int s, __attribute__((unused)) void *_args)
+{
+	TYPE *sub11;
+
+	sub11 = (TYPE *)STARPU_MATRIX_GET_PTR(descr[0]); 
+
+	unsigned long nx = STARPU_MATRIX_GET_NX(descr[0]);
+	unsigned long ld = STARPU_MATRIX_GET_LD(descr[0]);
+
+	unsigned long z;
+
+	switch (s) {
+		case 0:
+			for (z = 0; z < nx; z++)
+			{
+				TYPE pivot;
+				pivot = sub11[z+z*ld];
+				STARPU_ASSERT(pivot != 0.0);
+		
+				CPU_SCAL(nx - z - 1, (1.0/pivot), &sub11[z+(z+1)*ld], ld);
+		
+				CPU_GER(nx - z - 1, nx - z - 1, -1.0,
+						&sub11[(z+1)+z*ld], 1,
+						&sub11[z+(z+1)*ld], ld,
+						&sub11[(z+1) + (z+1)*ld],ld);
+			}
+			break;
+#ifdef STARPU_USE_CUDA
+		case 1:
+			for (z = 0; z < nx; z++)
+			{
+				TYPE pivot;
+				cudaMemcpy(&pivot, &sub11[z+z*ld], sizeof(TYPE), cudaMemcpyDeviceToHost);
+				cudaStreamSynchronize(0);
+
+				STARPU_ASSERT(pivot != 0.0);
+				
+				CUBLAS_SCAL(nx - z - 1, 1.0/pivot, &sub11[z+(z+1)*ld], ld);
+				
+				CUBLAS_GER(nx - z - 1, nx - z - 1, -1.0,
+						&sub11[(z+1)+z*ld], 1,
+						&sub11[z+(z+1)*ld], ld,
+						&sub11[(z+1) + (z+1)*ld],ld);
+			}
+			
+			cudaThreadSynchronize();
+
+			break;
+#endif
+		default:
+			STARPU_ABORT();
+			break;
+	}
+}
+
+void STARPU_LU(cpu_u11)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u11)(descr, 0, _args);
+}
+
+#ifdef STARPU_USE_CUDA
+void STARPU_LU(cublas_u11)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u11)(descr, 1, _args);
+}
+#endif// STARPU_USE_CUDA
+
+static struct starpu_perfmodel_t STARPU_LU(model_11) = {
+	.type = STARPU_HISTORY_BASED,
+#ifdef STARPU_ATLAS
+	.symbol = STARPU_LU_STR(lu_model_11_atlas)
+#elif defined(STARPU_GOTO)
+	.symbol = STARPU_LU_STR(lu_model_11_goto)
+#else
+	.symbol = STARPU_LU_STR(lu_model_11)
+#endif
+};
+
+starpu_codelet cl11 = {
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = STARPU_LU(cpu_u11),
+#ifdef STARPU_USE_CUDA
+	.cuda_func = STARPU_LU(cublas_u11),
+#endif
+	.nbuffers = 1,
+	.model = &STARPU_LU(model_11)
+};
+
+/*
+ *	U11 with pivoting
+ */
+
+static inline void STARPU_LU(common_u11_pivot)(void *descr[],
+				int s, void *_args)
+{
+	TYPE *sub11;
+
+	sub11 = (TYPE *)STARPU_MATRIX_GET_PTR(descr[0]); 
+
+	unsigned long nx = STARPU_MATRIX_GET_NX(descr[0]);
+	unsigned long ld = STARPU_MATRIX_GET_LD(descr[0]);
+
+	unsigned long z;
+
+	struct piv_s *piv = _args;
+	unsigned *ipiv = piv->piv;
+	unsigned first = piv->first;
+
+	switch (s) {
+		case 0:
+			for (z = 0; z < nx; z++)
+			{
+				TYPE pivot;
+				pivot = sub11[z+z*ld];
+
+				if (fabs((double)(pivot)) < PIVOT_THRESHHOLD)
+				{
+
+					/* find the pivot */
+					int piv_ind = CPU_IAMAX(nx - z, &sub11[z*(ld+1)], ld);
+
+					ipiv[z + first] = piv_ind + z + first;
+
+					/* swap if needed */
+					if (piv_ind != 0)
+					{
+						CPU_SWAP(nx, &sub11[z*ld], 1, &sub11[(z+piv_ind)*ld], 1);
+					}
+
+					pivot = sub11[z+z*ld];
+				}
+			
+				STARPU_ASSERT(pivot != 0.0);
+
+				CPU_SCAL(nx - z - 1, (1.0/pivot), &sub11[z+(z+1)*ld], ld);
+		
+				CPU_GER(nx - z - 1, nx - z - 1, -1.0,
+						&sub11[(z+1)+z*ld], 1,
+						&sub11[z+(z+1)*ld], ld,
+						&sub11[(z+1) + (z+1)*ld],ld);
+			}
+
+			break;
+#ifdef STARPU_USE_CUDA
+		case 1:
+			for (z = 0; z < nx; z++)
+			{
+				TYPE pivot;
+				cudaMemcpy(&pivot, &sub11[z+z*ld], sizeof(TYPE), cudaMemcpyDeviceToHost);
+				cudaStreamSynchronize(0);
+
+				if (fabs((double)(pivot)) < PIVOT_THRESHHOLD)
+				{
+					/* find the pivot */
+					int piv_ind = CUBLAS_IAMAX(nx - z, &sub11[z*(ld+1)], ld) - 1;
+	
+					ipiv[z + first] = piv_ind + z + first;
+
+					/* swap if needed */
+					if (piv_ind != 0)
+					{
+						CUBLAS_SWAP(nx, &sub11[z*ld], 1, &sub11[(z+piv_ind)*ld], 1);
+					}
+
+					cudaMemcpy(&pivot, &sub11[z+z*ld], sizeof(TYPE), cudaMemcpyDeviceToHost);
+					cudaStreamSynchronize(0);
+				}
+
+				STARPU_ASSERT(pivot != 0.0);
+				
+				CUBLAS_SCAL(nx - z - 1, 1.0/pivot, &sub11[z+(z+1)*ld], ld);
+				
+				CUBLAS_GER(nx - z - 1, nx - z - 1, -1.0,
+						&sub11[(z+1)+z*ld], 1,
+						&sub11[z+(z+1)*ld], ld,
+						&sub11[(z+1) + (z+1)*ld],ld);
+				
+			}
+
+			cudaThreadSynchronize();
+
+			break;
+#endif
+		default:
+			STARPU_ABORT();
+			break;
+	}
+}
+
+void STARPU_LU(cpu_u11_pivot)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u11_pivot)(descr, 0, _args);
+}
+
+#ifdef STARPU_USE_CUDA
+void STARPU_LU(cublas_u11_pivot)(void *descr[], void *_args)
+{
+	STARPU_LU(common_u11_pivot)(descr, 1, _args);
+}
+#endif// STARPU_USE_CUDA
+
+static struct starpu_perfmodel_t STARPU_LU(model_11_pivot) = {
+	.type = STARPU_HISTORY_BASED,
+#ifdef STARPU_ATLAS
+	.symbol = STARPU_LU_STR(lu_model_11_pivot_atlas)
+#elif defined(STARPU_GOTO)
+	.symbol = STARPU_LU_STR(lu_model_11_pivot_goto)
+#else
+	.symbol = STARPU_LU_STR(lu_model_11_pivot)
+#endif
+};
+
+starpu_codelet cl11_pivot = {
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = STARPU_LU(cpu_u11_pivot),
+#ifdef STARPU_USE_CUDA
+	.cuda_func = STARPU_LU(cublas_u11_pivot),
+#endif
+	.nbuffers = 1,
+	.model = &STARPU_LU(model_11_pivot)
+};
+
+/*
+ *	Pivoting
+ */
+
+static inline void STARPU_LU(common_pivot)(void *descr[],
+				int s, void *_args)
+{
+	TYPE *matrix;
+
+	matrix = (TYPE *)STARPU_MATRIX_GET_PTR(descr[0]); 
+	unsigned long nx = STARPU_MATRIX_GET_NX(descr[0]);
+	unsigned long ld = STARPU_MATRIX_GET_LD(descr[0]);
+
+	unsigned row;
+
+	struct piv_s *piv = _args;
+	unsigned *ipiv = piv->piv;
+	unsigned first = piv->first;
+
+	switch (s) {
+		case 0:
+			for (row = 0; row < nx; row++)
+			{
+				unsigned rowpiv = ipiv[row+first] - first;
+				if (rowpiv != row)
+				{
+					CPU_SWAP(nx, &matrix[row*ld], 1, &matrix[rowpiv*ld], 1);
+				}
+			}
+			break;
+#ifdef STARPU_USE_CUDA
+		case 1:
+			for (row = 0; row < nx; row++)
+			{
+				unsigned rowpiv = ipiv[row+first] - first;
+				if (rowpiv != row)
+				{
+					CUBLAS_SWAP(nx, &matrix[row*ld], 1, &matrix[rowpiv*ld], 1);
+				}
+			}
+
+			cudaThreadSynchronize();
+
+			break;
+#endif
+		default:
+			STARPU_ABORT();
+			break;
+	}
+}
+
+void STARPU_LU(cpu_pivot)(void *descr[], void *_args)
+{
+	STARPU_LU(common_pivot)(descr, 0, _args);
+}
+
+#ifdef STARPU_USE_CUDA
+void STARPU_LU(cublas_pivot)(void *descr[], void *_args)
+{
+	STARPU_LU(common_pivot)(descr, 1, _args);
+}
+
+#endif// STARPU_USE_CUDA
+
+static struct starpu_perfmodel_t STARPU_LU(model_pivot) = {
+	.type = STARPU_HISTORY_BASED,
+#ifdef STARPU_ATLAS
+	.symbol = STARPU_LU_STR(lu_model_pivot_atlas)
+#elif defined(STARPU_GOTO)
+	.symbol = STARPU_LU_STR(lu_model_pivot_goto)
+#else
+	.symbol = STARPU_LU_STR(lu_model_pivot)
+#endif
+};
+
+starpu_codelet cl_pivot = {
+	.where = STARPU_CPU|STARPU_CUDA,
+	.cpu_func = STARPU_LU(cpu_pivot),
+#ifdef STARPU_USE_CUDA
+	.cuda_func = STARPU_LU(cublas_pivot),
+#endif
+	.nbuffers = 1,
+	.model = &STARPU_LU(model_pivot)
+};

examples/cholesky_and_lu/lu/xlu_kernels.h

@@ -0,0 +1,46 @@
+/* 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.
+ */
+
+#ifndef __XLU_KERNELS_H__
+#define __XLU_KERNELS_H__
+
+#include <starpu.h>
+
+void STARPU_LU(cpu_pivot)(void *descr[], void *_args);
+void STARPU_LU(cpu_u11_pivot)(void *descr[], void *_args);
+void STARPU_LU(cpu_u11)(void *descr[], void *_args);
+void STARPU_LU(cpu_u12)(void *descr[], void *_args);
+void STARPU_LU(cpu_u21)(void *descr[], void *_args);
+void STARPU_LU(cpu_u22)(void *descr[], void *_args);
+
+#ifdef STARPU_USE_CUDA
+void STARPU_LU(cublas_pivot)(void *descr[], void *_args);
+void STARPU_LU(cublas_u11_pivot)(void *descr[], void *_args);
+void STARPU_LU(cublas_u11)(void *descr[], void *_args);
+void STARPU_LU(cublas_u12)(void *descr[], void *_args);
+void STARPU_LU(cublas_u21)(void *descr[], void *_args);
+void STARPU_LU(cublas_u22)(void *descr[], void *_args);
+#endif
+
+extern starpu_codelet cl11;
+extern starpu_codelet cl11_pivot;
+extern starpu_codelet cl12;
+extern starpu_codelet cl21;
+extern starpu_codelet cl22;
+extern starpu_codelet cl_pivot;
+
+#endif // __XLU_KERNELS_H__

examples/cholesky_and_lu/lu/xlu_pivot.c

@@ -0,0 +1,450 @@
+/* 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.
+ */
+
+#include "xlu.h"
+#include "xlu_kernels.h"
+
+#define TAG11(k)	((starpu_tag_t)( (1ULL<<60) | (unsigned long long)(k)))
+#define TAG12(k,i)	((starpu_tag_t)(((2ULL<<60) | (((unsigned long long)(k))<<32)	\
+					| (unsigned long long)(i))))
+#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 PIVOT(k,i)	((starpu_tag_t)(((5ULL<<60) | (((unsigned long long)(k))<<32)	\
+					| (unsigned long long)(i))))
+
+static unsigned no_prio = 0;
+starpu_data_handle xlu_pivot_dataA;
+starpu_data_handle *xlu_pivot_dataAp;
+struct timeval xlu_pivot_start;
+struct timeval xlu_pivot_end;
+struct timeval xlu_pivot_no_stride_start;
+struct timeval xlu_pivot_no_stride_end;
+
+
+/*
+ *	Construct the DAG
+ */
+
+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;
+}
+
+static void create_task_pivot(starpu_data_handle *dataAp, unsigned nblocks,
+					struct piv_s *piv_description,
+					unsigned k, unsigned i,
+			      starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = create_task(PIVOT(k, i));
+
+	task->cl = &cl_pivot;
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, i);
+	task->buffers[0].mode = STARPU_RW;
+
+	task->cl_arg = &piv_description[k];
+
+	/* this is an important task */
+	if (!no_prio && (i == k+1))
+		task->priority = STARPU_MAX_PRIO;
+
+	/* enforce dependencies ... */
+	if (k == 0) {
+		starpu_tag_declare_deps(PIVOT(k, i), 1, TAG11(k));
+	}
+	else 
+	{
+		if (i > k) {
+			starpu_tag_declare_deps(PIVOT(k, i), 2, TAG11(k), TAG22(k-1, i, k));
+		}
+		else {
+			starpu_tag_t *tags = malloc((nblocks - k)*sizeof(starpu_tag_t));
+			
+			tags[0] = TAG11(k);
+			unsigned ind, ind2;
+			for (ind = k + 1, ind2 = 0; ind < nblocks; ind++, ind2++)
+			{
+				tags[1 + ind2] = TAG22(k-1, ind, k);
+			}
+
+			/* perhaps we could do better ... :/  */
+			starpu_tag_declare_deps_array(PIVOT(k, i), (nblocks-k), tags);
+		}
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static struct starpu_task *create_task_11_pivot(starpu_data_handle *dataAp, unsigned nblocks,
+					unsigned k, struct piv_s *piv_description,
+					starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned))
+{
+	struct starpu_task *task = create_task(TAG11(k));
+
+	task->cl = &cl11_pivot;
+
+	task->cl_arg = &piv_description[k];
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, k);
+	task->buffers[0].mode = STARPU_RW;
+
+	/* this is an important task */
+	if (!no_prio)
+		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 void create_task_12(starpu_data_handle *dataAp, unsigned nblocks, unsigned k, unsigned j,
+			   starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx)
+{
+//	printf("task 12 k,i = %d,%d TAG = %llx\n", k,i, TAG12(k,i));
+
+	struct starpu_task *task = create_task(TAG12(k, j));
+	
+	task->cl = &cl12;
+
+	task->cl_arg = (void *)(task->tag_id);
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, k);
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = get_block(dataAp, nblocks, j, k);
+	task->buffers[1].mode = STARPU_RW;
+
+	if (!no_prio && (j == k+1)) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	/* enforce dependencies ... */
+#if 0
+	starpu_tag_declare_deps(TAG12(k, i), 1, PIVOT(k, i));
+#endif
+	if (k > 0) {
+		starpu_tag_declare_deps(TAG12(k, j), 2, TAG11(k), TAG22(k-1, k, j));
+	}
+	else {
+		starpu_tag_declare_deps(TAG12(k, j), 1, TAG11(k));
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_21(starpu_data_handle *dataAp, unsigned nblocks, unsigned k, unsigned i,
+				starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx)
+{
+	struct starpu_task *task = create_task(TAG21(k, i));
+
+	task->cl = &cl21;
+	
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, k); 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = get_block(dataAp, nblocks, k, i); 
+	task->buffers[1].mode = STARPU_RW;
+
+	if (!no_prio && (i == k+1)) {
+		task->priority = STARPU_MAX_PRIO;
+	}
+
+	task->cl_arg = (void *)(task->tag_id);
+
+	/* enforce dependencies ... */
+	starpu_tag_declare_deps(TAG21(k, i), 1, PIVOT(k, i));
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static void create_task_22(starpu_data_handle *dataAp, unsigned nblocks, unsigned k, unsigned i, unsigned j,
+				starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), 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;
+
+	task->cl_arg = (void *)(task->tag_id);
+
+	/* which sub-data is manipulated ? */
+	task->buffers[0].handle = get_block(dataAp, nblocks, k, i); /* produced by TAG21(k, i) */
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = get_block(dataAp, nblocks, j, k); /* produced by TAG12(k, j) */ 
+	task->buffers[1].mode = STARPU_R;
+	task->buffers[2].handle = get_block(dataAp, nblocks, j, i);  /* produced by TAG22(k-1, i, j) */
+	task->buffers[2].mode = STARPU_RW;
+
+	if (!no_prio &&  (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), TAG12(k, j), TAG21(k, i));
+	}
+	else {
+		starpu_tag_declare_deps(TAG22(k, i, j), 2, TAG12(k, j), TAG21(k, i));
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+/*
+ *	code to bootstrap the factorization 
+ */
+
+static double dw_codelet_facto_pivot(starpu_data_handle *dataAp,
+				     struct piv_s *piv_description,
+				     unsigned nblocks,
+				     starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), struct starpu_sched_ctx *sched_ctx,
+				     struct timeval *start)
+{
+	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_pivot(dataAp, nblocks, k, piv_description, get_block);
+
+		/* we defer the launch of the first task */
+		if (k == 0) {
+			entry_task = task;
+		}
+		else {
+		  starpu_task_submit_to_ctx(task, sched_ctx);
+		}
+
+		for (i = 0; i < nblocks; i++)
+		{
+			if (i != k)
+			  create_task_pivot(dataAp, nblocks, piv_description, k, i, get_block, sched_ctx);
+		}
+	
+		for (i = k+1; i<nblocks; i++)
+		{
+		  create_task_12(dataAp, nblocks, k, i, get_block, sched_ctx);
+		  create_task_21(dataAp, nblocks, k, i, get_block, sched_ctx);
+		}
+
+		for (i = k+1; i<nblocks; i++)
+		{
+			for (j = k+1; j<nblocks; j++)
+			{
+			  create_task_22(dataAp, nblocks, k, i, j, get_block, sched_ctx);
+			}
+		}
+	}
+
+	/* schedule the codelet */
+	gettimeofday(start, NULL);
+	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);
+	}
+
+	return 0;
+}
+
+starpu_data_handle get_block_with_striding(starpu_data_handle *dataAp,
+			unsigned nblocks __attribute__((unused)), unsigned j, unsigned i)
+{
+	/* we use filters */
+	return starpu_data_get_sub_data(*dataAp, 2, j, i);
+}
+
+
+void STARPU_LU(lu_decomposition_pivot)(TYPE *matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+
+	/* monitor and partition the A matrix into blocks :
+	 * one block is now determined by 2 unsigned (i,j) */
+	starpu_matrix_data_register(&xlu_pivot_dataA, 0, (uintptr_t)matA, ld, size, size, sizeof(TYPE));
+
+	/* We already enforce deps by hand */
+	starpu_data_set_sequential_consistency_flag(xlu_pivot_dataA, 0);
+
+	struct starpu_data_filter f;
+		f.filter_func = starpu_vertical_block_filter_func;
+		f.nchildren = nblocks;
+
+	struct starpu_data_filter f2;
+		f2.filter_func = starpu_block_filter_func;
+		f2.nchildren = nblocks;
+
+	starpu_data_map_filters(xlu_pivot_dataA, 2, &f, &f2);
+
+	unsigned i;
+	for (i = 0; i < size; i++)
+		ipiv[i] = i;
+
+	struct piv_s *piv_description = malloc(nblocks*sizeof(struct piv_s));
+	unsigned block;
+	for (block = 0; block < nblocks; block++)
+	{
+		piv_description[block].piv = ipiv;
+		piv_description[block].first = block * (size / nblocks);
+		piv_description[block].last = (block + 1) * (size / nblocks);
+	}
+
+#if 0
+	unsigned j;
+	for (j = 0; j < nblocks; j++)
+	for (i = 0; i < nblocks; i++)
+	{
+		printf("BLOCK %d %d	%p\n", i, j, &matA[i*(size/nblocks) + j * (size/nblocks)*ld]);
+	}
+#endif
+
+	dw_codelet_facto_pivot(&xlu_pivot_dataA, piv_description, nblocks, get_block_with_striding, sched_ctx, &xlu_pivot_start);
+}
+
+
+starpu_data_handle get_block_with_no_striding(starpu_data_handle *dataAp, unsigned nblocks, unsigned j, unsigned i)
+{
+	/* dataAp is an array of data handle */
+	return dataAp[i+j*nblocks];
+}
+
+void STARPU_LU(lu_decomposition_pivot_no_stride)(TYPE **matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	xlu_pivot_dataAp = malloc(nblocks*nblocks*sizeof(starpu_data_handle));
+
+	/* monitor and partition the A matrix into blocks :
+	 * one block is now determined by 2 unsigned (i,j) */
+	unsigned bi, bj;
+	for (bj = 0; bj < nblocks; bj++)
+	for (bi = 0; bi < nblocks; bi++)
+	{
+		starpu_matrix_data_register(&xlu_pivot_dataAp[bi+nblocks*bj], 0,
+			(uintptr_t)matA[bi+nblocks*bj], size/nblocks,
+			size/nblocks, size/nblocks, sizeof(TYPE));
+
+		/* We already enforce deps by hand */
+		starpu_data_set_sequential_consistency_flag(xlu_pivot_dataAp[bi+nblocks*bj], 0);
+	}
+
+	unsigned i;
+	for (i = 0; i < size; i++)
+		ipiv[i] = i;
+
+	struct piv_s *piv_description = malloc(nblocks*sizeof(struct piv_s));
+	unsigned block;
+	for (block = 0; block < nblocks; block++)
+	{
+		piv_description[block].piv = ipiv;
+		piv_description[block].first = block * (size / nblocks);
+		piv_description[block].last = (block + 1) * (size / nblocks);
+	}
+
+	dw_codelet_facto_pivot(xlu_pivot_dataAp, piv_description, nblocks, get_block_with_no_striding, sched_ctx, &xlu_pivot_no_stride_start);
+}
+
+void finish_lu_decomposition_pivot(unsigned nblocks)
+{
+	/* we wait the last task (TAG11(nblocks - 1)) and all the pivot tasks */
+	starpu_tag_t *tags = malloc(nblocks*nblocks*sizeof(starpu_tag_t));
+	unsigned ndeps = 0;
+
+	tags[ndeps++] = TAG11(nblocks - 1);
+
+	unsigned i, j;
+	for (j = 0; j < nblocks; j++)
+	{
+		for (i = 0; i < j; i++)
+		{
+			tags[ndeps++] = PIVOT(j, i);
+		}
+	}
+
+	/* stall the application until the end of computations */
+	starpu_tag_wait_array(ndeps, tags);
+//	starpu_task_wait_for_all();
+
+	gettimeofday(&xlu_pivot_end, NULL);
+
+	double timing = (double)((&xlu_pivot_end.tv_sec - &xlu_pivot_start.tv_sec)*1000000 + (&xlu_pivot_end.tv_usec - &xlu_pivot_start.tv_usec));
+
+	fprintf(stderr, "Computation took (in ms)\n");
+	fprintf(stderr, "%2.2f\n", timing/1000);
+
+	unsigned n = starpu_matrix_get_nx(xlu_pivot_dataA);
+	double flop = (2.0f*n*n*n)/3.0f;
+	fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
+
+	/* gather all the data */
+	starpu_data_unpartition(xlu_pivot_dataA, 0);
+}
+
+void finish_lu_decomposition_pivot_no_stride(unsigned nblocks)
+{
+	/* we wait the last task (TAG11(nblocks - 1)) and all the pivot tasks */
+	starpu_tag_t *tags = malloc(nblocks*nblocks*sizeof(starpu_tag_t));
+	unsigned ndeps = 0;
+
+	tags[ndeps++] = TAG11(nblocks - 1);
+
+	unsigned i, j;
+	for (j = 0; j < nblocks; j++)
+	{
+		for (i = 0; i < j; i++)
+		{
+			tags[ndeps++] = PIVOT(j, i);
+		}
+	}
+
+	/* stall the application until the end of computations */
+	starpu_tag_wait_array(ndeps, tags);
+//	starpu_task_wait_for_all();
+
+	gettimeofday(&xlu_pivot_no_stride_end, NULL);
+
+	double timing = (double)((&xlu_pivot_no_stride_end.tv_sec - &xlu_pivot_no_stride_start.tv_sec)*1000000 + (&xlu_pivot_no_stride_end.tv_usec - &xlu_pivot_no_stride_start.tv_usec));
+
+	fprintf(stderr, "Computation took (in ms)\n");
+	fprintf(stderr, "%2.2f\n", timing/1000);
+
+	unsigned n = starpu_matrix_get_nx(xlu_pivot_dataAp[0])*nblocks;
+	double flop = (2.0f*n*n*n)/3.0f;
+	fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
+
+	unsigned bi, bj;
+	for (bj = 0; bj < nblocks; bj++)
+	for (bi = 0; bi < nblocks; bi++)
+	{
+		starpu_data_unregister(xlu_pivot_dataAp[bi+nblocks*bj]);
+	}
+}