Author: bart
Date: 2008-03-09 13:41:26 +0000 (Sun, 09 Mar 2008)
New Revision: 7614
Log:
Added OpenMP test program.
Added:
trunk/exp-drd/tests/matinv_openmp.c
Modified:
trunk/exp-drd/tests/Makefile.am
Modified: trunk/exp-drd/tests/Makefile.am
===================================================================
--- trunk/exp-drd/tests/Makefile.am 2008-03-09 13:39:58 UTC (rev 7613)
+++ trunk/exp-drd/tests/Makefile.am 2008-03-09 13:41:26 UTC (rev 7614)
@@ -195,6 +195,7 @@
hg05_race2 \
hg06_readshared \
matinv \
+ matinv_openmp \
pth_barrier \
pth_broadcast \
pth_cond_race \
@@ -256,6 +257,10 @@
matinv_SOURCES = matinv.c
matinv_LDADD = -lpthread -lm
+matinv_openmp_SOURCES = matinv_openmp.c
+matinv_openmp_CFLAGS = -fopenmp
+matinv_openmp_LDADD = -lpthread -lm
+
pth_barrier_SOURCES = pth_barrier.c
pth_barrier_LDADD = -lpthread
Added: trunk/exp-drd/tests/matinv_openmp.c
===================================================================
--- trunk/exp-drd/tests/matinv_openmp.c (rev 0)
+++ trunk/exp-drd/tests/matinv_openmp.c 2008-03-09 13:41:26 UTC (rev 7614)
@@ -0,0 +1,293 @@
+/** Compute the matrix inverse via Gauss-Jordan elimination.
+ * This program uses OpenMP separate computation steps but no
+ * mutexes. It is an example of a race-free program on which no data races
+ * are reported by the happens-before algorithm (drd), but a lot of data races
+ * (all false positives) are reported by the Eraser-algorithm (helgrind).
+ */
+
+
+#define _GNU_SOURCE
+
+/***********************/
+/* Include directives. */
+/***********************/
+
+#include <assert.h>
+#include <math.h>
+#include <pthread.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+
+/*********************/
+/* Type definitions. */
+/*********************/
+
+typedef double elem_t;
+
+
+/********************/
+/* Local variables. */
+/********************/
+
+static int s_nthread;
+
+
+/*************************/
+/* Function definitions. */
+/*************************/
+
+/** Allocate memory for a matrix with the specified number of rows and
+ * columns.
+ */
+static elem_t* new_matrix(const int rows, const int cols)
+{
+ assert(rows > 0);
+ assert(cols > 0);
+ return malloc(rows * cols * sizeof(elem_t));
+}
+
+/** Free the memory that was allocated for a matrix. */
+static void delete_matrix(elem_t* const a)
+{
+ free(a);
+}
+
+/** Fill in some numbers in a matrix.
+ * @note It is important not to call srand() in this program, such that
+ * the results of a run are reproducible.
+ */
+static void init_matrix(elem_t* const a, const int rows, const int cols)
+{
+ int i, j;
+ for (i = 0; i < rows; i++)
+ {
+ for (j = 0; j < rows; j++)
+ {
+ a[i * cols + j] = rand() * 1.0 / RAND_MAX;
+ }
+ }
+}
+
+/** Print all elements of a matrix. */
+void print_matrix(const char* const label,
+ const elem_t* const a, const int rows, const int cols)
+{
+ int i, j;
+ printf("%s:\n", label);
+ for (i = 0; i < rows; i++)
+ {
+ for (j = 0; j < cols; j++)
+ {
+ printf("%g ", a[i * cols + j]);
+ }
+ printf("\n");
+ }
+}
+
+/** Copy a subset of the elements of a matrix into another matrix. */
+static void copy_matrix(const elem_t* const from,
+ const int from_rows,
+ const int from_cols,
+ const int from_row_first,
+ const int from_row_last,
+ const int from_col_first,
+ const int from_col_last,
+ elem_t* const to,
+ const int to_rows,
+ const int to_cols,
+ const int to_row_first,
+ const int to_row_last,
+ const int to_col_first,
+ const int to_col_last)
+{
+ int i, j;
+
+ assert(from_row_last - from_row_first == to_row_last - to_row_first);
+ assert(from_col_last - from_col_first == to_col_last - to_col_first);
+
+ for (i = from_row_first; i < from_row_last; i++)
+ {
+ assert(i < from_rows);
+ assert(i - from_row_first + to_row_first < to_rows);
+ for (j = from_col_first; j < from_col_last; j++)
+ {
+ assert(j < from_cols);
+ assert(j - from_col_first + to_col_first < to_cols);
+ to[(i - from_row_first + to_col_first) * to_cols
+ + (j - from_col_first + to_col_first)]
+ = from[i * from_cols + j];
+ }
+ }
+}
+
+/** Compute the matrix product of a1 and a2. */
+static elem_t* multiply_matrices(const elem_t* const a1,
+ const int rows1,
+ const int cols1,
+ const elem_t* const a2,
+ const int rows2,
+ const int cols2)
+{
+ int i, j, k;
+ elem_t* prod;
+
+ assert(cols1 == rows2);
+
+ prod = new_matrix(rows1, cols2);
+ for (i = 0; i < rows1; i++)
+ {
+ for (j = 0; j < cols2; j++)
+ {
+ prod[i * cols2 + j] = 0;
+ for (k = 0; k < cols1; k++)
+ {
+ prod[i * cols2 + j] += a1[i * cols1 + k] * a2[k * cols2 + j];
+ }
+ }
+ }
+ return prod;
+}
+
+/** Apply the Gauss-Jordan elimination algorithm on the matrix p->a starting
+ * at row r0 and up to but not including row r1. It is assumed that as many
+ * threads execute this function concurrently as the count barrier p->b was
+ * initialized with. If the matrix p->a is nonsingular, and if matrix p->a
+ * has at least as many columns as rows, the result of this algorithm is that
+ * submatrix p->a[0..p->rows-1,0..p->rows-1] is the identity matrix.
+ * @see http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
+ */
+static void gj(elem_t* const a, const int rows, const int cols)
+{
+ int i, j, k;
+
+ for (i = 0; i < rows; i++)
+ {
+ {
+ // Pivoting.
+ j = i;
+ for (k = i + 1; k < rows; k++)
+ {
+ if (a[k * cols + i] > a[j * cols + i])
+ {
+ j = k;
+ }
+ }
+ if (j != i)
+ {
+ for (k = 0; k < cols; k++)
+ {
+ const elem_t t = a[i * cols + k];
+ a[i * cols + k] = a[j * cols + k];
+ a[j * cols + k] = t;
+ }
+ }
+ // Normalize row i.
+ if (a[i * cols + i] != 0)
+ {
+ for (k = cols - 1; k >= 0; k--)
+ {
+ a[i * cols + k] /= a[i * cols + i];
+ }
+ }
+ }
+
+ // Reduce all rows j != i.
+#pragma omp parallel for
+ for (j = 0; j < rows; j++)
+ {
+ if (i != j)
+ {
+ const elem_t factor = a[j * cols + i];
+ for (k = 0; k < cols; k++)
+ {
+ a[j * cols + k] -= a[i * cols + k] * factor;
+ }
+ }
+ }
+ }
+}
+
+/** Matrix inversion via the Gauss-Jordan algorithm. */
+static elem_t* invert_matrix(const elem_t* const a, const int n)
+{
+ int i, j;
+ elem_t* const inv = new_matrix(n, n);
+ elem_t* const tmp = new_matrix(n, 2*n);
+ copy_matrix(a, n, n, 0, n, 0, n, tmp, n, 2 * n, 0, n, 0, n);
+ for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++)
+ tmp[i * 2 * n + n + j] = (i == j);
+ gj(tmp, n, 2*n);
+ copy_matrix(tmp, n, 2*n, 0, n, n, 2*n, inv, n, n, 0, n, 0, n);
+ delete_matrix(tmp);
+ return inv;
+}
+
+/** Compute the average square error between the identity matrix and the
+ * product of matrix a with its inverse matrix.
+ */
+static double identity_error(const elem_t* const a, const int n)
+{
+ int i, j;
+ elem_t e = 0;
+ for (i = 0; i < n; i++)
+ {
+ for (j = 0; j < n; j++)
+ {
+ const elem_t d = a[i * n + j] - (i == j);
+ e += d * d;
+ }
+ }
+ return sqrt(e / (n * n));
+}
+
+/** Compute epsilon for the numeric type elem_t. Epsilon is defined as the
+ * smallest number for which the sum of one and that number is different of
+ * one. It is assumed that the underlying representation of elem_t uses
+ * base two.
+ */
+static elem_t epsilon()
+{
+ elem_t eps;
+ for (eps = 1; 1 + eps != 1; eps /= 2)
+ ;
+ return 2 * eps;
+}
+
+int main(int argc, char** argv)
+{
+ int matrix_size;
+ int silent;
+ elem_t *a, *inv, *prod;
+ elem_t eps;
+ double error;
+ double ratio;
+
+ matrix_size = (argc > 1) ? atoi(argv[1]) : 3;
+ s_nthread = (argc > 2) ? atoi(argv[2]) : 3;
+ silent = (argc > 3) ? atoi(argv[3]) : 0;
+
+ eps = epsilon();
+ a = new_matrix(matrix_size, matrix_size);
+ init_matrix(a, matrix_size, matrix_size);
+ inv = invert_matrix(a, matrix_size);
+ prod = multiply_matrices(a, matrix_size, matrix_size,
+ inv, matrix_size, matrix_size);
+ error = identity_error(prod, matrix_size);
+ ratio = error / (eps * matrix_size);
+ if (! silent)
+ {
+ printf("error = %g; epsilon = %g; error / (epsilon * n) = %g\n",
+ error, eps, ratio);
+ }
+ if (ratio < 100)
+ printf("Error within bounds.\n");
+ else
+ printf("Error out of bounds.\n");
+ delete_matrix(prod);
+ delete_matrix(inv);
+ delete_matrix(a);
+
+ return 0;
+}
-------------------------------------------------------------------------
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse0120000070mrt/direct/01/
_______________________________________________
Valgrind-developers mailing list
Valgrind-developers@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/valgrind-developers
