Repaired conflict

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,362 @@
+/* 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 = 10240;//4096;
+static unsigned lu_nblocks = 10;
+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;
+
+TYPE *A;
+TYPE *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()
+{
+	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+lu_nblocks*bj]);
+	}
+}
+
+
+
+void copy_matrix_into_blocks()
+{
+	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()
+{
+	/* 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()
+{
+	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 lu_check_result()
+{
+	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);
+}
+
+double run_lu(struct starpu_sched_ctx *sched_ctx, int argc, char **argv)
+{
+  printf("enter lu\n");
+	lu_parse_args(argc, argv);
+
+	//	starpu_init(NULL);
+
+	//	starpu_helper_cublas_init();
+
+	init_matrix();
+
+	unsigned *ipiv;
+	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);
+
+	double gflops = -1;
+	/* 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();
+
+			gflops = 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 
+		{
+			struct timeval start;
+			struct timeval end;
+
+			gettimeofday(&start, NULL);
+
+			gflops = STARPU_LU(lu_decomposition_pivot)(A, ipiv, lu_size, lu_size, lu_nblocks, sched_ctx);
+	
+			gettimeofday(&end, NULL);
+
+			double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+			
+			unsigned n = lu_size;
+			double flop = (2.0f*n*n*n)/3.0f;
+			gflops = flop/timing/1000.0f;
+		}
+	}
+	else
+	{
+	  gflops = STARPU_LU(lu_decomposition)(A, lu_size, lu_size, lu_nblocks, sched_ctx);
+			printf("no pivot \n");
+	}
+
+	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);
+
+		lu_check_result();
+	}
+
+	//	starpu_helper_cublas_shutdown();
+
+	//	starpu_shutdown();
+
+	return gflops;
+}

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,257 @@
+/* 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;
+
+
+
+
+/*
+ *	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 double dw_codelet_facto_v3(starpu_data_handle dataA, unsigned lu_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;
+
+	for (k = 0; k < lu_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<lu_nblocks; i++)
+		{
+		  create_task_12(dataA, k, i, sched_ctx);
+		  create_task_21(dataA, k, i, sched_ctx);
+		}
+
+		for (i = k+1; i<lu_nblocks; i++)
+		{
+			for (j = k+1; j<lu_nblocks; j++)
+			{
+			  create_task_22(dataA, k, i, j, 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);
+	}
+
+
+
+	/* stall the application until the end of computations */
+	starpu_tag_wait(TAG11(lu_nblocks-1));
+	printf("lu finish waiting for %d nblocks\n", lu_nblocks-1);
+
+	gettimeofday(&end, NULL);
+
+	double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+
+	unsigned n = starpu_matrix_get_nx(dataA);
+	double flop = (2.0f*n*n*n)/3.0f;
+	return (flop/timing/1000.0f);
+}
+
+double STARPU_LU(lu_decomposition)(TYPE *matA, unsigned size, unsigned ld, unsigned lu_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));
+
+	/* We already enforce deps by hand */
+	starpu_data_set_sequential_consistency_flag(dataA, 0);
+
+	struct starpu_data_filter f;
+		f.filter_func = starpu_vertical_block_filter_func;
+		f.nchildren = lu_nblocks;
+		f.get_nchildren = NULL;
+		f.get_child_ops = NULL;
+
+	struct starpu_data_filter f2;
+		f2.filter_func = starpu_block_filter_func;
+		f2.nchildren = lu_nblocks;
+		f2.get_nchildren = NULL;
+		f2.get_child_ops = NULL;
+
+	starpu_data_map_filters(dataA, 2, &f, &f2);
+
+	double gflops = dw_codelet_facto_v3(dataA, lu_nblocks, sched_ctx);
+
+	/* gather all the data */
+	starpu_data_unpartition(dataA, 0);
+	return gflops;
+}

examples/cholesky_and_lu/lu/xlu.h

@@ -0,0 +1,117 @@
+/* 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 */
+};
+
+double STARPU_LU(lu_decomposition)(TYPE *matA, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx);
+double STARPU_LU(lu_decomposition_pivot_no_stride)(TYPE **matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx);
+double STARPU_LU(lu_decomposition_pivot)(TYPE *matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx);
+
+#endif // __XLU_H__

examples/cholesky_and_lu/lu/xlu_implicit.c

@@ -0,0 +1,166 @@
+/* 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 double dw_codelet_facto_v3(starpu_data_handle dataA, unsigned lu_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 < lu_nblocks; k++)
+	{
+	  create_task_11(dataA, k, sched_ctx);
+		
+		for (i = k+1; i<lu_nblocks; i++)
+		{
+		  create_task_12(dataA, k, i, sched_ctx);
+		  create_task_21(dataA, k, i, sched_ctx);
+		}
+
+		for (i = k+1; i<lu_nblocks; i++)
+		for (j = k+1; j<lu_nblocks; j++)
+		  create_task_22(dataA, k, i, j, sched_ctx);
+		  
+	}
+
+	/* stall the application until the end of computations */
+	//starpu_task_wait_for_all();
+	starpu_wait_for_all_tasks_of_sched_ctx(sched_ctx);
+
+	gettimeofday(&end, NULL);
+
+	double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
+
+	unsigned n = starpu_matrix_get_nx(dataA);
+	double flop = (2.0f*n*n*n)/3.0f;
+	return (flop/timing/1000.0f);
+}
+
+double STARPU_LU(lu_decomposition)(TYPE *matA, unsigned size, unsigned ld, unsigned lu_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 = lu_nblocks;
+		f.get_nchildren = NULL;
+		f.get_child_ops = NULL;
+
+	struct starpu_data_filter f2;
+		f2.filter_func = starpu_block_filter_func;
+		f2.nchildren = lu_nblocks;
+		f2.get_nchildren = NULL;
+		f2.get_child_ops = NULL;
+
+	starpu_data_map_filters(dataA, 2, &f, &f2);
+
+	double gflops = dw_codelet_facto_v3(dataA, lu_nblocks, sched_ctx);
+
+	/* gather all the data */
+	starpu_data_unpartition(dataA, 0);
+	starpu_data_unregister(dataA);
+	return gflops;
+}

examples/cholesky_and_lu/lu/xlu_implicit_pivot.c

@@ -0,0 +1,281 @@
+/* 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();
+	starpu_wait_for_all_tasks_of_sched_ctx(sched_ctx);
+
+	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);
+}
+
+
+double 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);
+
+	unsigned n = starpu_matrix_get_nx(dataA);
+	double flop = (2.0f*n*n*n)/3.0f;
+	double gflops = flop/timing/1000.0f;
+
+	/* gather all the data */
+	starpu_data_unpartition(dataA, 0);
+	starpu_data_unregister(dataA);
+	return gflops;
+}
+
+
+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];
+}
+
+double STARPU_LU(lu_decomposition_pivot_no_stride)(TYPE **matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	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);
+
+	unsigned n = starpu_matrix_get_nx(dataAp[0])*nblocks;
+	double flop = (2.0f*n*n*n)/3.0f;
+	double gflops = flop/timing/1000.0f;
+
+	for (bj = 0; bj < nblocks; bj++)
+	for (bi = 0; bi < nblocks; bi++)
+	{
+		starpu_data_unregister(dataAp[bi+nblocks*bj]);
+	}
+	return gflops;
+}

examples/cholesky_and_lu/lu/xlu_kernels.c

@@ -0,0 +1,605 @@
+/* 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:
+			status = cublasGetError();
+			if (STARPU_UNLIKELY(status != CUBLAS_STATUS_SUCCESS))
+				STARPU_ABORT();
+
+			if (STARPU_UNLIKELY((cures = cudaThreadSynchronize()) != cudaSuccess))
+				STARPU_CUDA_REPORT_ERROR(cures);
+
+
+			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;
+
+#ifdef STARPU_USE_CUDA
+	cublasStatus status;
+	cudaError_t cures;
+#endif
+
+	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);
+			}
+
+			status = cublasGetError();
+			if (STARPU_UNLIKELY(status != CUBLAS_STATUS_SUCCESS))
+				STARPU_ABORT();
+
+			if (STARPU_UNLIKELY((cures = cudaThreadSynchronize()) != cudaSuccess))
+				STARPU_CUDA_REPORT_ERROR(cures);
+
+			
+			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,421 @@
+/* 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;
+
+/*
+ *	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 lu_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, lu_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((lu_nblocks - k)*sizeof(starpu_tag_t));
+			
+			tags[0] = TAG11(k);
+			unsigned ind, ind2;
+			for (ind = k + 1, ind2 = 0; ind < lu_nblocks; ind++, ind2++)
+			{
+				tags[1 + ind2] = TAG22(k-1, ind, k);
+			}
+
+			/* perhaps we could do better ... :/  */
+			starpu_tag_declare_deps_array(PIVOT(k, i), (lu_nblocks-k), tags);
+		}
+	}
+
+	starpu_task_submit_to_ctx(task, sched_ctx);
+}
+
+static struct starpu_task *create_task_11_pivot(starpu_data_handle *dataAp, unsigned lu_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, lu_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 lu_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, lu_nblocks, k, k);
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = get_block(dataAp, lu_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 lu_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, lu_nblocks, k, k); 
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = get_block(dataAp, lu_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 lu_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, lu_nblocks, k, i); /* produced by TAG21(k, i) */
+	task->buffers[0].mode = STARPU_R;
+	task->buffers[1].handle = get_block(dataAp, lu_nblocks, j, k); /* produced by TAG12(k, j) */ 
+	task->buffers[1].mode = STARPU_R;
+	task->buffers[2].handle = get_block(dataAp, lu_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 lu_nblocks,
+				     starpu_data_handle (* get_block)(starpu_data_handle *, unsigned, unsigned, unsigned), 
+				     struct starpu_sched_ctx *sched_ctx)
+{
+	struct timeval start;
+	struct timeval end;
+
+	struct starpu_task *entry_task = NULL;
+
+	/* create all the DAG nodes */
+	unsigned i,j,k;
+
+	for (k = 0; k < lu_nblocks; k++)
+	{
+		struct starpu_task *task = create_task_11_pivot(dataAp, lu_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 < lu_nblocks; i++)
+		{
+			if (i != k)
+			  create_task_pivot(dataAp, lu_nblocks, piv_description, k, i, get_block, sched_ctx);
+		}
+	
+		for (i = k+1; i<lu_nblocks; i++)
+		{
+		  create_task_12(dataAp, lu_nblocks, k, i, get_block, sched_ctx);
+		  create_task_21(dataAp, lu_nblocks, k, i, get_block, sched_ctx);
+		}
+
+		for (i = k+1; i<lu_nblocks; i++)
+		{
+			for (j = k+1; j<lu_nblocks; j++)
+			{
+			  create_task_22(dataAp, lu_nblocks, k, i, j, get_block, sched_ctx);
+			}
+		}
+	}
+
+	/* we wait the last task (TAG11(lu_nblocks - 1)) and all the pivot tasks */
+	starpu_tag_t *tags = malloc(lu_nblocks*lu_nblocks*sizeof(starpu_tag_t));
+	unsigned ndeps = 0;
+
+	tags[ndeps++] = TAG11(lu_nblocks - 1);
+
+	for (j = 0; j < lu_nblocks; j++)
+	{
+		for (i = 0; i < j; i++)
+		{
+			tags[ndeps++] = PIVOT(j, i);
+		}
+	}
+
+	/* 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);
+	}
+
+	/* stall the application until the end of computations */
+	starpu_tag_wait_array(ndeps, tags);
+	printf("lu pivot  finish waiting for %d blocks \n", lu_nblocks);
+//	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 lu_nblocks __attribute__((unused)), unsigned j, unsigned i)
+{
+	/* we use filters */
+	return starpu_data_get_sub_data(*dataAp, 2, j, i);
+}
+
+
+double STARPU_LU(lu_decomposition_pivot)(TYPE *matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned lu_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));
+
+	/* We already enforce deps by hand */
+	starpu_data_set_sequential_consistency_flag(dataA, 0);
+
+	struct starpu_data_filter f;
+		f.filter_func = starpu_vertical_block_filter_func;
+		f.nchildren = lu_nblocks;
+
+	struct starpu_data_filter f2;
+		f2.filter_func = starpu_block_filter_func;
+		f2.nchildren = lu_nblocks;
+
+	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(lu_nblocks*sizeof(struct piv_s));
+	unsigned block;
+	for (block = 0; block < lu_nblocks; block++)
+	{
+		piv_description[block].piv = ipiv;
+		piv_description[block].first = block * (size / lu_nblocks);
+		piv_description[block].last = (block + 1) * (size / lu_nblocks);
+	}
+
+#if 0
+	unsigned j;
+	for (j = 0; j < lu_nblocks; j++)
+	for (i = 0; i < lu_nblocks; i++)
+	{
+		printf("BLOCK %d %d	%p\n", i, j, &matA[i*(size/lu_nblocks) + j * (size/lu_nblocks)*ld]);
+	}
+#endif
+
+	double timing;
+	timing = dw_codelet_facto_pivot(&dataA, piv_description, lu_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;
+	double gflops = flop/timing/1000.0f;
+
+	/* gather all the data */
+	starpu_data_unpartition(dataA, 0);
+	return gflops;
+}
+
+
+starpu_data_handle get_block_with_no_striding(starpu_data_handle *dataAp, unsigned lu_nblocks, unsigned j, unsigned i)
+{
+	/* dataAp is an array of data handle */
+	return dataAp[i+j*lu_nblocks];
+}
+
+double STARPU_LU(lu_decomposition_pivot_no_stride)(TYPE **matA, unsigned *ipiv, unsigned size, unsigned ld, unsigned lu_nblocks, struct starpu_sched_ctx *sched_ctx)
+{
+	starpu_data_handle *dataAp = malloc(lu_nblocks*lu_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 < lu_nblocks; bj++)
+	for (bi = 0; bi < lu_nblocks; bi++)
+	{
+		starpu_matrix_data_register(&dataAp[bi+lu_nblocks*bj], 0,
+			(uintptr_t)matA[bi+lu_nblocks*bj], size/lu_nblocks,
+			size/lu_nblocks, size/lu_nblocks, sizeof(TYPE));
+
+		/* We already enforce deps by hand */
+		starpu_data_set_sequential_consistency_flag(dataAp[bi+lu_nblocks*bj], 0);
+	}
+
+	unsigned i;
+	for (i = 0; i < size; i++)
+		ipiv[i] = i;
+
+	struct piv_s *piv_description = malloc(lu_nblocks*sizeof(struct piv_s));
+	unsigned block;
+	for (block = 0; block < lu_nblocks; block++)
+	{
+		piv_description[block].piv = ipiv;
+		piv_description[block].first = block * (size / lu_nblocks);
+		piv_description[block].last = (block + 1) * (size / lu_nblocks);
+	}
+
+	double timing;
+	timing = dw_codelet_facto_pivot(dataAp, piv_description, lu_nblocks, get_block_with_no_striding, sched_ctx);
+
+	unsigned n = starpu_matrix_get_nx(dataAp[0])*lu_nblocks;
+	double flop = (2.0f*n*n*n)/3.0f;
+	double gflops = flop/timing/1000.0f;
+
+	for (bj = 0; bj < lu_nblocks; bj++)
+	for (bi = 0; bi < lu_nblocks; bi++)
+	{
+		starpu_data_unregister(dataAp[bi+lu_nblocks*bj]);
+	}
+	return gflops;
+}