Subject: convolution example

From: L J Clarke <[email-masked]>
Date: Wed, 16 Jun 93 13:07:09 BST

Just for completeness here is David's convolution example cleaned up a
bit as promised.  I've really done little more than add a few comments
and localised the manipulation of communicators around the usage
thereof. Hope I didnt break it!

Best Wishes
Lyndon

----------------------------------------------------------------------
/* 
 *
 *  PROGRAM 2.
 *  Now we'll try something a little harder. In the next program we find the
 *  convolution, C, of two vectors, A and B, by first finding the FFT of A
 *  and B, evaluating their elementwise product in Fourier space, and then
 *  finding the inverse transform, i.e.,
 *  
 *  C = INVFFT ( FFT(A) * FFT(B) )
 *
 *  Of course, we can do this by having all processes first find FFT(A) and
 *  then FFT(B). The product in Fourier space can then be found with no
 *  communication being necessary. Finally, all processes cooperate to
 *  perform the inverse FFT. Program 2 doesn't do this. Instead, we split
 *  the initial group into two equally-sized subgroups, one of which
 *  evaluates FFT(A) and the other FFT(B) in parallel. Communication is
 *  then necessary between the initial group and the subgroups to 
 *  re-distribute the data before doing the Fourier space product in the
 *  initial group. Finally, the processes in the initial group evaluate
 *  the inverse FFT, giving the convolution, C. The point here is that
 *  communication between groups is required.
 *
 */

main ()
{
   group_t init_group, fft_group;
   comm_t  comm_fft, comm_merge_A, comm_merge_B, comm_inv;
   comm_handle_t handle1, handle2;
   rso_handle_t  rso;
   int     init_rank, init_size, key, fft_rank, fft_size,
           fftlen, local_len, nbytes;
   int     source_rank, dest_rank1, dest_rank2;
   double  F[100], A[100], B[100], C[100];

   /* input the fft length */
   ret = get_input2 (fftlen);

   /* find out about the initial group */
   init_group = mpi_initial_group ();
   init_rank  = mpi_rank (init_group);
   init_size  = mpi_group_size (init_group);

   /* split the group in half, can be done any old way */
   key = init_rank&1;
   ret = mpi_split_group (init_group, key, init_rank, fft_group);

   /* find out about this fft group */
   fft_rank = mpi_rank (fft_group);
   fft_size = mpi_rank (fft_group);

   /* work out how big merge messages will be */
   local_len  = fftlen/init_size;
   nbytes     = local_len*sizeof(double);

   /* work out sources and destinations of merge messages */
   source_rank = fft_rank;
   dest_rank1  = 2*source_rank;
   dest_rank2  = dest_rank1 + 1;

   /* do the fft */
   ret = mpi_safemake_comm (fft_group,  fft_group,  comm_fft);
   ret = initialize_fft (F, fftlen, comm_fft);
   ret = parallel_fft   (F, fftlen, comm_fft);
   ret = mpi_free(comm_fft);
   /* done */

   /* do the merge */
   if (key != 0) 
   {
      ret = mpi_safemake_comm (fft_group,  init_group, comm_merge_A);
      ret = mpi_isendc (handle1, &F[0], nbytes, MPI_DOUBLE,
                              0, dest_rank1, comm_merge_A);
      ret = mpi_isendc (handle2, &F[local_len], nbytes, MPI_DOUBLE,
                                      0, dest_rank2, comm_merge_A);
      ret = mpi_publish (comm_merge_A, "GROUP_A");
      ret = mpi_subscribe ("GROUP_B", comm_merge_B);
   }
   else 
   {
      ret = mpi_safemake_comm (fft_group,  init_group, comm_merge_B);
      ret = mpi_isendc (handle1, &F[0], nbytes, MPI_DOUBLE,
                              0, dest_rank1, comm_merge_B);
      ret = mpi_isendc (handle2, &F[local_len], nbytes, MPI_DOUBLE,
                                      0, dest_rank2, comm_merge_B);
      ret = mpi_publish (comm_merge_B, "GROUP_B");
      ret = mpi_subscribe ("GROUP_A", comm_merge_A);
   }
   ret = mpi_recvc (A, nbytes, MPI_DOUBLE, MPI_DONTCARE, 
                        source_rank, comm_merge_A, rso);
   ret = mpi_recvc (B, nbytes, MPI_DOUBLE, MPI_DONTCARE, 
                        source_rank, comm_merge_B, rso);
   ret = mpi_wait (handle1, rso);
   ret = mpi_wait (handle2, rso);
   ret = mpi_free(comm_merge_A);
   ret = mpi_free(comm_merge_B);
   /* done */

   /* do the product */
   ret = parallel_product (A, B, C, fftlen);
   /* done */

   /* do the inverse fft and output result */
   ret = mpi_safemake_comm (init_group, init_group, comm_inv);
   ret = parallel_invfft  (C, fftlen, comm_inv);
   ret = output_result    (C, fftlen, comm_inv); 
   ret = mpi_free(comm_inv);
   /* done */

   /* finished */
   exit (0);

}

----------------------------------------------------------------------

         /--------------------------------------------------------\
    e||) | Lyndon J Clarke    Edinburgh Parallel Computing Centre | e||) 
    c||c | Tel: 031 650 5021  Email: [email-masked]  | c||c 
         \--------------------------------------------------------/