[Threads] Solution 8 : Gram-Schmidt Algorithm

This material is part of the course Practical Introduction to Parallel Programming and meant for educational purposes only.

Task

Implement a parallel program that executes the Gram-Schmidt algorithm to make the columns of an \(N \times N\) matrix \(\mathbf{A}\) orthonormal through a series of dot products and vector subtractions.

The pseudo-code of the Gram-Schmidt algorithm reads as follows:

\begin{align*} &\text{do }j=1,\dots,N\\ &\quad \text{do }k=1,\dots,j-1\\ &\quad\quad c_k := \mathbf{A}_j^\top \cdot \mathbf{A}_k\\ &\quad \text{end do}\\ &\quad \mathbf{A}_j := \mathbf{A}_j - \sum_{k=1}^{j-1} c_k \cdot \mathbf{A}_k\\ &\quad \mathbf{A}_j := \frac{\mathbf{A}_j}{\| \mathbf{A}_j \|}\\ &\text{end do} \end{align*}

Note that \(\mathbf{A}_{j}\) denotes the \(j\)-th column of the matrix \(\mathbf{A}\) and superscript \(\top\) stands for taking the transpose.

The code below show a serial implementation of the Gram-Schmidt algorithm.

In the parallel implementation of the Gram-Schmidt algorithm, the matrix \(\mathbf{A}\) is distributed row-wise over the processors:

This means that each thread handles \(\frac{N}{p}\) rows, with \(p\) the number of threads.

Hints

Use the functions matrix_alloc and matrix_free to create and destroy a matrix with arbitrary dimensions.

Define a struct data_t that is to be passed to each thread. It should contain at least the matrix size, a pointer to the matrix, a pointer to a shared reduction variable, and the first and last row indices to be processed by a thread.

You must define a data_t instance for each thread!

What performance do you get? How could you improve the performance?

Paste this command into your terminal:

The password to connect is

Question 1: Gram-Schmidt Algorithm

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <pthread.h>

/* ---------------------------------------------------------------------
 * types and constants
 * ------------------------------------------------------------------ */

#define                   NTHREADS  5

struct                    reduction
{
  pthread_mutex_t           mutex;
  pthread_cond_t            cond;
  int                       nthreads;
  int                       cycle;
  int                       tcount[2];
  double                    result[2];
};

struct                    data
{
  int                       msize;
  int                       ifirst;
  int                       ilast;
  double**                  matrix;
  struct reduction*         red;
};

typedef struct reduction  reduction_t;
typedef struct data       data_t;

/* ---------------------------------------------------------------------
 * reduction_init()
 *
 *   Initialises a reduction variable. The parameter "nt" specifies
 *   the number of threads that will participate in a reduction
 *   operation.
 * ------------------------------------------------------------------ */

void                      reduction_init

  ( reduction_t*            red,
    int                     nt )

{
  assert ( red && (nt > 0) );

  pthread_mutex_init ( &red->mutex, NULL );
  pthread_cond_init  ( &red->cond,  NULL );

  red->nthreads  = nt;
  red->cycle     = 0;
  red->tcount[0] = 0;
  red->tcount[1] = 0;
  red->result[0] = 0.0;
  red->result[1] = 0.0;
}

/* ---------------------------------------------------------------------
 * reduction_destroy()
 *
 *   Destroys a reduction variable.
 * ------------------------------------------------------------------ */

void                      reduction_destroy

  ( reduction_t*            red )

{
  assert ( red );

  pthread_mutex_destroy ( &red->mutex );
  pthread_cond_destroy  ( &red->cond );

  red->nthreads = 0;
}

/* ---------------------------------------------------------------------
 * reduction_sum()
 *
 *   Calculates a global sum across multiple threads. The number
 *   of threads calling this function must be equal to the
 *   parameter "nt" specified in the call to the reduction_init()
 *   function.
 * ------------------------------------------------------------------ */

double                    reduction_sum

  ( reduction_t*            red,
    double                  val )

{
  assert ( red && (red->nthreads > 0) );

  /* Store the current cycle in a local variable. This is essential
   * as the last thread will change the cycle. */

  int     cycle = red->cycle;
  double  result;

  pthread_mutex_lock ( &red->mutex );

  if ( red->tcount[cycle] == 0 )
  {
    /* This is the first thread. */

    red->result[cycle]  = val;
  }
  else
  {
    red->result[cycle] += val;
  }

  red->tcount[cycle]++;

  if ( red->tcount[cycle] < red->nthreads )
  {
    /* Wait for the last thread. */

    while ( red->tcount[cycle] < red->nthreads )
    {
      pthread_cond_wait ( &red->cond, &red->mutex );
    }
  }
  else
  {
    /* This is the last thread. Prepare the next cycle and wake the
     * other threads. */

    red->cycle              = (red->cycle + 1) % 2;
    red->tcount[red->cycle] = 0;

    pthread_cond_broadcast ( &red->cond );
  }

  result = red->result[cycle];

  pthread_mutex_unlock ( &red->mutex );

  return result;
}

/* ---------------------------------------------------------------------
 * matrix_alloc()
 *
 *   Allocates an m * n matrix. The matrix must be deallocated by
 *   calling the function matrix_free()
 * ------------------------------------------------------------------ */

double**                  matrix_alloc

  ( int                     m,
    int                     n )

{
  double**  rows = NULL;
  double*   data = NULL;

  int       i;

  assert ( (m >= 0) && (n >= 0) );

  if ( (m == 0) || (n == 0) )
  {
    return rows;
  }

  rows = (double**) malloc ( m *     sizeof(double*) );
  data = (double*)  malloc ( m * n * sizeof(double)  );

  if ( ! rows || ! data )
  {
    fprintf ( stderr, "Out of memory.\n" );
    exit    ( 1 );
  }

  for ( i = 0; i < m; i++ )
  {
    rows[i] = &data[i * n];
  }

  return rows;
}

/* ---------------------------------------------------------------------
 * matrix_free()
 *
 *   Deallocates a matrix that has been allocated by the function
 *   matrix_alloc().
 * ------------------------------------------------------------------ */

void                      matrix_free

  ( double**                mat )

{
  if ( mat )
  {
    free ( mat[0] );
    free ( mat );
  }
}

/* ---------------------------------------------------------------------
 * thread_func()
 *
 *   The thread entry point. This executes the Gram-Schmidt
 *   algorithm in parallel.
 * ------------------------------------------------------------------ */

void*                     thread_func

  ( void*                   arg )

{
  data_t*       data   = (data_t*) arg;
  reduction_t*  red    = data->red;
  double**      matrix = data->matrix;
  double*       coeff  = NULL;
  const int     msize  = data->msize;
  const int     ifirst = data->ifirst;
  const int     ilast  = data->ilast;

  double        t;
  int           i, j, k;


  coeff = (double*) malloc ( msize * sizeof(*coeff) );

  for ( j = 0; j < msize; j++ )
  {
    for ( k = 0; k < j; k++ )
    {
      coeff[k] = 0.0;

      for ( i = ifirst; i < ilast; i++ )
      {
        coeff[k] += matrix[i][k] * matrix[i][j];
      }

      coeff[k] = reduction_sum ( red, coeff[k] );
    }

    for ( k = 0; k < j; k++ )
    {
      for ( i = ifirst; i < ilast; i++ )
      {
        matrix[i][j] -= coeff[k] * matrix[i][k];
      }
    }

    t = 0.0;

    for ( i = ifirst; i < ilast; i++ )
    {
      t += matrix[i][j] * matrix[i][j];
    }

    t = sqrt ( reduction_sum( red, t ) );

    if ( t <= 0.0 )
    {
      fprintf ( stderr, "Singular matrix.\n" );
      exit    ( 1 );
    }

    t = 1.0 / t;

    for ( i = ifirst; i < ilast; i++ )
    {
      matrix[i][j] *= t;
    }
  }

  free ( coeff );

  return NULL;
}

/* ---------------------------------------------------------------------
 * main()
 * ------------------------------------------------------------------ */

int                       main

  ( int                     argc,
    char**                  argv )

{
  const int    msize  = 200;
  double**     matrix = matrix_alloc ( msize, msize );

  pthread_t    threads[NTHREADS];
  data_t       data   [NTHREADS];
  reduction_t  red;
  double       t;
  int          i, j, k;


  reduction_init ( &red, NTHREADS );

  printf ( "Initialising the matrix ...\n" );

  for ( j = 0; j < msize; j++ )
  {
    for ( i = 0; i < msize; i++ )
    {
      matrix[i][j] = (double) abs ( i - j );
    }
  }

  printf ( "Executing the Gram-Schmidt algorithm ...\n" );

  j = 0;

  for ( i = 0; i < NTHREADS; i++ )
  {
    /* The number of rows to be processed by this thread: */

    k = (msize - j) / (NTHREADS - i);

    data[i].msize  = msize;
    data[i].ifirst = j;
    data[i].ilast  = j + k;
    data[i].red    = &red;
    data[i].matrix = matrix;

    pthread_create ( &threads[i], NULL, thread_func, &data[i] );

    j += k;
  }

  assert ( j == msize );

  for ( i = 0; i < NTHREADS; i++ )
  {
    pthread_join ( threads[i], NULL );
  }

  printf ( "Done.\n" );

  /* Check whether the first and last columns are
   * indeed orthogonal. */

  t = 0.0;

  for ( i = 0; i < msize; i++ )
  {
    t += matrix[i][0] * matrix[i][msize - 1];
  }

  if ( t < 1.0e-10 )
  {
    printf ( "Check (%g) passed.\n", t );
  }
  else
  {
    printf ( "Check (%g) FAILED.\n", t );
  }

  reduction_destroy ( &red );
  matrix_free       ( matrix );

  return 0;
}

Question 2: Command Line Argument

Specify the command line arguments that should be passed to your program when it is run. Arguments must be given in quotes and separated by commas. Leave this filed empty if you do not want to specify command line arguments.

Example: "1","int","arg=2" is interpreted as three arguments 1, int, and arg=2.

Author(s)	Erik-Jan Lingen, Matthias Möller
Deadline	No deadline
Submission limit	No limitation

Information

Sign in