reorder.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2006-2010 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: reorder.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "init.h"
#include "util/util.h"

int CheckPartition(GridData *gdata,int parts)
{
  MagparFunctionLogBegin;

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

  if (rank) {
    MagparFunctionLogReturn(0);
  }

  int *nvertloc,*neleloc;
  ierr = PetscMalloc((parts+1)*sizeof(int),&nvertloc);CHKERRQ(ierr);
  ierr = PetscMalloc((parts+1)*sizeof(int),&neleloc);CHKERRQ(ierr);

  for (int i=0;i<parts+1;i++) {
    nvertloc[i]=0;
    neleloc[i]=0;
  }

  /* calculate offset in global arrays for each processor */

  for (int i=0;i<gdata->n_vert;i++) {
    assert(gdata->vertnewproc[i]<parts);
    nvertloc[gdata->vertnewproc[i]+1]++;
  }
  for (int i=0;i<gdata->n_ele;i++) {
    assert(gdata->elenewproc[i]<parts);
    neleloc[gdata->elenewproc[i]+1]++;
  }
  for (int i=1;i<parts+1;i++) {
    nvertloc[i] += nvertloc[i-1];
    neleloc[i]  += neleloc[i-1];
  }

  PetscPrintf(PETSC_COMM_WORLD,"mesh partitioning:\n");
  for (int i=0;i<parts;i++) {
    PetscPrintf(PETSC_COMM_WORLD,
      "proc. %2i: %7i elements %7i vertices\n",
      i,
      neleloc[i+1]-neleloc[i],
      nvertloc[i+1]-nvertloc[i]
    );
  }

  assert(nvertloc[0]==0);
  assert(nvertloc[parts]==gdata->n_vert);
  assert(neleloc[0]==0);
  assert(neleloc[parts]==gdata->n_ele);

  ierr = PetscFree(nvertloc);CHKERRQ(ierr);
  ierr = PetscFree(neleloc);CHKERRQ(ierr);

  /* display matrix structure (if X11 compiled in), calculate bandwidth */
  /* TODO: only if X11 available for visualization: quite slow
     why? matviewstruct is fast, must be MatPermute then
  */

  IS isvert;
  ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_vert,gdata->vertl2g,&isvert);CHKERRQ(ierr);
  ierr = ISSetPermutation(isvert);CHKERRQ(ierr);

  Mat Atest,Btest;
  matcreateseqadj(gdata->n_vert,gdata->n_ele,gdata->elevert,&Atest);
  ierr = MatPermute(Atest,isvert,isvert,&Btest);CHKERRQ(ierr);

  ierr = ISDestroy(isvert);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"Adjacency matrix before partitioning:\n");
  ierr = matviewstruct(PETSC_COMM_SELF,Atest);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"Adjacency matrix after partitioning:\n");
  ierr = matviewstruct(PETSC_COMM_SELF,Btest);CHKERRQ(ierr);

  ierr = MatDestroy(Btest);CHKERRQ(ierr);
  ierr = MatDestroy(Atest);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


/* Permute data in gdata structure with new numbering */
int PermuteData(GridData *gdata,IS isvert,IS isele)
{
  MagparFunctionLogBegin;

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

  if (rank) {
    MagparFunctionLogReturn(0);
  }

  AO t_ao;

  /* permute vertex data */
  ierr = AOCreateBasicIS(isvert,PETSC_NULL,&t_ao);CHKERRQ(ierr);

  ierr = AOApplicationToPetscPermuteReal(t_ao,ND,gdata->vertxyz);CHKERRQ(ierr);
  ierr = AOApplicationToPetsc(t_ao,NV*gdata->n_ele,(PetscInt*)gdata->elevert);CHKERRQ(ierr);

#ifdef ADDONS
  if (gdata->vertcoupl!=PETSC_NULL) {
    ierr = AOApplicationToPetscPermuteInt(t_ao,1,gdata->vertcoupl);CHKERRQ(ierr);
    ierr = AOApplicationToPetsc(t_ao,gdata->n_vert,(PetscInt*)gdata->vertcoupl);CHKERRQ(ierr);
  }
#endif

/* more elegant than implementation with AOApplicationToPetscPermuteReal
   below but crashes with SEGV on single proc!
  ierr = VecPermute(gdata->M,isvert,PETSC_FALSE);CHKERRQ(ierr);
#ifdef ADDONS
  if (gdata->VH1!=PETSC_NULL) {
    ierr = VecPermute(gdata->VH1,isvert,PETSC_FALSE);CHKERRQ(ierr);
  }
  if (gdata->VH2!=PETSC_NULL) {
    ierr = VecPermute(gdata->VH2,isvert,PETSC_FALSE);CHKERRQ(ierr);
  }
#endif
*/

  PetscReal *ta_data;

  ierr = VecGetArray(gdata->M,&ta_data);CHKERRQ(ierr);
  ierr = AOApplicationToPetscPermuteReal(t_ao,ND,ta_data);CHKERRQ(ierr);
  ierr = VecRestoreArray(gdata->M,&ta_data);CHKERRQ(ierr);

#ifdef ADDONS
  if (gdata->VH1!=PETSC_NULL) {
    ierr = VecGetArray(gdata->VH1,&ta_data);CHKERRQ(ierr);
    ierr = AOApplicationToPetscPermuteReal(t_ao,ND,ta_data);CHKERRQ(ierr);
    ierr = VecRestoreArray(gdata->VH1,&ta_data);CHKERRQ(ierr);
  }
  if (gdata->VH2!=PETSC_NULL) {
    ierr = VecGetArray(gdata->VH2,&ta_data);CHKERRQ(ierr);
    ierr = AOApplicationToPetscPermuteReal(t_ao,ND,ta_data);CHKERRQ(ierr);
    ierr = VecRestoreArray(gdata->VH2,&ta_data);CHKERRQ(ierr);
  }

  /* renumber vertcoupl such that all are coupled to the vertex with min. id */
  if (gdata->vertcoupl!=PETSC_NULL) {
    for (int i=0;i<gdata->n_vert;i++) {
      if (gdata->vertcoupl[i]<0) continue;

      if (gdata->vertcoupl[i]>i) {
        int k;
        k=gdata->vertcoupl[i];
        for (int j=i+1;j<gdata->n_vert;j++) {
          if (gdata->vertcoupl[j]<0) continue;
          if (gdata->vertcoupl[j]==k) {
            gdata->vertcoupl[j]=i;
          }
          else if (gdata->vertcoupl[gdata->vertcoupl[j]]==i) {
            gdata->vertcoupl[j]=i;
          }
          else if (gdata->vertcoupl[gdata->vertcoupl[j]]==k) {
            gdata->vertcoupl[j]=i;
          }
        }
      }
      assert(gdata->vertcoupl[i]<i || gdata->vertcoupl[gdata->vertcoupl[i]]==i);
      if (gdata->vertcoupl[i]<i) {
        assert(gdata->vertcoupl[gdata->vertcoupl[gdata->vertcoupl[i]]]==gdata->vertcoupl[i]);
      }
    }
  }
#endif

  ierr = AOPetscToApplicationPermuteInt(t_ao,1,gdata->vertl2g);CHKERRQ(ierr);
  ascat(&gdata->vertnewproc,gdata->n_vert,1,gdata->vertnewproc,size,0);
  ierr = AODestroy(t_ao);CHKERRQ(ierr);


  /* permute element data */
  ierr = AOCreateBasicIS(isele,PETSC_NULL,&t_ao);CHKERRQ(ierr);
  ierr = AOApplicationToPetscPermuteInt(t_ao,NV,gdata->elevert);CHKERRQ(ierr);
  ierr = AOApplicationToPetscPermuteInt(t_ao,1,gdata->eleprop);CHKERRQ(ierr);

  ierr = AOPetscToApplicationPermuteInt(t_ao,1,gdata->elel2g);CHKERRQ(ierr);
  ascat(&gdata->elenewproc,gdata->n_ele,1,gdata->elenewproc,size,0);
  ierr = AODestroy(t_ao);CHKERRQ(ierr);


  /* newproc, l2g should be in ascending order now */
  for (int i=1; i<gdata->n_vert; i++)
  {
    assert(gdata->vertnewproc[i]>=gdata->vertnewproc[i-1]);
    assert(gdata->vertl2g[i]==i);
  }
  for (int i=1; i<gdata->n_ele; i++)
  {
    assert(gdata->elenewproc[i]>=gdata->elenewproc[i-1]);
    assert(gdata->elel2g[i]==i);
  }

  MagparFunctionLogReturn(0);
}


/* Reorder (Renumber) nodes to reduce matrix bandwidth */
int OptimizeBandwidth(GridData *gdata)
{
  PetscInt     dooptimizebw;
  PetscTruth   flg;

  MagparFunctionLogBegin;

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

  if (rank) {
    MagparFunctionLogReturn(0);
  }

  ierr = PetscOptionsGetInt(PETSC_NULL,"-optimizebw",&dooptimizebw,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    dooptimizebw=1;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -optimizebw not found, using default value %i (enabled)!\n",
      dooptimizebw
    );
  }
  if (dooptimizebw<0) {
    SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"-optimizebw invalid! Exiting!\n");
  }
  if (dooptimizebw==0) {
    PetscPrintf(PETSC_COMM_WORLD,"Matrix bandwidth optimization not performed.\n");
    MagparFunctionLogReturn(0);
  }

  Mat Atest;
  IS isvert,isrow;

  matcreateseqadj(gdata->n_vert,gdata->n_ele,gdata->elevert,&Atest);

  ierr = MatGetOrdering(Atest,MATORDERING_RCM,&isrow,&isvert);CHKERRQ(ierr);
  ierr = ISDestroy(isrow);CHKERRQ(ierr);
  ierr = MatDestroy(Atest);CHKERRQ(ierr);

  int iss;
  ierr = ISGetSize(isvert,&iss);CHKERRQ(ierr);
  assert(iss==gdata->n_vert);

#if PETSC_VERSION >= 300
  const PetscInt *pind;
#else
  PetscInt *pind;
#endif
  ierr = ISGetIndices(isvert,&pind);CHKERRQ(ierr);
  for (int i=0; i<gdata->n_vert; i++) {
    gdata->vertl2g[i]=pind[i];
  }
  ierr = ISRestoreIndices(isvert,&pind);CHKERRQ(ierr);

  for (int i=0; i<gdata->n_ele; i++) {
    gdata->elel2g[i]=i;
  }

  /* generate trivial mesh partitioning (just split according to vertex id) */
  for (int i=0;i<gdata->n_vert;i++) {
    gdata->vertnewproc[i]=distint(gdata->n_vert,size,i);
  }

  /* partition elements based on vertex data */
  for (int i=0;i<gdata->n_ele;i++) {
    /* trivial partitioning */
    gdata->elenewproc[i]=distint(gdata->n_ele,size,i);

    /* partition elements based on vertex data */
/* disabled to make sure that nothing gets renumbered
    gdata->elenewproc[i]=gdata->vertnewproc[gdata->elevert[NV*i+0]];
*/
  }

  /* Create a global reordering (mapping) of the elements
     by permuting the trivial mapping
  */
  ascat(&gdata->elel2g,gdata->n_ele,1,gdata->elenewproc,size,0);

  IS isele;
  ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_ele,gdata->elel2g,&isele);CHKERRQ(ierr);
  ierr = ISSetPermutation(isele);CHKERRQ(ierr);

  ierr = CheckPartition(gdata,size);CHKERRQ(ierr);
  ierr = PermuteData(gdata,isvert,isele);CHKERRQ(ierr);

  ierr = ISDestroy(isvert);CHKERRQ(ierr);
  ierr = ISDestroy(isele);CHKERRQ(ierr);


  MagparFunctionLogReturn(0);
}


int TrivialPartitioning(GridData *gdata)
{
  MagparFunctionLogBegin;

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

  if (rank) {
    MagparFunctionLogReturn(0);
  }

#if 1
  /* create trivial local->global mappings */
  for (int i=0; i<gdata->n_vert; i++) {
    gdata->vertl2g[i]=i;
  }
  for (int i=0; i<gdata->n_ele; i++) {
    gdata->elel2g[i]=i;
  }
#else
  /* DEBUG: revert numbering of nodes (for testing purposes)
     (activate PermuteDate below!)
   */
  for (int i=0; i<gdata->n_vert; i++) {
    gdata->vertl2g[i]=gdata->n_vert-1-i;
  }
  for (int i=0; i<gdata->n_ele; i++) {
    gdata->elel2g[i]=gdata->n_ele-1-i;
  }
#endif

  /* generate trivial mesh partitioning (just split according to vertex id) */
  for (int i=0;i<gdata->n_vert;i++) {
    gdata->vertnewproc[i]=distint(gdata->n_vert,size,i);
  }

  for (int i=0;i<gdata->n_ele;i++) {
    /* trivial partitioning */
    gdata->elenewproc[i]=distint(gdata->n_ele,size,i);

    /* partition elements based on vertex data */
/* disabled to make sure that nothing gets renumbered
    gdata->elenewproc[i]=gdata->vertnewproc[gdata->elevert[NV*i+0]];
*/
  }

  /* Create a global reordering (mapping) of the elements
     by permuting the trivial mapping
  */
  ascat(&gdata->elel2g,gdata->n_ele,1,gdata->elenewproc,size,0);

  /* Permute data according to new numbering */
  IS isvert,isele;
  ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_vert,gdata->vertl2g,&isvert);CHKERRQ(ierr);
  ierr = ISSetPermutation(isvert);CHKERRQ(ierr);
  ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_ele,gdata->elel2g,&isele);CHKERRQ(ierr);
  ierr = ISSetPermutation(isele);CHKERRQ(ierr);

  ierr = CheckPartition(gdata,size);CHKERRQ(ierr);
#if 0
  /* Apparently we don't have to do this for TrivialPartitioning (no renumbering) */
  ierr = PermuteData(gdata,isvert,isele);CHKERRQ(ierr);
#endif

  ierr = ISDestroy(isvert);CHKERRQ(ierr);
  ierr = ISDestroy(isele);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


/* Reorder/renumber nodes and elements */
int Reorder(GridData *gdata)
{
  MagparFunctionLogBegin;

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

  if (rank) {
    MagparFunctionLogReturn(0);
  }

  ierr = PetscMalloc(gdata->n_vert*sizeof(int),&gdata->vertl2g);CHKERRQ(ierr);
  ierr = PetscMalloc(gdata->n_ele*sizeof(int),&gdata->elel2g);CHKERRQ(ierr);
  ierr = PetscMalloc(gdata->n_vert*sizeof(int),&gdata->vertnewproc);CHKERRQ(ierr);
  ierr = PetscMalloc(gdata->n_ele*sizeof(int),&gdata->elenewproc);CHKERRQ(ierr);

  ierr = TrivialPartitioning(gdata);CHKERRQ(ierr);

  /* Optimize stiffness matrix bandwidth */
  ierr = OptimizeBandwidth(gdata);CHKERRQ(ierr);

  /* mesh partitioning */
  ierr = DataPartitionElementsSer(gdata);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}