ascat.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2009 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar 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 2 of the License, or
    (at your option) any later version.

    magpar 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 magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: ascat.c 2681 2009-07-31 04:30:53Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "util.h"

/* Scatters and gathers integer data */
int ascat(int **aindata,int n_items,int s_item,int *newproc,int insize,int send)
{
  MagparFunctionInfoBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  assert(insize>=1);
  /* this function only works for distributing data from first processor
     to other processors
     TODO: round robin below does not work anyway; thus, simplify this function!
  */
  if (!rank) {
    assert(n_items>0);
  }
  else {
    assert(n_items==0);
  }

  int     *ta_n_items;          /* number of items */
  ierr = PetscMalloc(insize*sizeof(int),&ta_n_items);CHKERRQ(ierr);

  int     *ta_cnt;              /* counter for items */
  ierr = PetscMalloc(insize*sizeof(int),&ta_cnt);CHKERRQ(ierr);

  int     *ta_getn_items;       /* number of items to be received */
  ierr = PetscMalloc(insize*sizeof(int),&ta_getn_items);CHKERRQ(ierr);

  int     **ta_procdata;       /* pointer into data arrays */
  ierr = PetscMalloc(insize*sizeof(int*),&ta_procdata);CHKERRQ(ierr);

  int     *indata;
  indata=*aindata;

  /* allocate memory for temporary array of sorted data;
     at least one item to avoid malloc's complaints about allocating zero bytes
  */
  int *t_data;
  ierr = PetscMalloc(n_items*s_item*sizeof(int)+1,&t_data);CHKERRQ(ierr);

  /* reset counters */
  for (int i=0; i<insize; i++)
  {
    ta_n_items[i]=0;
    ta_cnt[i]=0;
  }

  /* count number of items to be sent */
  for (int i=0; i<n_items; i++)
  {
    assert(newproc[i]>=0 && newproc[i]<insize);
    ta_n_items[newproc[i]]++;
  }

  /* set pointers to start of data in temporary array */
  ta_procdata[0]=t_data;
  for (int i=1; i<insize; i++)
  {
    ta_procdata[i] = ta_procdata[i-1] + s_item*ta_n_items[i-1];
  }

  /* check pointers for consistency */
  assert(ta_procdata[insize-1]+s_item*ta_n_items[insize-1]==t_data+n_items*s_item);


  /* copy data sorted into temporary array */
  for (int i=0; i<n_items; i++)
  {
    for (int j=0; j<s_item; j++)
    {
      *(ta_procdata[newproc[i]]+ta_cnt[newproc[i]]*s_item+j)=indata[i*s_item+j];
    }
    ta_cnt[newproc[i]]++;
  }

  if (!send) {
    ierr = PetscFree(*aindata);CHKERRQ(ierr);
    *aindata=t_data;

    ierr = PetscFree(ta_procdata);CHKERRQ(ierr);
    ierr = PetscFree(ta_n_items);CHKERRQ(ierr);
    ierr = PetscFree(ta_cnt);CHKERRQ(ierr);
    ierr = PetscFree(ta_getn_items);CHKERRQ(ierr);

    MagparFunctionInfoReturn(0);
  }

  /* sendrecv the number of items to be sent/received */
  int t_gettot=0;
  for (int i=0; i < insize-1; i++) {
    int pto,pfrom;
    pto   = (rank+i+1+insize) % insize;
    pfrom = (rank-i-1+insize) % insize;

#ifndef UNIPROC
    MPI_Status status;
    MPI_Sendrecv(&ta_n_items[pto],1,MPI_INT,pto,0,
                 &ta_getn_items[pfrom],1,MPI_INT,pfrom,0,PETSC_COMM_WORLD,&status);
#endif

    t_gettot += ta_getn_items[pfrom];
  }

  /* allocate array for final local data */
  /* at least one byte to avoid failure of malloc: Cannot malloc size zero! */
  /* may happen with bndfacvert, that no data are to be delivered */
  int *t_fdata;
  ierr = PetscMalloc(s_item*(t_gettot+ta_n_items[rank])*sizeof(int)+1,&t_fdata);CHKERRQ(ierr);

  ta_procdata[rank]=t_fdata;

  /* copy local data from original array into final array */
  for (int i=0; i<n_items; i++)
  {
    /* first processor copies its data to beginning of the array */
    assert(rank==0);
    if (newproc[i]==rank) {
      for (int j=0; j<s_item; j++)
      {
        *ta_procdata[rank]=indata[i*s_item+j];
        ta_procdata[rank]++;
      }
    }
  }

  /* check counter for consistency */
  assert(ta_procdata[rank]-t_fdata==s_item*ta_n_items[rank]);

  /* sendrecv data to/from all(!) processors */
  for (int i=0; i < insize-1; i++) {
    int pto,pfrom;
    pto   = (rank+i+1+insize) % insize;
    pfrom = (rank-i-1+insize) % insize;

#ifndef UNIPROC
    MPI_Status status;
    MPI_Sendrecv(ta_procdata[pto],s_item*ta_n_items[pto],MPI_INT,pto,1,
                 ta_procdata[rank],s_item*ta_getn_items[pfrom],MPI_INT,pfrom,1,PETSC_COMM_WORLD,&status);
#endif

    ta_procdata[rank] += s_item*ta_getn_items[pfrom];
  }

  /* check counter for consistency */
  assert(ta_procdata[rank]-t_fdata==s_item*(t_gettot+ta_n_items[rank]));

  ierr = PetscFree(t_data);CHKERRQ(ierr);
  if (n_items>0) {
    ierr = PetscFree(*aindata);CHKERRQ(ierr);
  }
  *aindata=t_fdata;

  t_gettot += ta_n_items[rank];

  ierr = PetscFree(ta_procdata);CHKERRQ(ierr);
  ierr = PetscFree(ta_n_items);CHKERRQ(ierr);
  ierr = PetscFree(ta_cnt);CHKERRQ(ierr);
  ierr = PetscFree(ta_getn_items);CHKERRQ(ierr);

  MagparFunctionInfoReturn(t_gettot);
}