Subject: Convolution still

From: David W. Walker <[email-masked]>
Date: Thu, 17 Jun 93 08:48:29 -0500

Dear Lyndon,
	The problem with forming the merge communicator arises because it must be
called by all processes with consistent arguments. Thus, all processes must call 
mpi_safemake_comm twice, once for each merge communicator. Please see the amended
version of Program 2 below. I've also changed the calls so they all return error codes.
This code would need to be further modified to cover
the general case of three arbitrary groups. Perhaps I'll do that when I have time -
perhaps while I'm sitting around in Rome airport next week waiting for some mad
European terrorist to blow me up.

David

----------------------------------------------------------------------
/* 
 *
 *  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, other_group;
   comm_t  comm_fft, comm_merge_A, comm_merge_B, comm_inv, other_comm;
   comm_handle_t handle1, handle2;
   rso_handle_t  rso;
   int     init_rank, init_size, key, fft_rank, fft_size,
           fftlen, local_len, nbytes, errno;
   int     source_rank, dest_rank1, dest_rank2;
   double  F[100], A[100], B[100], C[100];

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

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

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

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

   /* 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;

   /* create and publish communicator for FFTs */
   errno = mpi_safemake_comm (fft_group,  fft_group,  comm_fft);
   errno = mpi_publish (comm_fft, (key == 0 ? "even" : "odd"));

   /* do the FFT */
   errno = initialize_fft (F, fftlen, comm_fft);
   errno = parallel_fft   (F, fftlen, comm_fft);

   /* find out what the other FFT group's communicator is and extract group */
   errno = mpi_subscribe (other_comm, (key != 0 ? "even" : "odd"));
   errno = mpi_group (other_comm, other_group);
   /* free the FFT communicator */
   errno = mpi_free(comm_fft);

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

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

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

   /* finished */
   exit (0);