wrapperMPI.cpp

Go to the documentation of this file.

/*
HPhi  -  Quantum Lattice Model Simulator
Copyright (C) 2015 The University of Tokyo

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifdef __MPI
#include <mpi.h>
#endif
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "wrapperMPI.hpp"
#ifdef _OPENMP
#include <omp.h>
#endif
#include <cmath>
#include <complex>
#include "splash.hpp"
#include "global.hpp"
#include "common/setmemory.hpp"

void InitializeMPI(int argc, char *argv[]){
#ifdef __MPI
  int ierr;
  ierr = MPI_Init(&argc, &argv);
  ierr = MPI_Comm_size(MPI_COMM_WORLD, &nproc);
  ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
  if(ierr != 0) exitMPI(ierr);
#else
  nproc = 1;
  myrank = 0;
#endif
  if (myrank == 0) stdoutMPI = stdout;
  else stdoutMPI = fopen("/dev/null", "w");
  splash();

#pragma omp parallel default(none) shared(nthreads)
#pragma omp master
#ifdef _OPENMP
  nthreads = omp_get_num_threads();
#else
  nthreads=1;
#endif
  fprintf(stdoutMPI, "\n\n#####  Parallelization Info.  #####\n\n");
  fprintf(stdoutMPI, "  OpenMP threads : %d\n", nthreads);
  fprintf(stdoutMPI, "  MPI PEs : %d \n\n", nproc);
}/*void InitializeMPI(int argc, char *argv[])*/
void FinalizeMPI(){
#ifdef __MPI
  int ierr;
  ierr = MPI_Finalize();
  if (ierr != 0) fprintf(stderr, "\n  MPI_Finalize() = %d\n\n", ierr);
#endif
  if (myrank != 0) fclose(stdoutMPI);
}
void exitMPI(
  int errorcode
)
{
  fflush(stdout);
#ifdef __MPI
  fprintf(stdout,"\n\n #######  [HPhi] You DO NOT have to WORRY about the following MPI-ERROR MESSAGE.  #######\n\n");
  int ierr;
  ierr = MPI_Abort(MPI_COMM_WORLD, errorcode);
  ierr = MPI_Finalize();
  if (ierr != 0) fprintf(stderr, "\n  MPI_Finalize() = %d\n\n", ierr);
#endif
  exit(errorcode);
}/*void exitMPI*/
FILE* fopenMPI(
  const char* FileName,
  const char* mode
){
  FILE* fp;

  if (myrank == 0) fp = fopen(FileName, mode);
  else fp = fopen("/dev/null", "w");

  return fp;
}/*FILE* fopenMPI*/
char* fgetsMPI(
  char* InputString,
  int maxcount,
  FILE* fp
){
  int inull;
  char *ctmp;

  ctmp = InputString;
  inull = 0;
  if (myrank == 0) {
    ctmp = fgets(InputString, maxcount, fp);
    if (ctmp == NULL){
      inull = 1;
    }
    
    while(*InputString == '\n' || strncmp(InputString, "#", 1)==0){
      ctmp = fgets(InputString, maxcount, fp);
      if (ctmp == NULL){
        inull=1;
        break;
      }
    }
  }
#ifdef __MPI
  MPI_Bcast(InputString, maxcount, MPI_CHAR, 0, MPI_COMM_WORLD);
  MPI_Bcast(&inull, 1, MPI_INT, 0, MPI_COMM_WORLD);
#endif
  if (myrank != 0 && inull == 1) {
    ctmp = NULL;
  }

  return ctmp;
}/*char* fgetsMPI*/
void BarrierMPI(){
#ifdef __MPI
  MPI_Barrier(MPI_COMM_WORLD);
#endif
}/*void BarrierMPI()*/
long int MaxMPI_li(
  long int idim
){
#ifdef __MPI
  int ierr;
  ierr = MPI_Allreduce(MPI_IN_PLACE, &idim, 1,
    MPI_LONG, MPI_MAX, MPI_COMM_WORLD);
  if(ierr != 0) exitMPI(-1);
#endif
  return(idim);
}/*long int MaxMPI_li*/
double MaxMPI_d(
  double dvalue
){
#ifdef __MPI
  int ierr;
  ierr = MPI_Allreduce(MPI_IN_PLACE, &dvalue, 1,
    MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
  if(ierr != 0) exitMPI(-1);
#endif
  return(dvalue);
}/*double MaxMPI_d*/
std::complex<double> SumMPI_dc(
  std::complex<double> norm
){
#ifdef __MPI
  int ierr;
  ierr = MPI_Allreduce(MPI_IN_PLACE, &norm, 1,
    MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
  if(ierr != 0) exitMPI(-1);
#endif
  return(norm);
}/*std::complex<double> SumMPI_dc*/
double SumMPI_d(
  double norm
){
#ifdef __MPI
  int ierr;
  ierr = MPI_Allreduce(MPI_IN_PLACE, &norm, 1,
    MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD);
  if(ierr != 0) exitMPI(-1);
#endif
  return(norm);
}/*double SumMPI_d*/
void SumMPI_dv(
  int nnorm,
  double *norm
) {
#ifdef __MPI
  int ierr;
  ierr = MPI_Allreduce(MPI_IN_PLACE, norm, nnorm,
    MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD);
  if (ierr != 0) exitMPI(-1);
#endif
}/*void SumMPI_dv*/
void SumMPI_cv(
  int nnorm,
  std::complex<double> *norm
) {
#ifdef __MPI
  int ierr;
  ierr = MPI_Allreduce(MPI_IN_PLACE, norm, nnorm,
    MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
  if (ierr != 0) exitMPI(-1);
#endif
}/*void SumMPI_cv*/
long int SumMPI_li(
  long int idim
){
#ifdef __MPI
  int ierr;
  ierr = MPI_Allreduce(MPI_IN_PLACE, &idim, 1,
    MPI_LONG, MPI_SUM, MPI_COMM_WORLD);
  if(ierr != 0) exitMPI(-1);
#endif
  return(idim);
}/*long int SumMPI_li*/
int SumMPI_i(
  int idim
) {
#ifdef __MPI
  int ierr;
  ierr = MPI_Allreduce(MPI_IN_PLACE, &idim, 1,
                       MPI_INT, MPI_SUM, MPI_COMM_WORLD);
  if(ierr != 0) exitMPI(-1);
#endif
  return(idim);
}/*int SumMPI_i*/
long int BcastMPI_li(
  int root,
  long int idim
) {
  long int idim0;
  idim0 = idim;
#ifdef __MPI
    MPI_Bcast(&idim0, 1, MPI_LONG, root, MPI_COMM_WORLD);
#endif
  return(idim0);
}/*long int BcastMPI_li*/
double NormMPI_dc(
  long int idim,
  std::complex<double> *_v1
){
  double dnorm =0;
  long int i;

  dnorm = 0.0;
#pragma omp parallel for default(none) private(i) \
shared(_v1, idim) reduction(+:dnorm)
    for (i = 1; i <= idim; i++) 
      dnorm += real(conj(_v1[i])*_v1[i]);
 
#ifdef __MPI
  dnorm = SumMPI_d(dnorm);
#endif
  return dnorm;
}/*double NormMPI_dc*/
void NormMPI_dv(
  long int ndim,
  int nstate,
  std::complex<double> **_v1,
  double *dnorm
) {
  long int idim;
  int istate;

  for (istate = 0; istate < nstate; istate++) dnorm[istate] = 0.0;
  for (idim = 1; idim <= ndim; idim++) {
    for (istate = 0; istate < nstate; istate++) {
      dnorm[istate] += real(conj(_v1[idim][istate])*_v1[idim][istate]);
    }
  }
  SumMPI_dv(nstate, dnorm);
  for (istate = 0; istate < nstate; istate++) dnorm[istate] = sqrt(dnorm[istate]);
}/*double NormMPI_cv*/
std::complex<double> VecProdMPI(
  long int ndim,
  std::complex<double> *v1,
  std::complex<double> *v2
){
  long int idim;
  std::complex<double> prod, *prod_thr;
  int mythread;

  prod_thr = cd_1d_allocate(nthreads);
#pragma omp parallel default(none) shared(v1,v2,ndim,prod,prod_thr) private(idim,mythread)
  {
#ifdef _OPENMP
    mythread = omp_get_thread_num();
#else
    mythread = 0;
#endif
#pragma omp for
    for (idim = 1; idim <= ndim; idim++) 
      prod_thr[mythread] += conj(v1[idim]) * v2[idim];
  }
  prod = 0.0;
  for (mythread = 0; mythread < nthreads; mythread++)
    prod += prod_thr[mythread];
  free_cd_1d_allocate(prod_thr);

  prod = SumMPI_dc(prod);

  return(prod);
}/*std::complex<double> VecProdMPI*/
void MultiVecProdMPI(
  long int ndim,
  int nstate,
  std::complex<double> **v1,
  std::complex<double> **v2,
  std::complex<double> *prod
) {
  long int idim;
  int istate;

  for (istate = 0; istate < nstate; istate++) prod[istate] = 0.0;
  for (idim = 1; idim <= ndim; idim++) {
    for (istate = 0; istate < nstate; istate++) {
      prod[istate] += conj(v1[idim][istate])*v2[idim][istate];
    }
  }
  SumMPI_cv(nstate, prod);
}/*void MultiVecProdMPI*/
void SendRecv_cv(
  int origin,
  long int nMsgS,
  long int nMsgR,
  std::complex<double> *vecs,
  std::complex<double> *vecr
) {
#ifdef __MPI
  int ierr, two31m1 = 2147483647, modMsg, nMsgS2, nMsgR2;
  long int nMsg, nnMsg, iMsg, sMsgR, sMsgS;
  MPI_Status statusMPI;

  if (nMsgS > nMsgR) nMsg = nMsgS;
  else nMsg = nMsgR;
  nnMsg = nMsg / two31m1;
  modMsg = nMsg % two31m1;
  if (modMsg != 0) nnMsg += 1;

  sMsgS = 0;
  sMsgR = 0;
  for (iMsg = 0; iMsg < nnMsg; iMsg++) {
    nMsgS2 = nMsgS / nnMsg;
    nMsgR2 = nMsgR / nnMsg;
    if (iMsg < nMsgS % nnMsg) nMsgS2 += 1;
    if (iMsg < nMsgR % nnMsg) nMsgR2 += 1;

    ierr = MPI_Sendrecv(&vecs[sMsgS], nMsgS2, MPI_DOUBLE_COMPLEX, origin, 0,
                        &vecr[sMsgR], nMsgR2, MPI_DOUBLE_COMPLEX, origin, 0,
                        MPI_COMM_WORLD, &statusMPI);
    if (ierr != 0) exitMPI(-1);

    sMsgS += nMsgS2;
    sMsgR += nMsgR2;
  }
#endif
}/*void SendRecv_cv*/
void SendRecv_iv(
  int origin,
  long int nMsgS,
  long int nMsgR,
  long int *vecs,
  long int *vecr
) {
#ifdef __MPI
  int ierr, two31m1 = 2147483647, modMsg, nMsgS2, nMsgR2;
  long int nMsg, nnMsg, iMsg, sMsgR, sMsgS;
  MPI_Status statusMPI;

  if (nMsgS > nMsgR) nMsg = nMsgS;
  else nMsg = nMsgR;
  nnMsg = nMsg / two31m1;
  modMsg = nMsg % two31m1;
  if (modMsg != 0) nnMsg += 1;

  sMsgS = 0;
  sMsgR = 0;
  for (iMsg = 0; iMsg < nnMsg; iMsg++) {
    nMsgS2 = nMsgS / nnMsg;
    nMsgR2 = nMsgR / nnMsg;
    if (iMsg < nMsgS % nnMsg) nMsgS2 += 1;
    if (iMsg < nMsgR % nnMsg) nMsgR2 += 1;

    ierr = MPI_Sendrecv(&vecs[sMsgS], nMsgS2, MPI_LONG, origin, 0,
                        &vecr[sMsgR], nMsgR2, MPI_LONG, origin, 0,
                        MPI_COMM_WORLD, &statusMPI);
    if (ierr != 0) exitMPI(-1);

    sMsgS += nMsgS2;
    sMsgR += nMsgR2;
  }
#endif
}/*void SendRecv_iv*/
long int SendRecv_i(
  int origin,
  long int isend
) {
#ifdef __MPI
  int ierr;
  MPI_Status statusMPI;
  long int ircv;
  ierr = MPI_Sendrecv(&isend, 1, MPI_LONG, origin, 0,
                      &ircv,  1, MPI_LONG, origin, 0,
                      MPI_COMM_WORLD, &statusMPI);
  if (ierr != 0) exitMPI(ierr);
  return ircv;
#else
  return isend;
#endif
}/*void SendRecv_i*/