mpi_gs.cpp

Go to the documentation of this file.

#undef USE_ROS

#ifdef USE_ROS
#include "ros/ros.h"
#endif

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include <Eigen3/Array>
#include <Eigen3/SVD>
#include <Eigen3/LU>
#include <iostream>
#include <sstream>

void get_max_avg(unsigned int m, double* delta, double &max, double &avg)
{
    avg = 0;
    for (unsigned int j=0; j<m ; j++)
    {
      if (delta[j] > max) max = delta[j];
      avg += delta[j];
    }
    avg = avg/(double)m;
}

void mat_print(FILE* file,unsigned int m,unsigned int n, double *a)
{
  if (!file) return;
  // print matrix
  for (unsigned int j=0;j<n;j++)
  {
    for (unsigned int i=0;i<m;i++) fprintf(file,"%12.5f\t",a[i+j*m]);
    fprintf(file,"\n");
  }
}

void mat_print(FILE* file,Eigen3::MatrixXf a)
{
  if (!file) return;
  unsigned int m = a.cols();
  unsigned int n = a.rows();
  // print matrix
  for (unsigned int j=0;j<n;j++)
  {
    for (unsigned int i=0;i<m;i++) fprintf(file,"%12.5f\t",a(i,j));
    fprintf(file,"\n");
  }
}

void mat_print(FILE* file,Eigen3::VectorXf v)
{
  if (!file) return;
  unsigned int m = v.size();
  // print vector
  for (unsigned int i=0;i<m;i++) fprintf(file,"%12.5f\t",v(i));
  fprintf(file,"\n");
}

void gs(unsigned int m, double* a, double* x, double* b,unsigned int* order,double* rq ,double* delta,FILE* file)
{
  // performs one iteration of gs
  for (unsigned int i = 0; i < m; i++)
  {
    unsigned int index = order[i];
    delta[index] = b[index]/a[index+index*m];
    for (unsigned int j = 0; j < m; j++)
    {
      delta[index] -= a[j+index*m]*x[j]/a[index+index*m];
    }
    x[index] += delta[index];
    if (delta[index] == 0)
      rq[index] = 0; //-1.0*(double)(i+1); //0.5/(double)(i+1);
    else
      rq[index] = fabs(delta[index])/(double)(i+1); // -1.0*(double)i; //

    rq[index] = (double)i;
  }
/*
  if (file)
  {
    fprintf(file,"x: ");
    mat_print(file,m,1,x);
    fprintf(file,"d: ");
    mat_print(file,m,1,delta);
    fprintf(file,"r: ");
    mat_print(file,m,1,rq);
  }
*/
}

int main(int argc, char *argv[]) {
  int numprocs, rank, namelen;
  char processor_name[MPI_MAX_PROCESSOR_NAME];

  MPI_Init(&argc, &argv);
  //MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  //MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
  rank = MPI::COMM_WORLD.Get_rank();
  numprocs = MPI::COMM_WORLD.Get_size();

  MPI_Get_processor_name(processor_name, &namelen);

  std::ostringstream log_file_stream;
  log_file_stream << "/tmp/mpi_gs_" << rank << ".x";
  FILE* log_file;
  log_file = fopen(log_file_stream.str().c_str(),"w");
  if (log_file) fprintf(log_file,"Process %d on %s out of %d\n", rank, processor_name, numprocs);
  printf("Process %d on %s out of %d\n", rank, processor_name, numprocs);

  // init random seed
  srand(time(NULL));

  // construct a mxm square matrix, note, this should be consistent across all mpi nodes
  const double matrix_row_sum_coeff = 2.0;
  const unsigned int m = 1000;
  const int num_iterations = 50;
  double* mat_a = new double[m*m];
  double* mat_x = new double[m];
  double* mat_b = new double[m];
  double* mat_rq = new double[m];
  double* mat_delta = new double[m];
  unsigned int* mat_order = new unsigned int[m]; // rearrange order of matrix

  for (unsigned int j=0;j<m;j++)
  {
    mat_rq[j] = 0; // recursive quality
    mat_order[j] = j; // access order
    mat_b[j] = (double)rand()/(double)RAND_MAX; // rhs
    for (unsigned int i=0;i<m;i++)
    {
      mat_x[i] = 0; // solution, preset to 0
      mat_a[i+j*m] = matrix_row_sum_coeff*(double)rand()/(double)RAND_MAX/(double)m; // strides in i first
    }
    if (j<m) mat_a[j+j*m] = 1.0; // ensure diagonal of A is non-zero, set to 1
  }

#if 1
  // shift n-times based on rank
  for (int s=0;s<rank;s++) // shift rank-times
  {
    unsigned int tmp = mat_order[0];
    for (unsigned int j=1;j<m;j++) mat_order[j-1] = mat_order[j];
    mat_order[m-1] = tmp;
  }

  // shuffle
  for (unsigned int j=1;j<m;j++)
  {
    unsigned int tmpj = rand() % j; // between 0 and j-1
    //if (rank==20) printf("swap %d and %d\n",j,tmpj);
    unsigned int tmp = mat_order[j];
    mat_order[j] = mat_order[tmpj];
    mat_order[tmpj] = tmp;
  }
#endif

  // write A matrix
  if (rank == 0)
  {
    std::ostringstream mat_a_file_stream;
    mat_a_file_stream << "/tmp/mpi_gs_A.x";
    FILE* mat_a_file;
    mat_a_file = fopen(mat_a_file_stream.str().c_str(),"w");
    if (mat_a_file) fprintf(mat_a_file,"Matrix:\n");
    if (mat_a_file) mat_print(mat_a_file,m,m,mat_a);
  }

  double init_time;

  // receiving buffers
  MPI::Request x_send_request[numprocs];
  MPI::Request x_recv_request[numprocs];
  MPI::Request rq_send_request[numprocs];
  MPI::Request rq_recv_request[numprocs];
  MPI::Status x_send_status[numprocs];
  MPI::Status x_recv_status[numprocs];
  MPI::Status rq_send_status[numprocs];
  MPI::Status rq_recv_status[numprocs];
  double* buf_x = new double[numprocs*m];  // this fails on single machine with 21+ nodes
  double* buf_rq = new double[numprocs*m];  // this fails on single machine with 21+ nodes

  if (rank == 0) init_time = MPI::Wtime();
  for (int iteration=0;iteration<num_iterations;iteration++)
  {
    double gs_time;
    if (rank == 0) gs_time = MPI::Wtime();
    gs(m,mat_a,mat_x,mat_b,mat_order,mat_rq,mat_delta,log_file);
    if (rank == 0) if (log_file) fprintf(log_file,"rq(gs) %d time: %f",iteration,MPI::Wtime()-gs_time);

    // round of exchange based on rq
    // synchronize solutions across nodes with latest updates
    double comm_time;
    if (rank == 0) comm_time = MPI::Wtime();
    for (int i=0; i < numprocs; i++)
    {
      if (i != rank)
      {
        // send solution and rq to all other nodes
        x_send_request[i] = MPI::COMM_WORLD.Isend(mat_x, m, MPI_DOUBLE , i, 1000+1000*(iteration+1)+i);
        rq_send_request[i] = MPI::COMM_WORLD.Isend(mat_rq, m, MPI_DOUBLE , i, 2000+1000*(iteration+1)+i);
        // receive solution and rq from all other nodes
        x_recv_request[i] = MPI::COMM_WORLD.Irecv(&buf_x[i*m], m, MPI_DOUBLE , i, 1000+1000*(iteration+1)+rank);
        rq_recv_request[i] = MPI::COMM_WORLD.Irecv(&buf_rq[i*m], m, MPI_DOUBLE , i, 2000+1000*(iteration+1)+rank);
      }
    }
    for (int i=0; i < numprocs; i++)
    {
      if (i != rank)
      {
        x_send_request[i].Wait(x_send_status[i]);  // don't need to wait for send request as they are implicitly enforced by receivers
        rq_send_request[i].Wait(rq_send_status[i]);  // don't need to wait for send request as they are implicitly enforced by receivers

        x_recv_request[i].Wait(x_recv_status[i]);
        rq_recv_request[i].Wait(rq_recv_status[i]);
      }
    }
    if (rank == 0) if (log_file) fprintf(log_file," rq(comm) %d time: %f",iteration,MPI::Wtime()-comm_time);

    for (int i=0; i < numprocs; i++)
    {
      if (i != rank)
      {
        // update solution based on rq
        for (unsigned int j=0; j<m ; j++)
        {
          //if (buf_rq[i*m+j] < mat_rq[j])
          if (buf_rq[i*m+j] == m-1)
          {
            mat_x[j] = buf_x[i*m+j];
            mat_rq[j] = buf_rq[i*m+j];
            //mat_x[j] = 0.5*(mat_x[j] + buf_x[i*m+j]);
            //mat_rq[j] = 0.5*(mat_rq[j] + buf_rq[i*m+j]);
          }
        }
      }
    }
    double tmp_avg = 0;
    double tmp_max = 0;
    get_max_avg(m, mat_delta,tmp_max,tmp_avg);
    if (log_file) fprintf(log_file," rq convergence: %d max: %f avg:%f \n",iteration,tmp_max,tmp_avg);
  }
  if (rank == 0) if (log_file) fprintf(log_file,"rq time: %f\n",MPI::Wtime()-init_time);
  if (log_file) mat_print(log_file,m,1,mat_x);

  delete buf_x;
  delete buf_rq;

  // re compute without rq
  for (unsigned int j=0; j<m ; j++) mat_x[j] = 0.0; // reset solution and try again
  if (rank == 0) init_time = MPI::Wtime();
  for (int iteration=0;iteration<num_iterations;iteration++)
  {
    gs(m,mat_a,mat_x,mat_b,mat_order,mat_rq,mat_delta,log_file);

    double tmp_avg = 0;
    double tmp_max = 0;
    get_max_avg(m, mat_delta,tmp_max,tmp_avg);
    if (log_file) fprintf(log_file,"gs convergence: %d max: %f avg:%f \n",iteration,tmp_max,tmp_avg);

  }
  if (rank == 0) if (log_file) fprintf(log_file,"gs time: %f\n",MPI::Wtime()-init_time);
  if (log_file) mat_print(log_file,m,1,mat_x);

  // re compute with eigen3
  Eigen3::MatrixXf e_a; e_a = Eigen3::MatrixXf::Zero(m,m);
  Eigen3::VectorXf e_x; e_x = Eigen3::VectorXf::Zero(m);
  Eigen3::VectorXf e_b; e_b = Eigen3::VectorXf::Zero(m);
  for (unsigned int j=0;j<m;j++)
  {
    for (unsigned int i=0;i<m;i++)
    {
      //e_x(i) = mat_x[i];
      e_a(j,i) = mat_a[i+j*m];
    }
    e_b(j) = mat_b[j];
  }

  if (rank == 0)
  {
    init_time = MPI::Wtime();
    e_x = e_a.fullPivLu().solve(e_b);
    if (log_file) fprintf(log_file,"lu time: %d %f\n",rank,MPI::Wtime()-init_time);
    if (log_file) mat_print(log_file,e_x);
  }

  if (rank == 0)
  {
    init_time = MPI::Wtime();
    Eigen3::JacobiSVD<Eigen3::MatrixXf> svd(e_a, Eigen3::ComputeThinU | Eigen3::ComputeThinV);
    std::cout << "Its singular values are:" << std::endl << svd.singularValues() << std::endl;
    e_x = svd.solve(e_b);
    if (log_file) fprintf(log_file,"svd time: %d %f\n",rank,MPI::Wtime()-init_time);
    if (log_file) mat_print(log_file,e_x);
  }
  

  if (log_file) fclose(log_file);

  delete mat_a;
  delete mat_x;
  delete mat_b;
  delete mat_rq;
  delete mat_delta;
  delete mat_order;

#ifdef USE_ROS
  setenv("ROS_ROOT","/wg/stor5/wgsim/hsu/projects/cturtle/ros",1);
  setenv("ROS_MASTER_URI","http://wgsg0:11311",1);
  // test string buffer
  char uri_buf[256];

  std::string ros_master_uri;
  if (rank == 0)
  {
    ros::init(argc,argv,"mpi_node",ros::init_options::NoSigintHandler);
    ros_master_uri = ros::master::getURI();
    ros_master_uri = "http://akb:11311";
    //ROS_INFO("rank %d master uri %s",rank,ros_master_uri.c_str());
    for (int i=1; i < numprocs; i++)
      MPI::COMM_WORLD.Send(ros_master_uri.c_str(), ros_master_uri.size(), MPI_CHAR , i, 100);
    // ros::master::V_TopicInfo topics;
    // bool succ = ros::master::getTopics(topics);
    // printf("rank %d has %d topics %d",rank,(int)topics.size(),(int)succ);
    // printf("final: rank %d master uri %s",rank,ros::master::getURI().c_str());
  }
  else
  {
    MPI::COMM_WORLD.Recv(uri_buf, 256, MPI_CHAR , 0, 100);
    ros_master_uri = std::string(uri_buf);
  }
#endif




  MPI_Finalize();
}