[MPI] Solution 5 : 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 process stores \(\frac{N}{p}\) rows, with \(p\) the number of processes.

Hints

Use the functions alloc_matrix and free_matrix to create and destroy a matrix with arbitrary dimensions.

Use a global reduction operation to implement a parallel dot product.

What performance do you get?

Cole este comando no seu terminal:

A senha para conectar é

Questão 1: Gram Schmidt Algorithm

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


/* ---------------------------------------------------------------------
 * alloc_matrix()
 * ------------------------------------------------------------------ */


double**                  alloc_matrix

  ( 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;
}


/* ---------------------------------------------------------------------
 * free_matrix()
 * ------------------------------------------------------------------ */


void                      free_matrix

  ( double**                mat )

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


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


int                       main

  ( int                     argc,
    char**                  argv )

{
  const int  N = 1000;
  double**   a = alloc_matrix ( N, N );

  double     c[N];
  double     d[N];
  double     t;
  double     wt0;
  int        my_rank, proc_count;
  int        i, j, k, m;


  MPI_Init      ( &argc, &argv );
  MPI_Comm_rank ( MPI_COMM_WORLD, &my_rank );
  MPI_Comm_size ( MPI_COMM_WORLD, &proc_count );

  if ( proc_count > N )
  {
    fprintf ( stderr, "Too many processes.\n" );
    return 1;
  }

  wt0 = MPI_Wtime ();
  m   = N / proc_count;

  if ( my_rank == 0 )
  {
    printf ( "Initializing the matrix ...\n" );
  }

  for ( j = 0; j < N; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      a[i][j] = (double) abs ( (my_rank * m + i) - j );
    }
  }

  if ( my_rank == 0 )
  {
    printf ( "Executing the Gram-Schmidt algorithm ...\n" );
  }

  for ( j = 0; j < N; j++ )
  {
    for ( k = 0; k < j; k++ )
    {
      c[k] = 0.0;

      for ( i = 0; i < m; i++ )
      {
        c[k] += a[i][k] * a[i][j];
      }
    }

    MPI_Allreduce ( c, d, j, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );

    for ( k = 0; k < j; k++ )
    {
      for ( i = 0; i < m; i++ )
      {
        a[i][j] -= d[k] * a[i][k];
      }
    }

    t = 0.0;

    for ( i = 0; i < m; i++ )
    {
      t += a[i][j] * a[i][j];
    }

    MPI_Allreduce ( &t, d, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );

    t = sqrt ( d[0] );

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

    t = 1.0 / t;

    for ( i = 0; i < m; i++ )
    {
      a[i][j] *= t;
    }
  }

  if ( my_rank == 0 )
  {
    printf ( "Ready in %g seconds.\n", (MPI_Wtime() - wt0) );
  }

  MPI_Finalize ();

  free_matrix ( a );
  return 0;
}

Questão 2: Number of MPI processes

Specify the number of MPI processes (1-64) that should be used for running your program (mpirun -np <nproc>)

Questão 3: 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.

Autores	Erik-Jan Lingen, Matthias Möller
Prazo de entrega	Sem prazo
Limite de submissão	No limitation

Informações

Entrar